{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 06. Force Matching \n",
    "\n",
    "In this tutorial, we show how to use force matching to fit a potential energy function to forces from a simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# disable GPU. Remove this if you've compiled HOOMD for GPU\n",
    "import os\n",
    "os.environ['CUDA_VISIBLE_DEVICES'] = '-1'\n",
    "\n",
    "\n",
    "import hoomd, hoomd.htf as htf, hoomd.md\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import tensorflow as tf\n",
    "import math"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Building the model graph\n",
    "\n",
    "We build an LJ potential with unknown, trainable parameters (`epsilon`, `sigma`) which start out at 0.9 and 1.1. We then obtain forces from our potential and the simulation. Force matching is used to modify the LJ potential until the forces agree. We use Keras layers to implement the trainable parameters. We'll make our starting parameters $\\epsilon=1.2$ and $\\sigma=0.9$. The goal is train them to reach $\\sigma, \\epsilon = 1.0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "class LJLayer(tf.keras.layers.Layer):\n",
    "    def __init__(self, sig, eps):\n",
    "        super().__init__(self, name='lj')\n",
    "        self.start = [sig, eps]\n",
    "        self.w = self.add_weight(\n",
    "            shape=[2],\n",
    "            initializer=tf.constant_initializer([sig, eps]),\n",
    "            constraint=tf.keras.constraints.NonNeg())\n",
    "\n",
    "    def call(self, r):\n",
    "        r6 = tf.math.divide_no_nan(self.w[1]**6, r**6)\n",
    "        energy = self.w[0] * 4.0 * (r6**2 - r6)\n",
    "        # divide by 2 to remove double count\n",
    "        return energy / 2.\n",
    "    \n",
    "class TrainableLJ(htf.SimModel):\n",
    "    def setup(self):\n",
    "        self.lj = LJLayer(0.9, 1.2)\n",
    "\n",
    "    def compute(self, nlist):\n",
    "        # get r\n",
    "        r = htf.safe_norm(tensor=nlist[:, :, :3], axis=2)\n",
    "        p_energy = self.lj(r)\n",
    "        energy = tf.reduce_sum(input_tensor=p_energy, axis=1)\n",
    "        forces = htf.compute_nlist_forces(nlist, energy)\n",
    "        return forces, self.lj.w, energy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running the simulation\n",
    "\n",
    "The simulation consists of LJ particles whose sigma and epsilon values are 1.0. Note that when we compile our Keras model, we prevent the loss from looking at our other outputs by using `None`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HOOMD-blue 2.5.2 DOUBLE HPMC_MIXED TBB SSE SSE2 SSE3 \n",
      "Compiled: 04/30/2019\n",
      "Copyright (c) 2009-2019 The Regents of the University of Michigan.\n",
      "-----\n",
      "You are using HOOMD-blue. Please cite the following:\n",
      "* J A Anderson, C D Lorenz, and A Travesset. \"General purpose molecular dynamics\n",
      "  simulations fully implemented on graphics processing units\", Journal of\n",
      "  Computational Physics 227 (2008) 5342--5359\n",
      "* J Glaser, T D Nguyen, J A Anderson, P Liu, F Spiga, J A Millan, D C Morse, and\n",
      "  S C Glotzer. \"Strong scaling of general-purpose molecular dynamics simulations\n",
      "  on GPUs\", Computer Physics Communications 192 (2015) 97--107\n",
      "-----\n",
      "HOOMD-blue is running on the CPU\n",
      "notice(2): Group \"all\" created containing 256 particles\n",
      "notice(2): Force mode is FORCE_MODE.hoomd2tf \n",
      "notice(2): Starting TensorflowCompute \n",
      "notice(2): completed reallocate\n",
      "notice(2): -- Neighborlist exclusion statistics -- :\n",
      "notice(2): Particles with 0 exclusions             : 256\n",
      "notice(2): Neighbors included by diameter          : no\n",
      "notice(2): Neighbors excluded when in the same body: no\n",
      "** starting run **\n",
      "Time 00:00:08 | Step 1000 / 1000 | TPS 115.593 | ETA 00:00:00\n",
      "Average TPS: 115.585\n",
      "---------\n",
      "-- Neighborlist stats:\n",
      "128 normal updates / 10 forced updates / 0 dangerous updates\n",
      "n_neigh_min: 5 / n_neigh_max: 17 / n_neigh_avg: 11.4609\n",
      "shortest rebuild period: 7\n",
      "-- Cell list stats:\n",
      "Dimension: 8, 8, 1\n",
      "n_min    : 1 / n_max: 7 / n_avg: 4\n",
      "** run complete **\n"
     ]
    }
   ],
   "source": [
    "# run the simulation\n",
    "N = 256\n",
    "hoomd.context.initialize('--mode=cpu')\n",
    "\n",
    "model = TrainableLJ(64, output_forces = False)\n",
    "model.compile('Adam', loss=['MeanSquaredError', None, None])\n",
    "tfcompute = htf.tfcompute(model)\n",
    "rcut = 7.5\n",
    "system = hoomd.init.create_lattice(\n",
    "    unitcell=hoomd.lattice.sq(a=4.0),\n",
    "    n=[int(math.sqrt(N)), int(math.sqrt(N))])\n",
    "nlist = hoomd.md.nlist.cell(check_period=1)\n",
    "lj = hoomd.md.pair.lj(r_cut=rcut, nlist=nlist)\n",
    "lj.pair_coeff.set('A', 'A', epsilon=1, sigma=1)\n",
    "hoomd.md.integrate.mode_standard(dt=0.005)\n",
    "hoomd.md.integrate.nve(group=hoomd.group.all(\n",
    "        )).randomize_velocities(kT=2, seed=2)\n",
    "tfcompute.attach(nlist, train=True, r_cut=rcut, save_output_period=5)\n",
    "hoomd.run(1e3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting training progress"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAvvklEQVR4nO3deXxU1f3/8deZZLIvZINsQBL2ECBAWBQQKi6ogGJdK+5LtbW2tlptv/5Uqlbtl2+tttZqrbt1ww0V3BBFUYGw72uAhCX7vs/M+f1xhgAhkBCS3OTm83w88iBz7525nzs3vOfOueeeq7TWCCGE6PocVhcghBCibUigCyGETUigCyGETUigCyGETUigCyGETfhateLo6GidlJRk1eqFEKJLWrlyZYHWOqapeZYFelJSEpmZmVatXgghuiSl1J7jzZMmFyGEsAkJdCGEsAkJdCGEsIlm29CVUi8A04E8rXVaE/OvAu4BFFAO3Ka1XtvWhQoh7Ke+vp6cnBxqamqsLqXTCQgIIDExEafT2eLntOSk6EvAP4BXjjM/C5istS5WSp0HPAeMa3EFQohuKycnh9DQUJKSklBKWV1Op6G1prCwkJycHJKTk1v8vGabXLTWS4CiE8z/Xmtd7H34I5DY4rULIbq1mpoaoqKiJMwbUUoRFRV10t9c2roN/UZgYRu/phDCxiTMm9aa96XNAl0p9RNMoN9zgmVuUUplKqUy8/PzW7Weg6U1zPloI/VuTysrFUIIe2qTQFdKDQeeBy7UWhcebzmt9XNa6wytdUZMTJMXOjVrTXYJLy7dzZNfbm9ltUII0bbmz5/PY489BsCDDz7I3LlzLanjlK8UVUr1Ad4DrtZabzv1kk5sWlosl2Uk8s+vdzB5UAxjkiLbe5VCCHFCM2fOZObMmVaX0fwRulLqDeAHYJBSKkcpdaNS6lal1K3eRe4HooB/KqXWKKXa/Xr++2cMJTEiiDvfWkN5TX17r04IYXOvvfYaY8eOJT09nZ///Oe43W5CQkK48847GTp0KFOnTuVQM/FTTz1Famoqw4cP54orrgDgpZde4vbbbz/mddesWcP48eMZPnw4s2bNorjY9B+ZMmUK99xzD2PHjmXgwIF8++23bbIdzR6ha62vbGb+TcBNbVJNC4X4+/LE5elc+q/veWD+Rv56WXpHrl4I0Q7mfLSRTfvL2vQ1U+PDeGDG0BMus3nzZt566y2WLl2K0+nkF7/4Ba+//jqVlZVkZGTwxBNP8Kc//Yk5c+bwj3/8g8cee4ysrCz8/f0pKSk54Wtfc801/P3vf2fy5Mncf//9zJkzh7/97W8AuFwuli9fzoIFC5gzZw5ffvnlKW9vl71SdHTfCG4/cwDvrdrHx+v2W12OEKKLWrRoEStXrmTMmDGkp6ezaNEidu3ahcPh4PLLLwdg9uzZfPfddwAMHz6cq666itdeew1f3+MfE5eWllJSUsLkyZMBuPbaa1myZEnD/IsvvhiA0aNHs3v37jbZFstGW2wLvzqzP0u25fM/729gdN8I4sIDrS5JCNFKzR1JtxetNddeey2PPvroUdMfeuihox4f6kb4ySefsGTJEj766CMeeeQR1q9f36r1+vv7A+Dj44PL5WrVazTWZY/QAZw+Dp64PJ16t4ffvb0Wj0dbXZIQoouZOnUq8+bNIy8vD4CioiL27NmDx+Nh3rx5APz3v/9l4sSJeDwesrOz+clPfsLjjz9OaWkpFRUVTb5ueHg4ERERDe3jr776asPRenvp0kfoAMnRwdw/PZV731vPf77L4uYzUqwuSQjRhaSmpvLwww9zzjnn4PF4cDqdPP300wQHB7N8+XIefvhhevbsyVtvvYXb7Wb27NmUlpaiteaOO+6gR48ex33tl19+mVtvvZWqqipSUlJ48cUX23VblNbWHNVmZGTotrrBhdaan7+6kq+35jP/VxMYHBvWJq8rhGhfmzdvZsiQIVaX0aSQkJDjHn13lKbeH6XUSq11RlPLd+kml0OUUjx68TBCA3y55931uKXpRQjRDdki0AGiQvy5f0Yqa7NLeOn73VaXI4To4qw+Om8N2wQ6wMwR8fxkUAxzP9tKdlGV1eUIIUSHslWgK6V4eNYwHAr+54MNWHV+QAghrGCrQAdI6BHI76cNZsm2fOavlQuOhBDdh+0CHWD2+L4MSwjn0QVbqKprmw77QgjR2dky0H0cigdmpHKwrIZnvt5pdTlCiC7mpptuYtOmTVaXcdJsGegAGUmRzBwRz3NLdnGgtNrqcoQQXcjzzz9Pamqq1WWcNNsGOsDd5w7CozVPLZKbYQghmlZZWckFF1zAiBEjSEtL46233mLKlCkcuvDxP//5DwMHDmTs2LHcfPPNDcPkXnfdddx2222MHz+elJQUvv76a2644QaGDBnCdddd1/D6t912GxkZGQwdOpQHHnigXbely1/6fyK9I4O4alxfXv1xDzdNSqFfTIjVJQkhjmfhvXCwdQNdHVfsMDjvsRMu8umnnxIfH88nn3wCmFESn3nmGQD279/PQw89xKpVqwgNDeXMM89kxIgRDc8tLi7mhx9+YP78+cycOZOlS5fy/PPPM2bMGNasWUN6ejqPPPIIkZGRuN1upk6dyrp16xg+fHjbbqeXrY/QAW4/sz/+vg7+LkfpQogmDBs2jC+++IJ77rmHb7/9lvDw8IZ5y5cvZ/LkyURGRuJ0Orn00kuPeu6MGTNQSjFs2DB69erFsGHDcDgcDB06tGFI3LfffptRo0YxcuRINm7c2K5t87Y+QgeIDvFn9vi+PP/tLn591kCSo4OtLkkI0ZRmjqTby8CBA1m1ahULFizgvvvuY+rUqS1+7qEhcB0OR8Pvhx67XC6ysrKYO3cuK1asICIiguuuu46ampo234aG9bbbK3ciN09Kwenj4OnFO6wuRQjRyezfv5+goCBmz57N3XffzapVqxrmjRkzhm+++Ybi4mJcLhfvvvvuSb12WVkZwcHBhIeHk5uby8KFC9u6/KPY/ggdICbUn5+N68OrP+zhrnMGERseYHVJQohOYv369dx99904HA6cTifPPPMMd911FwAJCQn88Y9/ZOzYsURGRjJ48OCjmmSaM2LECEaOHMngwYPp3bs3EyZMaK/NAGwyfG5LZBdVMfl/F3PblH7cfe7gDluvEOL4OvPwuYdUVFQQEhKCy+Vi1qxZ3HDDDcyaNatD1t0th89tid6RQZyd2ov/LttLdZ3b6nKEEF3Egw8+SHp6OmlpaSQnJ3PRRRdZXdJxdYsml0Oun5DMZxtz+WDNPq4c28fqcoQQXcDcuXOtLqHFus0ROsC45EhS48J4cWmWjMQoRCch/xeb1pr3pVsFulKK6ycksS23gqU7Cq0uR4huLyAggMLCQgn1RrTWFBYWEhBwch04ulWTC8CMEfE8/ukWXliaxcQB0VaXI0S3lpiYSE5ODvn5+VaX0ukEBASQmJh4Us/pdoEe4PThZ+P68tSi7WQXVdE7MsjqkoTotpxOJ8nJyVaXYRvdqsnlkMvH9EYpeCcz2+pShBCizXTLQE/oEcikATG8szIHt0fa7oQQ9tAtAx3g8ozeHCitYcl2absTQthDtw30s1J7Ehnsx9srpNlFCGEP3TbQ/X19uHhkAl9syqWgotbqcoQQ4pR120AHc3LU5dG8v2qf1aUIIcQp69aBPqBXKKP69ODNFXvlwgYhRJfXrQMd4LKM3uzMr2RtTqnVpQghxCnp9oF+3rA4/HwdfLBaml2EEF1bs4GulHpBKZWnlNpwnPmDlVI/KKVqlVJ3tX2J7Ss80MlZQ3ry0dr91Ls9VpcjhBCt1pIj9JeAaSeYXwTcAXSdMSYbuSg9gcLKOr7bXmB1KUII0WrNBrrWegkmtI83P09rvQKob8vCOtKUQT3pEeTkfWl2EUJ0YR3ahq6UukUplamUyuxMo6v5+TqYPjyOzzcdpKLWZXU5QgjRKh0a6Frr57TWGVrrjJiYmI5cdbNmjUygpt7DZxsOWl2KEEK0Srfv5XLIqD4R9I4MlGYXIUSXJYHupZRiVnoCS3cWkFtWY3U5Qghx0lrSbfEN4AdgkFIqRyl1o1LqVqXUrd75sUqpHOC3wH3eZcLat+z2MTM9Aa3h43UHrC5FCCFOWrN3LNJaX9nM/IPAyd0nqZPq3zOEtIQw5q/Zx40T5S4qQoiuRZpcGrlwRAJrc0rJKqi0uhQhhDgpEuiNTB8Rh1Iwf81+q0sRQoiTIoHeSFx4IOOSI/lw7T4ZgVEI0aVIoDfhwvQEduVXsnF/mdWlCCFEi0mgN+G8tFicPooP10ifdCFE1yGB3oQeQX5MHtiT+Wv34/ZIs4sQomuQQD+OC9PjyS2rZXnWccclE0KITkUC/TjOGtKLID8f5q+VZhchRNcggX4cgX4+nDs0lgXrD1LrcltdjhBCNEsC/QRmpsdTWl3Pkm1y4wshROcngX4CE/tHExnsxzuZ2VaXIoQQzZJAPwGnj4Orx/fl8025rMkusbocIYQ4IQn0Ztx8RgrRIX48umCzXDkqhOjUJNCbEeLvyx1TB7Asq4iPZFhdIUQnJoHeAj8b24eRfXpw3/vrOVgqN78QQnROEugt4Ovj4K+XpVPv1tz/4QaryxFCiCZJoLdQcnQwt5yRwuebcmWsdCFEpySBfhKuGt8Hp4/i5e93W12KEEIcQwL9JPQMDWDG8HjeycymrKbe6nKEEOIoEugn6boJSVTWuflgtYzxIoToXCTQT9LwxB4MjQ/jrRVy9agQonORQG+FK8b0ZuP+MjbsK7W6FCGEaCCB3goz0xPw93Xw5oq9VpcihBANJNBbITzQyXlpsXy09gB1Lo/V5QghBCCB3mqHhtb9bke+1aUIIQQggd5qE/vHEB7o5KO1Mr6LEKJzkEBvJT9fB+elxfL5xoPU1MsdjYQQ1pNAPwUzRsRTWefmy825DdPWZJfw2MItvL86R4JeCNGhfK0uoCsbnxJFfHgAb2fmMH14PJv2lzH7+WVU1LoAWLj+IM9ePRqlFDX1brKLqugTFUSty0NhRR29IwLx9em4z9SaejcBTp8OW58QomNJoJ8CH4fiktGJ/H3xDlbuKeKXr68mxN+XT38ziY/XHeCxhVv4fx9uoKC8jm+25VNd78ahwOO9T0aA08GUgT25aVIyqfFh1Ls0uwoqyCqoxOXWxIT6Ex3iz+YDZXy8/gBBTh/69QzmzMG9SO/dAx+HalGdHo/mzws28/x3WQyND+OGCclcPCoBpY7//FqXm2e/2UVlrYthieFcMCzuuMtrrSmuqifIz+eoD4yaejcr9xSTHB1MfI9AANweTUFFLXUuDwk9AnE02gaPR6MUJ6xNCNE0ZdVdeDIyMnRmZqYl625L2UVVTPrLYvx8HPj7OnjnttMYHBuG1pqbX8nky8159Az155yhvRiR2IO9RVUE+vkQHezPxv2lvL96H2U1rmbXkxQVhNPHYcLeo4kO8eOMATGM7NOD9N4RDI4LxdnE0X5eWQ1zPtrEJ+sPcMGwOLIKKtl0oIyzhvTi8jG9GZsUSXiQ86jnaK25e9465q3Mwc/HQZ3bw9Xj+/LAjNSjvlHU1LuZv2Y/L32/m00HygAY2CuEc1Jjqah18fG6AxRU1DZMn5YWxwer97G3qAqAyGA/Jg2IZsqgGDL6RpK5p4gH528iwOlg0oAY7j53EL3CAiirqeeF77L4YlMuOcXVXDwqgVvOSCEu3HxIFFbU8tWWPPr1DGFYQniT70Nj5TX1hPj7nvCDY39JNXUuD0nRwc2+nhAdRSm1Umud0eQ8CfRTN/v5Zfy4q5CXbxjLhP7RDdNr6t3syq9kcGzoMUeih1TUuli8JY+9RVU4fRQp0SEkRQfj7+ugoKKWvPJaokP8GNUnAqUUpdX1fLMtny825fLDzgIKKuoAcPoo+kQGkRwdQp/IIAL9HOzMq2TJ9nzq3R5+e/Ygbp2cgtbwwtIs/u/zbVTXu1EKBseGMTYpgjHJkfSLCeHl73fz5opsfj11AL+eOoDHP93Cs0t2MTY5kkcuSsOj4f3V+3hzxV5KquoZ1CuUC0fG4/FoFm/NZ+WeYoL8fBiTFMmVY/uwv6SaeStz2HSgjLSEMC4d3Runj4PM3UV8sy2fwsq6hvdjdN8IEnoE8tnGg/j5OkiLD2fLwTKKq+oZmxxJTIg/n208SKDTh1un9COnuJoPVu+j2nu+IiLIycwR8Vw8KpHhieFHBXZNvZu/frGNtzOzKamqJzEikOnD47lkdCL9e4YctV8+33iQ299YTZ3Lw+DYUManRJESE0x4oJMhcWH0jwk57j5tzOPRHCyrafjWFeh3cs1e23PLeXfVPoor6xiTHMkFw+IaXiO7qIq9RVWMT4k66htbvdtDVkElSVHB+PnKqTI7kUBvZ/nlteSX15IaH9ah69Vak1NczZrsEjYdKCMrv5JdBRXsK66mut5NQkQgE/pFc9uUfvSNOvoos6bezdrsEpZlFbEsq5BVe0oaQtHXoZg9vi/3T09tCK33V+dw3/sbqKwzyzgUnJMay7WnJzE+JfKY4PT3dRw1TWvNgdIa4sIDjpru8Wg2HShjTXYJTh/FJaN74+NQ7C6o5PFPt1BQUUvP0ABum9KPtIRwAPYWVnH3vLUsyyoi2M+Hc4fGcs3pSewvqWbB+gN8vimXOpeHPpFBTBoQTVpCOJW1Lt5ckc2OvAouGB7HkNhQMvcU8932Alwezei+EVw0MoEhsaF8uuEgL36/m7SEcKYPi2PRFnOT8Jr6wxeRRYf4cc7QWM4YEE3fqGDKa1xkF1VRWFmL2wMerXG5NVV1LhZsOEB2UTUASkFiRCADeobSv2fIUT9hAcd+U/rXN7t4/NMt+DoUwf6+lFabD6LrJySTX17Li0uzqHV56BsVxO/PHcz5w2Ipq3Fx8yuZLM8qItDpw+VjenPXuYMI8W++hXXrwXL+9c1OvtiUy7CEcC4bk8iM4fFNnuvRWvPmimw+WXeA+B4BnDm4F+ek9mr4m8krr2HZriL8fR2c1i+K0AAnlbUuPlyzn2255fj7mm9i41IiW/StqrFDnQ5ae15Ia011vZtAp0+XauKTQBfNqnd72Li/jC0Hypg4IJrEiKBjlskuquLrrXkE+/syLiWKBG+7uBU8Hs3uwkp6RwYdEwal1fUsXH+ALzfn8uOuooaT1KlxYdx17kDOHNyrYdn88lreW5XDW5nZ7Mo3Ny5xKLh4VCJzZg4l2BuC9W4PJVX1FFXWsS6nhK+35bN4Sx5VdSfuyaQUjEuO5PxhcQQ6fThQWsP2vAp25FWwM7/iqCuNY8MCGsK9V1gAy7IK+XprPtOHx/HgzKFEBvnx465C7p+/kR15FQCcPyyWs1N78dySLDYfKCM1LoyCilqKq+r4zVkDySqo5N1VOcSGBXD9hCTOGBhDdIg/EUF+x5yD+XprHr94fRU+SjF1SE/W5pSSVVBJSnQw109M5qL0eEK9Hzr55bX89u01fLu9gJToYIqr6iiuqicpKogzBsawI6+C73cWNry2r0Mxqm8Eu/IrKaioJcjPB5dbU+f2EBnsx/ThcVyYnsCIxHB8fRxorSmrdlFQWUthRR2FFbUUVNZRVFFHgNNBcVU9r/+4h/JaF4kRgVxzWl+mDY3D4YDcsloWrj/Awg0H6RXmz7iUKH42tg+9I83fdEFFLQ99vIkF6w9Q79YkRwdzXlos5w+LY2h8WEO478qv4PNNuTgUpCWEc1pKFEqphtocDvNh4utQx/1AyCuvISu/kup6NxFBfiREBBIV7HdKHyAS6KLb8ng0+0rM0fGh/9BN0Vqzp7CKbbnlDOwV2qJ28zqXh3U5JeSV1xLi70tCRCA9Q/3xdTjwcSh8HArHCU7wuj2a7KIqduRVsD2vgu155ez0/l5V5yahRyCXZiRyx5kDjmreOXRiOcjPpyFgXW4PLy7dzRebc4kO8WP2+L6c3s80/63cU8TjC7eyfHdRw2s4lDmHERHkR0SwH0WVdezIqyA1LoyXrh9Dz7AAtNZ8vimXpxZtZ+P+Mvx8HUzoF0VseACLNudRWl3P/5ueylXj+uD2aBZsOMg7mdms2F1ETKg/l47uzZRBMVTXuflmWz5LtufTI9CPO88eyKg+Paiud/Pt9gLmr93Pl5tyqXV5cPoowgP9KKmqw+U5fjYpBeelxTIkNowfdhUe9eEB5gNkyqAYyqpdrNxbjNaaUX0iiArxY+mOQmpdbq4Y04eeof4syyrih12FuD2axIhAxqdEkV9ey7fb8zmyhPjwAMICneSW1VBcdfh+CA4F/r4+OH1Uw373cSg82nzwNRbgdHDb5P78+qwBJ/jrOr5TCnSl1AvAdCBPa53WxHwFPAmcD1QB12mtVzVXlAS6EE07dATY+GT1qdqRV87WgxXkl9dQUFFHYWUdxZV1FFXVEeLvy+i+EVxzWt+GD4kj61mTXcL8tfv5Zls+FTUu4sIDeOynwxkSd2wzo9ujT/hB1pTymnq+2pLH5gPlFFfWERXiR1SIP9EhfkQG+xEVbH6PCPaj1uWhzmWO7A9Zm13CttxyPNqcp0iLD6dnWABgTm6/tSKbRVtyKamqZ0K/aG6clMzAXqENzy+qrOPLTbks2HCADfvKiAr2Y/KgGG6amIy/rw+LtuSyaHMedW4P0SF+pESHoJRp9ql1eRpqcns0bq3xeDRaQ/+eIQyOCyXIz5eiyjpyiqvYV1zN2ORIzhkae7K7EDj1QD8DqABeOU6gnw/8ChPo44AntdbjmitKAl0IIU7eiQK92TMRWuslQNEJFrkQE/Zaa/0j0EMpFde6UoUQQrRWW/RnSgCOvH1PjnfaMZRStyilMpVSmfn5MkqhEEK0pQ7toKq1fk5rnaG1zoiJienIVQshhO21RaDvA3of8TjRO00IIUQHaotAnw9co4zxQKnWWgYJF0KIDtbspWNKqTeAKUC0UioHeABwAmit/wUswPRw2YHptnh9exUrhBDi+JoNdK31lc3M18Av26wiIYQQrSKj9gghhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE20KNCVUtOUUluVUjuUUvc2Mb+vUmqRUmqdUuprpVRi25cqhBDiRJoNdKWUD/A0cB6QClyplEpttNhc4BWt9XDgT8CjbV2oEEKIE2vJEfpYYIfWepfWug54E7iw0TKpwFfe3xc3MV8IIUQ7a0mgJwDZRzzO8U470lrgYu/vs4BQpVRU4xdSSt2ilMpUSmXm5+e3pl4jb3PrnyuEEDbVVidF7wImK6VWA5OBfYC78UJa6+e01hla64yYmJjWrWn16/DM6bDrm1OpVwghbKclgb4P6H3E40TvtAZa6/1a64u11iOB//FOK2mrIo+SeiFE9Yd3b4Ly3HZZhRBCdEUtCfQVwAClVLJSyg+4Aph/5AJKqWil1KHX+gPwQtuWeQT/ELj0Zagth3dvBM8xXwSEEKJbajbQtdYu4HbgM2Az8LbWeqNS6k9KqZnexaYAW5VS24BewCPtVK/RKxUumAu7v4WvH2vXVQkhRFehtNaWrDgjI0NnZmae2ou8fxusfQNmvwv9p7ZNYUII0YkppVZqrTOamte1rxS9YC7EDIb3boGyA1ZXI4QQlurage4XDJe9DPXV5iSptKcLIboxX6sLOGUxg+D8v8CHv4Qf/gETfm11RUIIYQ4wd39nfgq2QV0laA9oNwy9GEZf2+ar7PqBDpB+FWz7FBY9BP3OhNhhVlckhOgOtIbKfMjfakK7cAeU7IXSbCjaDbWloBwQkQwBYaB8zGNPfbuUY49AVwqmPwnZp8G7N8MtX4MzwOqqhBB243bBzkWw+SNviG+FmtLD853B0KMP9OgNCaMheTIMONs0D3cAewQ6QHAUXPhPeP2nsOhPMO3PVlckhOjK8rfBsmcgJ9M0lfiFQO4GqKuAgHCIHQ5pl0D0QIgZaP4NSzAHmBaxT6ADDDgLMm6EH/8Jwy6BhFFWVySE6Kzc9bBvJbhqzRG0XzA4g6C2DDZ+AN//HRy+0Hss+AZATQmMuAJSfgIDzgFfP6u34Bj2CnSAsx6ALZ/Ax7+BmxeDw8fqioQQnYnHA9/9FZY+Zdq4jyf1Qjh/LoT07LjaTpH9Aj0gHKY9CvOuhxXPw7ifW12REMIK1SVQtt80gfToY47AD24wTbLbP4NBF0D6lRAYAXVVUF9peqI4nJAyGUJjrd6Ck2a/QAcYOgtWv2Z6vQyZAWHxVlckhOgotRWw+M+m/Vt7Dk93OE3vEt9AmPa4OdizsL27Pdgz0JUyV5E+PR4+/YO5+EgIYS8l2ZCz3Iy6GpFk2rh3LoatC8yJy1HXmiNtraFkjxnQLyga0n8GQZFWV98u7BnoAJEpMOm38PWjsHcZ9BlndUVCiFPl8UBFLqx/B756GNy1R88PjIS0i2HUdZA42pISrWTfQAc4/VeQ+QJ8cT/c8Kntvl4JYXuuWti3CvYshb0/mIOzunIzb/B0mHwPhMZB8W7wC4KYIeDo2iOanAp7B7pfMEz5g+nxsnUBDL7A6oqEEGC6DJbmmKYRHz/TLdA3AOqroGiXCe8935s+4IeOwmOGwPBLoddQMyhf3wmHD9JCWnkHNJuxd6ADjLwavn/KjJs+6Hw5ShfCSnVV5v/jiv9AZd7xl1M+EDcCxt4MfU+H3uPNxYPihOwf6D6+cMbd8MFtsHUhDD7f6oqE6D4q8qAoy/QuqasyXQZz18OAc00PtIAw8LjMiKn11eZbdWgcJGaAf6jV1Xc59g90gGGXwTd/gW8eg0HnyVG6EO2tIg++eRxWvQLuusPTAyPgqnlmfBPR5rpHoPv4wsTfwEe/NkNZJk+yuiIh7GH/alj+b9PzJDDCjCroqoaVr5j28JFXmZOXPn7g44ToQdJ00o66R6ADDL/cfN378RkJdCHawtq3YP6vwBlo+oHnbTHdCVGQMgXO/1+IHmBxkd1L9wl0ZyBk3ABL5kLhTojqZ3VFQnROteWw4T3ze1g89Ew1l8EfGhfJVQdf/ckMXpU0CS59+fBRt8fTrbsNWq37BDrAmJvgu7/B8ufgvMetrkaIzmfdO/DpvVBV0GiGMuMkBYSZ+/d66mHsLXDun01TyiES5pbqXoEeGgtpPzXjvPzkj+YPVAhhbtzw/VOwaA70HgdXvmF6m5TshbxNUFkA1UVmwKvwBEiaCP3Psrpq0Uj3CnSA034B696EVa/C6bdbXY0Q7WP1a7D4UXOCsmcqDJlpbs3o8DEnMCtyTU+U8oPm9+zlJrCHXgyznj081neP3pA0wdptES3W/QI9bgT0nQjLnoXxt8l46aJrqymDA2sguCdE9Td/z2vfgA9vN325ew2F3Uth4d1NPFlBcDSExJpuhENmmCFlpdmky+p+gQ7m6rN3rjUjsw2Qr42ii8rfBv+91IxjAuAXapoRy3LMZfFXzTPjmxwabbBgu1kupCeE9DIjD/p0zwiwq+65NwedD0FRsOplCXTR+dXXwA9/hx1fmSFiYwaZS+O3LjSBfcmLZmyUvT+YO9CfeZ+5J8ChG6UrZboVRiRZuBGiI3TPQPf1gxFXwrJ/QUW+DOwjOq+iXfD6pVC4AxLHQI++pr3bVWOGiT3jbojoa5Ydcbm1tQrLdc9ABxh1DfzwD9PeOOEOq6sR4lhl++GVC80deGa/B/2nWl2R6OS6b6DHDDLds1a9YsZNl/FdhJVctWY42dIc0+skdyOsfAk8brh2PiSMsrpC0QV030AHc5T+4S8hexn0GW91NcLu6qqgYJs5OVmwDYqzTD/v4j1QcbDRwgqGTIfJ90JsmiXliq6newd66kWw8F5zlC6BLtpL1hL4/h+wa/HhkQeVA8ITTZt4/7PMXenDE02/79A40xNFLnwTJ6l7B7p/CAz7qRlkaNpj5rJmIdpCVZH55rf+HdjwLoTGw5ibzYFD9ECITAZff6urFDbTvQMdIP0q01a55WNzN3AhWqO2whyJb/8cdi4yTSkAfiEw8bcw+fdmgDgh2pEEeuIY0z933VsS6OLEakrNfS4Ld5gjcE+9GQMlf7OZ7q4zF/ekTDYDwcWPNLdOO3QZvRDtTAJdKTNW+jd/MaPIhcVZXZHoDDwe+HYuHFhrLqcvzTG/e1xmvsMXHE4z0mBYAoz7OQw4RwJcWKpFga6UmgY8CfgAz2utH2s0vw/wMtDDu8y9WusFbVtqOxp2mbld1oZ5pguj6N60hs/vgx+fNuOjoMy44Kf9EgZOg55DzN15hOhkmg10pZQP8DRwNpADrFBKzddabzpisfuAt7XWzyilUoEFQFI71Ns+ovtDwmjT7CKBbn9Fu2D7F6aPd2AEBPYw/wb0MEfci/9sPtzH3WpOlss1CqKLaMkR+lhgh9Z6F4BS6k3gQuDIQNfAoS4i4cD+tiyyQwy/HBb+HnI3Qa9Uq6sR7aGuEubfYcL6RBxOmPIHOOP3EuaiS2lJoCcA2Uc8zgHGNVrmQeBzpdSvgGCgyRGvlFK3ALcA9OnT52RrbV9DL4ZP/2CO0s+eY3U14mRobU5U1pSZ0QODY8xIgr5+pi3cUw87v4Iv50D+Fph0F4y6GvzDzGBX1cXenxJz4jNpormSWIgupq1Oil4JvKS1/j+l1GnAq0qpNK2158iFtNbPAc8BZGRk6DZad9sIiTFjZax/B6Y+IGNCd0ZamwD2Dz+8fyry4ZM7YfNHxy7v8D18EhPMxTuz5x19p52gyHYtWYiO1JJA3wf0PuJxonfakW4EpgForX9QSgUA0UBeWxTZYYZfDu/eCHu/N0dpovPYtxK+eAB2f2uaRBJGmTvxrH/HdBec8keITzdjolQVQGUh1FeZi3d8nKZr6pCZR9//UgibaUmgrwAGKKWSMUF+BdC4w/ZeYCrwklJqCBAA5LdloR1i0HngDIL18yTQrVKRZ8Y28fWD2OGgPbD0b/DVw2YM+8n3QH21uTnJqpfNuN+T75EmEiFoQaBrrV1KqduBzzBdEl/QWm9USv0JyNRazwd+B/xbKXUn5gTpdVrrztWk0hJ+wSbUN30I5/+vHM11tM0fw7zrD493EpYIteVQW2rOccx48ujhGdz1so+EOEKL2tC9fcoXNJp2/xG/bwLscSfZtEvM2Bs7F8PAc6yuxn48btj9nblEvrbMnLTUbhPcWxeappQzfm/uvLN1gTkqH3iuuctU4x4nEuZCHEWuFG2s/1Qzyt2GeRLobeHQiUyH07Rtv32tuamxb4B5n5WPOXnpF2TuwDP9CfAPNc8deZWVlQvR5UigN+brb06ebXzfjF/tF2R1RV2T1pD5H/jiQagrN9McvuAMhoueMUMXy3srRJuSQG/KsEtg9auw/TNz0k00z+2CvE3mhGVdBSz/N2xbCClToP/Z4K41PU/G3AhR/ayuVghbkkBvStIkCOllertIoJ+Yux6+/T9Y+hTUVx6eHhhh+vNP+I306Reig0igN8XhY4I880Vz5aDcOaZpuRvh/Vvh4DpIvdA0VQX0MBfzpEyW8b+F6GAS6MeTdgks+5fpStf45Fx3CfnqYtj7o2lGGXzB4TvsuF3w/ZOw+FHzPlz+GgyZYW2tQggJ9ONKzDCXim+Yd3SgZ30Lr1wIV/wXBk2zrr72dnA9vDDNtIcDBPeE4ZdBeG9Y/RrkrjcnNi/4PwiOtrRUIYQhgX48SkHaT03bcEW+GevF4zHjZGs3fP2o6R9tx9H4aivgnevN7dOufNNc6LP837DsWTPQVc9UuOQF8/4IIToNCfQTSbsEvnsCNn0AY2+Gje+ZPtT9ppr7Ru5cdPRAT51ZzkrzIRQQDkkTYPgVTXcbLN4NH/zCjF547XxInmSm959qmpoq8k0vFTt+kAnRxSmrrtDPyMjQmZmZlqy7xbSGZyaAqwZu+Ayem2J6b9z0Jfx9FEQkw/WfWF1l8w6sg5eng4+fuaCnNNscfQdHA8oMaOWqNv/WV5t55/8vpF9pdeVCiEaUUiu11hlNzZMj9BNRCqb92bSZP3sGlB+AS18CZwBk3ABfPWQGkoro2zH1uOpMk4dfcMuWd9ebnjqL5pjeJzcsNG3ge3805wZqyszgV84A8A00//qHwYgroUfvZl9eCNG5SKA3J2UKjLzaXGg07jboPcZMH3apCfT178AZd7Xf+rWGdW/D4kegZC+gzRF0aCyExEJoL9M9UPmY7pYOXzNeSmWeOYFbUwL9zoQZTx0O6b6nmR8hhK1IoLfEtEfNWNsjjhg1OKIv9DndhO2k37VNm7LWpnfJ1oWmF0nZfnO3+Ypcc8/T9KvMsLLluVBxEMoPwr5VpqlEu02QazegDg9qlfZTczd6afMWwvYk0FvCPxTG3HTs9OGXwsd3woG1JvABNn4APz5jAnTSXTDAe9K0Is/cnMHjNiM57v7W3IAhLME8t77GNIWU7gUURA+E8ARzh/k+p5kPE7niUghxAhLop2LoLFh4r+mXHZ9ujpbfu9m0U3tc8PbVMOtZ2PEFrH3z8Djf/mGmd0xgD9ObZMsC0y7eKxUm3w0Dp0FITyu3TAjRBUmgn4rACHPJ+7q3YeJvzNCwIb1MLxh3Hfz7TBPqvgGmHX7kbHNHpIi+clm8EKLNSaCfqtHXwfq34d9ToboIbvj08I2Hr/7A3MhhxBVyNaUQot1JoJ+qvqdD1AAo3A7T/2ZOXh4SM9D8CCFEB5BAP1VKwYy/wcEN5mhdCCEsIoHeFpImmh8hhLCQ9IMTQgibkEAXQgibkEAXQgibkEAXQgibkEAXQgibkEAXQgibkEAXQgibkEAXQgibsOwWdEqpfGBPK58eDRS0YTldRXfcbtnm7kG2ueX6aq1jmpphWaCfCqVU5vHuqWdn3XG7ZZu7B9nmtiFNLkIIYRMS6EIIYRNdNdCfs7oAi3TH7ZZt7h5km9tAl2xDF0IIcayueoQuhBCiEQl0IYSwiS4X6EqpaUqprUqpHUqpe62up60opXorpRYrpTYppTYqpX7tnR6plPpCKbXd+2+Ed7pSSj3lfR/WKaVGWbsFraeU8lFKrVZKfex9nKyUWubdtreUUn7e6f7exzu885MsLbyVlFI9lFLzlFJblFKblVKn2X0/K6Xu9P5db1BKvaGUCrDjflZKvaCUylNKbThi2knvW6XUtd7ltyulrm3p+rtUoCulfICngfOAVOBKpVSqtVW1GRfwO611KjAe+KV32+4FFmmtBwCLvI/BvAcDvD+3AM90fMlt5tfA5iMePw48obXuDxQDN3qn3wgUe6c/4V2uK3oS+FRrPRgYgdl22+5npVQCcAeQobVOA3yAK7Dnfn4JmNZo2kntW6VUJPAAMA4YCzxw6EOgWVrrLvMDnAZ8dsTjPwB/sLqudtrWD4Gzga1AnHdaHLDV+/uzwJVHLN+wXFf6ARK9f+RnAh8DCnP1nG/jfQ58Bpzm/d3Xu5yyehtOcnvDgazGddt5PwMJQDYQ6d1vHwPn2nU/A0nAhtbuW+BK4Nkjph+13Il+utQROof/MA7J8U6zFe9XzJHAMqCX1vqAd9ZBoJf3d7u8F38Dfg94vI+jgBKttcv7+Mjtathm7/xS7/JdSTKQD7zobWZ6XikVjI33s9Z6HzAX2AscwOy3ldh7Px/pZPdtq/d5Vwt021NKhQDvAr/RWpcdOU+bj2vb9DNVSk0H8rTWK62upQP5AqOAZ7TWI4FKDn8FB2y5nyOACzEfZvFAMMc2S3QL7b1vu1qg7wN6H/E40TvNFpRSTkyYv661fs87OVcpFeedHwfkeafb4b2YAMxUSu0G3sQ0uzwJ9FBK+XqXOXK7GrbZOz8cKOzIgttADpCjtV7mfTwPE/B23s9nAVla63ytdT3wHmbf23k/H+lk922r93lXC/QVwADv2XE/zImV+RbX1CaUUgr4D7BZa/3XI2bNBw6d5b4W07Z+aPo13jPl44HSI77WdQla6z9orRO11kmYffmV1voqYDFwiXexxtt86L24xLt8lzqS1VofBLKVUoO8k6YCm7DxfsY0tYxXSgV5/84PbbNt93MjJ7tvPwPOUUpFeL/dnOOd1jyrTyC04oTD+cA2YCfwP1bX04bbNRHzVWwdsMb7cz6m7XARsB34Eoj0Lq8wPX52AusxPQgs345T2P4pwMfe31OA5cAO4B3A3zs9wPt4h3d+itV1t3Jb04FM777+AIiw+34G5gBbgA3Aq4C/Hfcz8AbmPEE95tvYja3Zt8AN3u3fAVzf0vXLpf9CCGETXa3JRQghxHFIoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE1IoAshhE38f9pHI8rAS3WTAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "time = range(0, 1000, 5)\n",
    "epsilon = tfcompute.outputs[0][:,1]\n",
    "sigma = tfcompute.outputs[0][:,0]\n",
    "plt.figure()\n",
    "plt.plot(time, epsilon, label = 'epsilon')\n",
    "plt.plot(time, sigma, label = 'sigma')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting the potential"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = np.linspace(0.2, 4, 1000)\n",
    "start_model = TrainableLJ(64)\n",
    "start_vals = htf.compute_pairwise(start_model, r)\n",
    "end_vals = hoomd.htf.compute_pairwise(model, r)\n",
    "true_pot = 2 * (r**-12 - r**-6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAewAAAFpCAYAAABeVxsLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAABAq0lEQVR4nO3deXxU9b3/8dd39slkX0iAAAEFcWER415bKy5ordbaulRtpbVUrbX7rV1+2uX2XtvaTVtrvdalrfturdJiRa3iAiKgIPuWQAhJyD5JZvv+/piASLPMZA6Bwffz8cgjs5x8zjlOzJvvcr7HWGsRERGR/ZtrXx+AiIiIDE6BLSIikgUU2CIiIllAgS0iIpIFFNgiIiJZQIEtIiKSBTIObGPMGGPMfGPMCmPMcmPMV/vYxhhjbjbGrDXGLDPGzMh0vyIiIh8kHgdqxIBvWmsXG2PygDeNMfOstSt22+ZMYGLv17HAH3q/i4iISAoybmFba+ustYt7H7cD7wKj99jsXODPNuk1oNAYMzLTfYuIiHxQODqGbYypAo4EXt/jrdFAzW7Pa/nPUBcREZF+ONElDoAxJhd4FPiatbZtiDXmAHMAQqHQUZMnT3bq8GQINuxYiSsRZ1zp4RnXqmkOE+6Jc0hFngNHJiJyYHrzzTcbrbVlfb3nSGAbY7wkw/pea+1jfWyyBRiz2/PK3tfex1p7O3A7QHV1tV20aJEThydDNO/xz+Le+AqnfD3zz+G/HlnKS6sbee17Mx04MhGRA5MxZlN/7zkxS9wAfwLetdb+qp/NngI+2ztb/Dig1Vpbl+m+Ze86LWcMp4S7HKnl87iIxBOO1BIR+SByooV9InAZ8LYxZknva98DxgJYa28DngHOAtYCYWC2A/uVvazVWKI2RqkDtXxuN5GYAltEZKgyDmxr7cuAGWQbC3w5033J8Lq+5S1qRxTxqLVgBvyIB6UWtohIZhybdCYHngvzJ9O58S1IxMDtzaiWz22IxBJYazEZhr+IHLii0Si1tbV0d3fv60PZqwKBAJWVlXi9qf9tVWBLv07IHQfhLohHMg9sT3K6RDRu8XkU2CLSt9raWvLy8qiqqjpg/3FvraWpqYna2lrGjx+f8s9pLXHpVwsJtnjcEI9mXGtnYKtbXEQG0t3dTUlJyQEb1gDGGEpKStLuRVBgS79+27SQS0dWOBPY7t7A1sQzERnEgRzWOw3lHBXY0q9zi6fwnR3NyS7xDHl3dYkrsEUku/zmN78hHA6n/XN33303W7dudew4FNjSr+l5VczqDENCLWwR+eAaSmDH43EFtgyfZhtjvdfj6Bh2jwJbRPZjnZ2dfOxjH2PatGkcccQR/OhHP2Lr1q189KMf5aMf/SgAV111FdXV1Rx++OHccMMNu362qqqK73znO8yYMYP777+fRYsWcckllzB9+nS6ujJfhEqzxKVf9zYs5PbRI1ka6xn4QvsU+D1qYYtIen70t+Ws2DqkW1P067BR+dzw8f7vjzB37lxGjRrF3//+dwBaW1u56667mD9/PqWlyWWkfvrTn1JcXEw8HmfmzJksW7aMqVOnAlBSUsLixYsBuOOOO7jpppuorq525NjVwpZ+nV46nZ83NDkyhq1Z4iKSDaZMmcK8efP4zne+w7///W8KCgr+Y5uHHnqIGTNmcOSRR7J8+XJWrFix670LL7xwrx2bWtjSr0l545jUGYZEPONaXrcmnYlIegZqCe8tkyZNYvHixTzzzDP84Ac/YObM99+waMOGDdx0000sXLiQoqIiLr/88vddnhUKhfbasamFLf1qSfSw0uclFst87EWTzkQkG2zdupWcnBwuvfRSvv3tb7N48WLy8vJob28HoK2tjVAoREFBAfX19Tz77LP91tr955ygFrb065nGJfzv6JG82NNCcYa1fBrDFpEs8Pbbb/Ptb38bl8uF1+vlD3/4A6+++iqzZs1i1KhRzJ8/nyOPPJLJkyczZswYTjzxxH5rXX755Vx55ZUEg0FeffVVgsFgRsemwJZ+nVR2JOWv3ELIZP5rolniIpINzjjjDM4444z3vVZdXc1XvvKVXc/vvvvuPn9248aN73t+/vnnc/755zt2bAps6deYvNGMCXeBtRnX8mkMW0QkIxrDln61xnp42+ejO9qZcS11iYuIZEaBLf16tfldPjO6gtpwQ8a1dFmXiEhmFNjSr6PKpvL7bdsZ6c5sogRolriISKY0hi39KguVU9bVDRmvc6YucRGRTKmFLf1qT0R40++nzYEx7J0Lp6hLXERkaBTY0q+VbRu4fFQ5K8OZ321GXeIiki1uvvlmDj30UIqKirjxxhuHXCc3N9fBo1KXuAxgUvGh3F5XzyFV+RnXcrkMXrdRC1tE9nu33norzz33HJWVlfv6UN5HLWzpV0GwlOO7eyhw6NfE53aphS0i+7Urr7yS9evXc+aZZ/LrX/+aa665BkiuWnbttddywgknMGHCBB555BEAOjo6mDlzJjNmzGDKlCk8+eSTe+3Y1MKWfoUTEZYFgxzc00qpA/V8HgW2iKTh2etg29vO1qyYAmf238192223MXfuXObPn8/TTz/9vvfq6up4+eWXWblyJeeccw6f+tSnCAQCPP744+Tn59PY2Mhxxx3HOeecgzGZT9bdk1rY0q+6zjq+WFHGoq7Mx7AB/B63AltEstYnPvEJXC4Xhx12GPX19QBYa/ne977H1KlTOfXUU9myZcuu95ymFrb0a1TuKO5q7GBC6X/eD3Yo/F4XPbHMb9UpIh8QA7SE9wW/37/rse1dsvnee++loaGBN998E6/XS1VV1ftut+kkBbb0K+gJUp1wg0ONYr/HpZt/iMgBpbW1lREjRuD1epk/fz6bNm3aa/tSYEu/IvEIrwf8TIi2MdqBej4FtogcYC655BI+/vGPM2XKFKqrq5k8efJe25cCW/rVGe3k6nwX1/Vs5xIH6vk9bnWJi8h+b+dtMi+//HIuv/xy4D9vqdnR0QFAaWkpr776ap91dm7jFAW29CvPl8dfw35G+525+N/vcdETVQtbRGQoFNjSL4/LwzQThETm98OGZGC3d8ccqSUi8kGjy7pkQM97LWuj7Y7U0mVdIiJDp8CWAX3D3cozttWRWrqsS0Rk6NQlLgN6wFRSGnXmmkJd1iUiMnQKbBnQZG8+dDkz0zE5S1yBLSIyFOoSlwG9SBdv2y5HaiVniatLXET2Xy0tLdx66637+jD6pMCWAf00toUHPBFHaiXHsNXCFpH9V3+BHYvt+ytcFNgyoNtDU/hahzOtYp/bTSxhieme2CKyn7ruuutYt24d06dP5+ijj+akk07inHPO4bDDDmPjxo0cccQRu7a96aab+OEPfwjAunXrmDVrFkcddRQnnXQSK1eudPzYFNgyoCpfIWWxHkdq+b3JX7eIAltEUjR77myeWPsEANFElNlzZ/O3dX8DoCvWxey5s5m7YS4A7ZF2Zs+dzXObngOgubuZ2XNn80LNCwA0djUOur8bb7yRgw46iCVLlvCLX/yCxYsX89vf/pbVq1cP+HNz5szhlltu4c033+Smm27i6quvHtoJD0CTzmRAL8fbcLkTnOBALb8nGdg90QQ5PgcKiojsZccccwzjx48fcJuOjg4WLFjApz/96V2v9fQ409DZnQJbBnRbeC3BkNehwHYDaBxbRFJ216y7dj32urzvex70BN/3PM+X977nRYGi9z0vDZamvf9QKLTrscfjIZF47+/XzttoJhIJCgsLWbJkSdr106EucRnQTSUn8r8Ng3cjpWJnC1urnYnI/iovL4/29r5XdywvL2f79u00NTXR09PD008/DUB+fj7jx4/n4YcfBpL3yl66dKnjx6YWtgyoIlAEsRjEY+DO7Ndl5xi2VjsTkf1VSUkJJ554IkcccQTBYJDy8vJd73m9Xq6//nqOOeYYRo8e/b5bad57771cddVV/Pd//zfRaJSLLrqIadOmOXpsCmwZ0Ks9DbTnBDk93pN5YKtLXESywH333dfve9deey3XXnvtf7w+fvx45s6duzcPS13iMrCH21dza1EBxDO/FnvXpDO1sEVE0qYWtgzo/406lcTylyDmYGDrntgiImlTYMuAivyFkEhAPPNLFHy7WtgKbBGRdDnSJW6MudMYs90Y804/759sjGk1xizp/breif3K3vdmdz2P5IUcamHvHMNWl7iI9M9au68PYa8byjk6NYZ9NzBrkG3+ba2d3vv1Y4f2K3vZvNZV/KqoyJkxbK9a2CIysEAgQFNT0wEd2tZampqaCAQCaf2cI13i1tqXjDFVTtSS/cs1Y87kyoWPONIlrjFsERlMZWUltbW1NDQ07OtD2asCgQCVlZVp/cxwjmEfb4xZCmwFvmWtXb7nBsaYOcAcgLFjxw7joUl/cgP5yTFsdYmLyDDwer2DLgX6QTVcl3UtBsZZa6cBtwBP9LWRtfZ2a221tba6rKxsmA5NBvJOuI578vOIxzK/J7a6xEVEhm5YAtta22at7eh9/AzgNcakv6irDLvXW9dxU0kRkWg441p+zRIXERmyYQlsY0yFMcb0Pj6md79Nw7FvycylVWfy2sYaAonMJ4D43ApsEZGhcmQM2xhzP3AyUGqMqQVuALwA1trbgE8BVxljYkAXcJE9kKcAHkD8vjz81kIi8zFsYwx+j0tj2CIiQ+DULPGLB3n/d8DvnNiXDK814W28UJDPhd1t5DtQz+9xaZa4iMgQaC1xGdDKjs3cXFxIc6TFkXo+j1td4iIiQ6DAlgGdOfZ03tywmbHuHEfqqUtcRGRotJa4DMjjCyYfxKOO1PN7XWphi4gMgVrYMqCarkZuLSxga7czk/r9HrfGsEVEhkCBLQOq62rkD0UFbHFoDFtd4iIiQ6MucRlQ9ahjWbqxFldloSP1/B4XEXWJi4ikTS1sGZDLuHC5/Y7c/APA73XTrcAWEUmbAlsG1NTVxG+L8ljV3ehIveR12OoSFxFJlwJbBtQeaefukJ8NkVZH6gW8broV2CIiadMYtgyoqqCKt3ZYKCxypF7Q66Jbs8RFRNKmFrYMzuODmDNj2AGvm27NEhcRSZsCWwbUHevmpiAsjDpzHXbA66YrosAWEUmXAlsGlLAJHvLGWB3rcKRewJtcS1w3axMRSY/GsGVAOd4c3oiPBFfQkXoB73v3xA543Y7UFBH5IFALWwbn8Ts3hu1JhrRmiouIpEeBLYO6ybTxXNy5y7oAuhTYIiJpUWDLoP5p23kXZ1rYQV/yV06XdomIpEdj2DKof+bOgNpFjtRSl7iIyNCohS2D8wQcvQ4b1CUuIpIuBbYM6paeGh72xhyptTOw1cIWEUmPusRlUItiLRzsdua66V2XdWkMW0QkLQpsGdQ95afCml+BtWBMRrXUwhYRGRp1icvgPAGwCUhk3i0e1Bi2iMiQKLBlUHe2r+S2wnyIdWdc670WtrrERUTSoS5xGdSaaCvdvt47dvnzMqq1cwxbXeIiIulRYMug/nfcObB0HkS7Mq6ly7pERIZGXeIyOE8g+d2Ba7H9np2zxBXYIiLpUGDLoB7esYyflBQ5MoZtjCHgddEd0xi2iEg6FNgyqPpYmLU+r6OrnXVF1MIWEUmHAlsGdc2Ec7mnbrsjLWxIXtqlSWciIulRYMvgdo1hZz7pDJItbHWJi4ikR4Etg5rXtIyvjCgl7sAscUhOPFMLW0QkPQpsGVR7IsI2j4dYtNORegF1iYuIpE2BLYP65ISP8fDWbfgTznRjawxbRCR9CmwZ3K4xbGcmnQW8Li1NKiKSJgW2DGrhjhVcWV5GfVeTI/UCXrdWOhMRSZMCWwYVNYY2l4tYNOxIPXWJi4ikT2uJy6BOqPwIJ9TVwySfI/X8Xre6xEVE0qQWtgzO5QK3z9ExbK0lLiKSHgW2DGpN8xquKC9leVe9I/U0hi0ikj4FtgzKYIgYF/G4M2uJB71uYglLNK5ucRGRVGkMWwZ1cNHB/LnDBaV+R+rl+N67J7bXrX8zioikQn8tJTUev4Nj2L2BrTt2iYikTIEtg9rRvYPPhWLMjzQ4Um9nCzuswBYRSZkCWwblNm48xoUrHnWk3nuBHXOknojIB4EjgW2MudMYs90Y804/7xtjzM3GmLXGmGXGmBlO7FeGR4G/gD+ZUXwk7syUh6AvWUdd4iIiqXOqhX03MGuA988EJvZ+zQH+4NB+Zbh4/ODQ7TXVJS4ikj5HAtta+xKwY4BNzgX+bJNeAwqNMSOd2LcMj8sStdxHuyO1gl4FtohIuoZrDHs0ULPb89re197HGDPHGLPIGLOoocGZCU7ijEKXj2DcmTHn9y7r0hi2iEiq9qtJZ9ba26211dba6rKysn19OLKbW3Kncl6XU5POkmPYamGLiKRuuAJ7CzBmt+eVva9JtvD4IebMGHbQp+uwRUTSNVyB/RTw2d7Z4scBrdbaumHatzjgKx1v86scZ35dNOlMRCR9jlynY4y5HzgZKDXG1AI3AF4Aa+1twDPAWcBaIAzMdmK/MnwqvLmURHrAWjAmo1petwuv2yiwRUTS4EhgW2svHuR9C3zZiX3JvvH9ER+Ct5+HWA94AxnXC3rddGnhFBGRlO1Xk85kP+bNSX6Phh0pl+Pz6BabIiJpUGBLSm7Y/hLfKitxdPEUdYmLiKROt9eUlIwJlCWvm3YosIM+t2aJi4ikQS1sSckVY07jKy2tDnaJq4UtIpIOBbakxhtMfneshe0hrDFsEZGUKbAlJTfX/JOLRpU718LWLHERkbRoDFtSMjavks7uiKNj2OoSFxFJnVrYkpJPjDuD7+5odqyFrUlnIiLpUWBLahwew87xqoUtIpIOBbak5KGa5/jI2NGEe1odqZfjc9MVjZNIWEfqiYgc6BTYkpIxBeM5tTMMsW5H6gV7b7HZHVMrW0QkFZp0Jik5vvIjHN/UDHFnZnbvfseunffHFhGR/qmFLalxucATcHTSGeie2CIiqVJgS0oWbFnACaPLWBHe5kg93RNbRCQ9CmxJSXmonI/3WPId7xLX4ikiIqnQ4KGk5KDCg/huNAAJZ+rtHLdWC1tEJDVqYUvqvEHHrsPO9ScDu6NHLWwRkVQosCUlNW01HBNo45nIdkfqhXoDu1OBLSKSEgW2pCTfn88FJp9xDl03HfInx7A71SUuIpISjWFLSgr8BXzLNwa6ahypl6sWtohIWtTCltR5g1inrsP2ujFGgS0ikioFtqTEWstRnYv4nafHkXrGGEI+jyadiYikSIEtKTHGMDtnAkd1O7OWOCTHsdXCFhFJjQJbUnZN4VRO6OxwrF7I76GzR5PORERSocCWlFlPkGisGxLOrJ6S61eXuIhIqhTYkrJP1D3DdWUlEHNm8ZSQz6MucRGRFCmwJWWXlBzJGZ1hx1Y7C6mFLSKSMgW2pOyCEcdwergLIp2O1Mv1u+nUzT9ERFKiwJaURT1+wsY42sLWpDMRkdQosCVl3974OJeNLHewha0ucRGRVCmwJWUfH3kil7a1Q8SZS7tCfg+RWIJo3KF7doqIHMC0lrikbOboD0FHp2Mt7J137Ar3xCnI0b8dRUQGor+SkrIet49ml8u5FrYveceuDk08ExEZlAJbUvaHdY8xc+xoR7vEQTcAERFJhQJbUnby2I/yraZmbI8zgb3zFpuaeCYiMjiNYUvKpo88luntHeDQLTbVwhYRSZ1a2JKyHhtjuz+XWE+bI/VC/uQYtgJbRGRwCmxJ2bMbnmXmqGLqupscqfdel7gWTxERGYwCW1I2vWw614cNBdGII/XUJS4ikjqNYUvKqgqqqCIfoj2O1NOkMxGR1KmFLSmLJqJs8fkJR5wZww543fjcLtq6o47UExE5kCmwJWXrW9Yzy7OdV6MtjtXMC3ho71YLW0RkMApsSdnI3JH82DOGQyPOjGGDAltEJFUKbElZvi+f83LGMarHmdtrAuQHvbR1qUtcRGQwCmxJmbWWTW5oijlz8w/Y2cJWYIuIDMaRwDbGzDLGrDLGrDXGXNfH+5cbYxqMMUt6v65wYr8yvCyWjze/zAN+61jNPL9XXeIiIinI+LIuY4wb+D1wGlALLDTGPGWtXbHHpg9aa6/JdH+y77iMixvLTmJi7UMQi4DHl3HN/KBHs8RFRFLgRAv7GGCttXa9tTYCPACc60Bd2Q+dVTyFidEoRJ3pFs8LqIUtIpIKJwJ7NFCz2/Pa3tf2dL4xZpkx5hFjzBgH9iv7wCYbocbjgYhTge0hHIkTiyccqScicqAarklnfwOqrLVTgXnAPX1tZIyZY4xZZIxZ1NDQMEyHJun4Vs3T/Ly40MHA9gJa7UxEZDBOBPYWYPcWc2Xva7tYa5ustTvXs7wDOKqvQtba26211dba6rKyMgcOTZz2rfHnckVrG0ScuSd2fiA5jaKtS4EtIjIQJwJ7ITDRGDPeGOMDLgKe2n0DY8zI3Z6eA7zrwH5lHzi2bDrTeiKOt7A18UxEZGAZzxK31saMMdcA/wDcwJ3W2uXGmB8Di6y1TwHXGmPOAWLADuDyTPcr+8aWeJhWn5fDepxtYWvimYjIwBy5W5e19hngmT1eu363x98FvuvEvmTfum3j33mtvIx5DnWJq4UtIpIarXQmabl00oX8T0MT9Dhzx678oFrYIiKp0P2wJS2HlE+H7h7oaXek3s4WtpYnFREZmFrYkpbt0XbeCOYQ62pxpF6exrBFRFKiwJa0zNv8HF+oKKWje4cj9bxuF0GvW3fsEhEZhAJb0jJz7Ezu7DCEesKO1dQ9sUVEBqcxbElLRaiCCneBYwunQDKwNUtcRGRgamFLWlp7WlkQ8NHa3eJYzcIcH63qEhcRGZACW9KypnkNXzLbeTfa4ljNwqCXlrACW0RkIApsScvk4sn8JecIDu9ybgy7IMerFraIyCAU2JKWXF8u00OV5HU7s3AKQFGOj5ZwxLF6IiIHIgW2pCWaiPJSooON8S6w1pGahUEvnZE4kZjuiS0i0h8FtqQlYRN8uellnsvxO3bHrsKc5Gpn6hYXEemfAlvS4nf7uW/8xZzX3uHYeuIFOT4AWrvULS4i0h8FtqRtStEkShIJcGgcuzCYbGFrpriISP8U2JK2V7vredPvd6yFvbNLXIEtItI/rXQmaftNzbOUFOZzlFOBHUx2ibdoDFtEpF9qYUvafnHkN7mhcYdzXeKhnS1sjWGLiPRHLWxJ29jiSRCPO9Ylnuf34HYZdYmLiAxALWxJ2zudtTwbynGshW2MoSDopUWzxEVE+qXAlrQ9XfM8PykpdqyFDVpPXERkMApsSdsXp87h0cYOx1rYoPXERUQGozFsSVtJsAR8BeDkLTaDXho71CUuItIftbAlbTVtNTyUG6Q93OhYzaIcH82aJS4i0i8FtqRtxY4V/MQbZlt3k2M1i0I+mjsV2CIi/VFgS9pOGn0SzwenMT7s3Bh2cchHZyROdzTuWE0RkQOJAlvSluPNoSxUjqerxbGapbnJ1c6a1MoWEemTAlvS1hHp4N5YA2vjHZBw5h7WxSE/ADs08UxEpE8KbElbOBbmxpa3WOz3QU+rIzVLelvYjZ09jtQTETnQKLAlbaXBUl48/KvJe2J3NTtSsySUDGy1sEVE+qbrsCVtLuOiOK8y+cSpwM5Ndok3qYUtItIntbBlSB5uXsaLwYBjgR3yufF5XJp0JiLSDwW2DMndNfP4e24Iws4EtjGG0pCPJnWJi4j0SYEtQ/LQqXdwY0OTYy1sgOJcHzvUwhYR6ZPGsGVIQvmjkw8cDOySkJ+mDo1hi4j0RS1sGZKX6hZwZ0mZw4Ht0w1ARET6ocCWIVmwdQF35Qaha4djNUvUJS4i0i8FtgzJN6u/yUuRQmfHsEN+uqJxwpGYYzVFRA4UCmwZEq/LiwmWQNjZFjZAY7ta2SIie1Jgy5Csa1nHrzydNHQ5d0/s8vwAANvbux2rKSJyoFBgy5DUddZxb3Q79T3OdYmX5ydXO9verpniIiJ7UmDLkJw46kTenPBZjuhshUjYkZoj8pIt7Po2tbBFRPakwJYhMcZAaETySWeDIzWLcrx43Yb6NrWwRUT2pMCWIYnGo/y84VVeCQag05lxbGMMI/ICbFcLW0TkPyiwZUjcLjePNyxitc8Lndsdqzsi368xbBGRPiiwZUhcxsWrZz/C7NZ2x7rEAcrzAhrDFhHpgwJbhi5Ulvze4WwLW4EtIvKfHAlsY8wsY8wqY8xaY8x1fbzvN8Y82Pv+68aYKif2K/vWn1c/zB3FpY6NYUPyWuy27hjd0bhjNUVEDgQZB7Yxxg38HjgTOAy42Bhz2B6bfQFottYeDPwa+Fmm+5V97+3Gt1maE3K0S3xEXu+12JopLiLyPk7cXvMYYK21dj2AMeYB4FxgxW7bnAv8sPfxI8DvjDHGWmsd2L/sI7/4yC9g7VKHJ531Xovd3s3YkhzH6soHnLVgE5CIg41j41FIxDFYbDxGONKOxxj8xo1NxGjoaiTk8hFy+4jHo2zu2EKxN0SBJ4doPMLq9hrK/YWUenOJJCK807aRsYFSSr25dMd6WNqxkQmBMsq8eYRjXSxu38AhwXLKvLm0R7tY2L6eKaFRlHlyaY2Gea19HUfmjGaEN5emaDsL2jZwXGgMZZ4cGqLtvNy+gZNyx1HqzqEu2sZL7RuYmVtFiTvA5p4WXuzcxOmhCRQZPxsirbwU3sxZOQdR4PKzPtLCgu5azg4eRJ7Lx9rIDl6L1PHxwARyjY9V0R280VPHeYGJBI2Ld6M7WBSt51P+ifgwvBNtYnF8Oxd6J+LBxbJYA0vjjVzkmYgbeCveyNvxJi71TsRay+J4AysSLXzGfRBYyxuJBtbaNi4y48FaXrMNbKSdCxkPWF5hO1tsmAsYiwVesvXU082nGQPAfOppoodPJSoByzxTTxtRzo+PwgJzXdvoIsF5iZFYC8+4txHDck68HIC/ueoBODteDliecG/Da118LF6GBR51byOEm1mxUgAe8tZRmPByWqwUsDzgraPM+jglWgIkuM9Xx6iEn5NjxVgLhSdfw9gjThiWX2MnAns0ULPb81rg2P62sdbGjDGtQAnwvr5UY8wcYA7A2LFjHTg02etCZbBjvWPldq52tq1V49hZJZGAWBdEuyHWRaS7nVi0gxwLxLqoad9CPBamyp0LsW4Wtq3DFYtwlLcI4hGe7FhL0MLp3lKIR7grvJFCazjPVQCJGDdFtzDaGi6O50A8wnWuZg5JGGZ3A/EIXwglqI7EuaqjB2ycc0tzmRnu5trWdkjEObmygvM6OvhqcysAM6rGMLu1jWubW7HAcePHcnVzC1e1tNFjDDOrxvDVHS1c0dpG2BjOqRrDt5qa+VxbO+0uFxeNq+S7jTv4THsHzW43nxs7mhsam/hUeycNHjdXjBnNfzc0cW5HJ9u8Hq6qHMXPtzdyZmeYrT4vXx09kl/XN3BquIsan49vja7gd9u2M6Krm81+H98bVcEf67ZT1t3N+oCf60eWc2ddPaXdPawLBvjvihFMXvIQpT0R1uUE+UV5GdXLnqAiEmVTTpBbyss4ZfnTVEaj1IRyuHVEKWet+DuFsRg1uSH+WFbCJ999hpJYnNq8XP5UWsylq+ZSGk8wLz+Xu0uKuWLVPylIJJhbkMdfi4v48sp/EbSWpwvzuK+4kK+uegEv8ERRHo8U5vHtNf8G4PGifJ4uyOHba1/BYthcnMs/83L45qbXsMCmkjxeDPn52uZFAGwqyeWNoJdra9/EAhtLQ7zj93LNlrewGDaUhljvc3Pl1mUAbCgLUedx8cW6dwDDxrIgLW4XxduSbcQNI4L0GCiuTz5fXx7EAMX17+56nmMtxdvfxWJYXxGgOG4pbliJBdaNDDAqailuXAnA2pFBEpEEJU2rks9HBfB1JyjesQaAbdvPZ7jSymTayDXGfAqYZa29ovf5ZcCx1tprdtvmnd5tanufr+vdpt/Bz+rqarto0aKMjk32rjfq3uChF77LDbUbyfv2OkdqtnVHmfrDf/LdMyfzpY8c5EhNGUAiAT2tEN6BDTdjulugq5nN7ZtoDDcyw50LPR280LmRLZFWLknkQE8H9ySaqE308P2WDoh0cn1RiHq3mz/WJ4dHvlReRofLxb11ydbN5ytGkADu3pbsjblsZDl+a7lj23YwLi4aVU5xAm5t6QG3j0uKfFQmDD+LBMHt5Yv+MAdbH1+3xcTw8F/u7Rxsc7jYlhOxHn7pqmF8IpeTYxVEE/BnzybGxPI4oruMngQ8G9xMRSSfseESIgnDa3mbKAoXUhIuJpIwrC6qoaCriILuYqLWRU1BDXndJQR7iogCDXlbCfYU448UEiNBW2gbvp5ifLFcYlg6g/X4ooUEE/kYj6Ur0EhOvJiAycO4E3T5mwgligm488AVp8u9gxxK8LuCJFwJelwt5LiK8bmCWBOny7SR6yrC6w5gTYKI6SDkLsTj8oMrQY/pJuTOw+v2YU2ChCtCwJWDz+MlQRxcMXyuAF6XhwRxLHG8Lj8utxtLHEjgcXnxuDxAAmssHpcXlys5SmqMxePy4HIZDGAMuIzB5TLJ78bgMsm1E3b/7jLmvW13vbfz/fdqGZLb0fs8+f3977PrMbDnz/SzPea9n9n5PgPts/e1ndu/V/u958PNGPOmtba6r/ecaGFvgd6+i6TK3tf62qbWGOMBCoAmB/Yt+1BrpJVViS7ae1rIi8fAnfmvU37AS17Aw5aWLgeO8APIWuhuhY7t2PY6TGcDdNSzoWUtq9o3MysCdDXzZKyJF1w9/LpuG2D5ZVEhT+WFeHFz8n/du0qKmJ+Twws1W8Gfx7+K83nN5+KSaCn482hy57PNxKDydPCFOKSrhgriMONE8AQ4v2MDEZcLTjkePAG+3FEDHh+UTSfm8vGNlm20RXz821TQFI4zs6Odjm7Lj7strV1Rgl0R1ndFOT0SpSUcpbUryrxYgj/scbq/2eP5rwG/x0Wu30OO382/fB5ygm6CPje1OW4ay9z4vS5GedwEvC58Xjf5HjeVXhcBb/I1f+97fq+bgCe5vc/twut24XUbvG4XPk/yucdtdr3ndu2bP/DyweFEYC8EJhpjxpMM5ouAz+yxzVPA54BXgU8Bz2v8OvudNu40TmuogTVfT45j549ypO7owiBbFdh9i/VAay20bCbcvA5/ax3u1lqWt65jXk89VzXU44/3cG9+LjcVF/HqploC1jKvqIhbCvP4aE8x/mAJnZ4ctiWaSZx0Ka6cYqqjTQSjzXDaRRAo5LJoK58gASOPAWO4Ph7F7XKDSbbAvrHHYV0CROMJalq7qWvtpsvdxdaWbn74bjdNnREa2yfQ1NlDY8cGmsMR3vu///3/tg/53BTm+MgPeikIehhfGqIw6KMgx0tB0Et+wENuwEPI5yHk7/3yuXc9zvG58bp1taocmDIO7N4x6WuAfwBu4E5r7XJjzI+BRdbap4A/AX8xxqwFdpAMdTkQ5I9Ofm/b6lhgVxYFqW3+AAe2tdBeB42roXENNK5ha+MK/tm5kY83bKUkEedvoRy+N6KUp2u3MS6nnHWFxdwTNJxffRljCicwmW4+17WV+BmXQeE4PmljnBptx5M3DlxuPsP7/1X9kd6vnSbscUhet5dYPMGW1jAbmzrZ2BRmU2Mntc1d1LV2sbW1m8aOHvb8Z3jI56Y0z09prp/xpSGqq4opzfVTluujJDf5enHIR1GOl/ygV2ErMgAnWthYa58Bntnjtet3e9wNfNqJfcn+IxwN84P1D3FWTpBT27YAfQ67pG10YZDXN+xwpNZ+LxaBhpWwbRltWxfjrl9OqH4FqxJd/KCshO817eBI62NL2Th+mWOYfPRnKRl5DFN9fr7asYHQ+ZdCbgVnJqJ8DFeyFQwc1fu1UylQSllKh9TWHWVNfTur6ztYta2djU2dbG4KU9McJhp/L5H9HheVRUFGFQY5pCKPkQVBRhUGqCgIMqogwMjCILl+R/7EiAgOBbZ8MPndfjZ0NdDqciVb2A4ZVRikvTtGW3eU/IDXsbr7hbY6qHmNxKbXWLPlFUINa6iMdLPZ4+FjY0bxYwo4b+oFFBeOoaRpAea0L0LVqUxLRHk51kWBvwCAccAVu5X1utL/75RIWDY2dfL2llaWb21j1bZ21tS3s3W3Gfo5PjfjS0McOjKfWUdUMK4kh3ElIapKQozI8+PSuK3IsFFgy5C5XW4e/8ST8NMKaNtznuHQjS4KArC1pYv8iiwP7PZ6WD8fu/Z5Vm55FdtRx2GRKFFPkIsqy/jsIR/i64d8hsryKXxj27+ZUvkRKDqYMuA2vrarjM/tw+f2ZXQo9W3dLN7UzNLaVpbVtvD2llbau2PJ+h4XB5flcuyEEiaV5zGpPJdJ5XmMLgwqlEX2EwpsyYwxybFrB1vYowuTgb2luYvJFfmO1R0W8RhsfhXW/JOWdf+ipnkNUyIRTE4p36go5OBRx3HL8T/CXzGV39a9yqSiSRCqwAXMLpvk6KHUNod5ff0OXt/QxOsbdrCpKQyA122YXJHPOdNGMbWygKmVhUwckYtH48ci+zUFtmTkD0v/wPb8ADc4Gdi9LeysubQrHoUNL8KKp2hb9TT5nU3g8nLd2PFsKZjE3z56K5RP4camtxkZGgk5IwD4cOWHHT2Mzp4Yr6xtZP6qBl5a3bDrv19B0MvRVcVcdtw4qquKmVyRR8DrdnTfIrL3KbAlIz2xHsLeADQ41yVeGvLj87jYsj/PFLcWtiyGJX+Fdx6D7hbuLi7l1hF5vHDkjeQcchZXt63HbdxQejgA08qmOX4Ym5o6mbeinhdWNfDGhh1E4gly/R5OPLiEL540nmMnlHBIeZ66tUUOAApsycjXjvoaNLfAut8lV81yZd6t6nIZxhbnsLGpM+NajutsgiX3wpJ7qWlew23FJVx50IcZM/UzHF1Ywee3vU588pngy2Vq2dS9cgg1O8L8/e06/r6sjre3JJfanFSey+wTqzj5kBEcNa4In0fd2yIHGgW2ZC5/NCSiEG6E3BGOlKwqCbGhcT8K7O0r4bVbiSx7kI5EhOJR1Xhn/pAXNz7IzGNnM2bsKRwOHF5+5F7ZfUs4whNvbeHxt7awtDYZ0tPGFPL9sw5l1hEVjCnWjVJEDnQKbMnIW9vf4qbNj/NTj4eq1lrHAnt8aQ7/XtNAImH3bXfuhn/Dy7+Cdc+T8AQ4f9w4Diufwc9Ov40KYP6xVw3pkqpUJBKWBeuaeHBRDf9Yvo1ILMHho/K57szJfGzKSIW0yAeMAlsy4nf7yQkUEDEGWjbB6BmO1K0qDdETS1DX1r1r1viw2vgyzP9fOja/wovFFXzslB/gOurzzN76YnLiWK+9EdYdPTEeWljD3Qs2snlHmIKgl4uPHsMFR4/h8FEFju9PRLKDAlsycljJYfzfabfDokrYscGxuuNLQgBsbOwc3sDeugTm/T/Y8BLklvPw0Rfwq4ZXOXzap6kKlfDJiZ/ca7uubQ5zz4KNPPBGDe09MY4aV8Q3T5/EGYdXaFa3iCiwxQH+PMgphWbnAruqNBnYGxo7OfHgUsfq9qu9Hp7/Mfate5lfVEbph7/G1JOu4wISHNO6kaqCqr226w2Nndzy/BqeXJK8NO6sKSP5wofGM31M4V7bp4hkHwW2ZOwrz3+F8aWlfMPBFnZFfgC/x8XGvT3xLB6D126FF38GsR4ix13J/7QtZIani597g4SAw3svy3LaxsZObnl+LU8s2YLXbbj8hCo+/6Hx+2YIQET2ewpsyVhFTgUlwVLYttGxmi6X2fszxbe9DU9eQ6xuCc8cfBxnz/od/tKJ3NG6kcq8yr2224b2Hn793GoeXFiDx2WYfUIVcz4ygRF5gb22TxHJfgpsydj3j/s+dP0PrHoleb9mj9+RugeNCLF8a5sjtd4n1gMv/QJe/jUEi3h+5rf5/voHKeyu48NM3Gvd393ROHe+soFb56+jOxrnsuPGcfVHD1JQi0hKFNjijKLxgIWWzVA60ZGSh5Tn8+w72whHYuT4HPpVbVwLj8zGbltG/ZRPUnHmLzktWMSfDj6dY0Ye48w++vCP5dv4ydMrqG3u4tRDy/nuWZM5qCx3r+1PRA48Wg5JMvb0+qeZufJW2lzG0Znih1TkYS2sru/IvJi18Na98McPQ2sNvzz+Ei6KbaDV7cYYs9fCeltrN3P+vIgv/eVNcv0e7r3iWO74XLXCWkTSpha2ZKw8p5zjy48ltvod2LHesbqTK/IAWLWtLbMZ09EuePrrsPR+qDoJPnk7H4+1U77tDfJ8ec4c7B7iCcu9r2/i53NXEUskuO7MyXzhQ+Px6o5YIjJECmzJ2NEVR3N0eTW8/hA0rnKs7pjiHAJeFyu3tQ+9SOsWePAS2PoWc6svZlPldL6UP4pDgEOKD3HsWHdXsyPMNx9ayhsbd3DSxFJ++okpjC3RqmQikhkFtjjDGOyIQzHb33WspNtlmFSex+r6IQb25tfgwcsgGoaL7uO1HQvZUPcqn5/6hb2yQpm1lsff2sL1Ty7HADd9ehrnzxiNMbpTlohkToEtjjjvyfM4JuTmu+vfSY4XOxRSh5Tn8fzK7Vhr0wu+ZQ/BE1cTLayk8+L7KKw8mu/FT8Ng9kpYt4QjfP+Jd/j7sjqOqSrmlxdM01rfIuIoDaiJI2aOncnU4sOgpxXatjpW9/BR+TR1Rqhr7U79hxbcAo99EcYex38ddiJz3voF0UQUn9uH1+18WC+taeFjN7/MP97Zxn/NOoT75xynsBYRx6mFLY645shrYOMr8Mr/wfZ3oWC0I3Wnjy0CYElNC6MGWwEskYB//gBe+z0cfh6c90c+ue0NGrsa91oX+H1vbOZHT62gLM/Po1edwDQtJyoie4kCWxwTLT0YN+DavgImnupIzUNH5uFzu1hS08JZU0b2v2E8Bk9cBW8/xI6jZ7N6+oUc5/FzUuVJjhzHnroicb7/xNs8tngLH5lUxm8unE5RyLdX9iUiAuoSF4fM3TiX6sdOp7ZgJNQvd6yu3+PmsFH5LNnc0v9GifiusGbm9dyUH+QbL32TtsheWCUN2N7Wzaf/uIDH39rC10+dxF2XH62wFpG9Ti1sccSkoknMmTqHwPJ/wdbFjtaePqaQBxfWEIsn8Ox5HXMiAU9duyusOembfKenlQsOuZB8X76jxwGwdns7n7tzIc3hCH/6XDWnTC53fB8iIn1RC1scMaFgAl+e/mVGjDkBGldDeIdjtY8cW0hXNP6f12NbC898E5b8lfaTvsEf83OJJWIU+AuYPmK6Y/vfaeHGHZz/h1fpiSV4cM7xCmsRGVYKbHFMLBGjsXxy8smWNx2re+z4EgAWrGt870VrYe51sOhO+NDXeW7sVG5behsrd6x0bL+7e/btOi6543VKcn08fvUJTKks2Cv7ERHpj7rExTFXPncl3ZFO/mpcULsQJp7mSN2KggAHj8jl5bVNzPnwQcmwnnc9vH4bHPdlmHkD5xnDkSOO3Ct32rrz5Q385O8rmDG2iDs+W63xahHZJxTY4piLJ19MNBGFmhqoecPR2h86uJQHFm6mJxbH/9L/woKbsdVf4DelZZzfXsPY/LGOh3UiYfmfZ97ljpc3cMbh5fz2oiMJeN2O7kNEJFXqEhfHzBw7k1lVs2DssckWdiziWO0TDy6lO5qg/m//nbyX9YzPUvfhb/D42seZXzPfsf3s1B2N85UH3uKOlzdw+QlV3HrJUQprEdmn1MIWx1hrqeuswzf2WEoX3gGbF8CEkx2pfdyEYq7yPs3YpffB1Ivg7N8yyuXisXMfoyRQ4sg+dmoNR/niXxbxxoYdfP+sQ7nipPFaD1xE9jm1sMUxXbEuZj06i4ej28DtgzXzHKudt+QOvuO+j+dcH2L5h67m4bWPAlAaLHU0TGubw5x/2wKWbG7h5ouP5IsfnqCwFpH9ggJbHJPjzeF/TvofzjjobKj6EKz5pzOFF/4J5l5HbcVMrgzP4dYl93PHsjsIR8PO1O+1fGsrn7x1AfVt3dzz+WM4Z9ooR+uLiGRCgS2OOnvC2UwomAATz0hej92wOrOCb/0V/v4NmDSLvEv+gnF7qYhcyj1n3kOO17kbbPx7TQMX3PYqHpfh0atO4PiDnO1mFxHJlAJbHNUd62bB1gU0HnQyGDe89ZehF1v6ADx5DRx0ColP3cXDG+/nw4fk8tjirRT6yhw75kffrGX2XQsZU5zDY1efyKTyPMdqi4g4RYEtjtrasZUvzfsS/255Fw45E5beD/Fo+oWW3A+PXwnjT4IL72VF23p+/9bvOeSgDTSHozy9rC7jY7XWcvO/1vDNh5dy3IQSHr7yeCoKAhnXFRHZGxTY4qjxBeP542l/5PSq0+HIy6CzAd59Kr0iS+5P3sxjwkfg4gfBl8MRpUfw6LmP8q3jL+PgEbnc+fIGEgk75OPsicX55kNL+dW81XxyxmjuvPxo8gLO34JTRMQpCmxxlDGGE0adQMgbSq50VnoIvPCz5B21UvHabbuF9QMsaVnNom2LgOR65S6Xi6tPPogVdW38bdnWIR3jjs4Il93xBo+9tYVvnT6JX356Gj6P/lcQkf2b/kqJ41p7Wvnrir9S07EVPvo9aFwFb9w+8A8l4jD3ezD3OzD5Y3DxA1hPgN8s/g0/ff2nxHcL/E9MH83ho/L5+dxVdPbE0jq2ZbUtnPv7l1lSm7xs65pTJuqyLRHJCgpscVx3rJufLfwZC7YugMPOhYmnw3M/hNpFff9Aay3c83F47fdw7JVwwZ/BG8QYw82n3MzNp9yM2/XeKmMul+GGjx/O1tYu/t8T72Dt4F3j1lruemUD5/9hAYkEPDjnOF22JSJZxaTyx25fqK6utosW9fMHXvZ7dR11jMwdmXzS2Qj/dwp0NcPZv4bDzwOXG9q3wZt3w4JbktuddRNMT65H/viaxzl/4vnvC+o9/Xrean77rzV87dSJfHVm/y3l1fXt/ODxd3hj4w5OPXQEN316GoU5uoGHiOx/jDFvWmur+3pPS5PKXrErrAFCpTD7GXjwUnj0C/D0N8AXgvbeMejJZ8PpP4HiCQD8a/O/+MlrP6Eyt5ITRp/Q7z6unTmR2uYufvPcGtY1dHLdmZMZXRgEki3qZbWt3LNgI08t3UpuwMPPzp/CBdVj1AUuIllJLWzZK6y1/HzhzynLKePzR3w++WIiDiufhvUvQrQLSg+GQ89Nft/D0oalTCubNuh+EgnLLc+v5ffz1xKJJzh4RC4hv4faHWGaOiPk+NxcdPRYrjnlYIp1W0wR2c+phS3DzhjDlo4teF27XSrlcifHtA87t8+feWLtE1SXV1OZV5lSWENyPPurp07kkzNG89TSrSytaSEciXPK5BEcO6GE0w4tpyBHl2uJSPZTC1v2mlgihseV2r8JW3taOfvxszlt3Glcf/z1e/nIRET2T3uthW2MKQYeBKqAjcAF1trmPraLA2/3Pt1srT0nk/1KdtgZ1uta1hH0BBmV+5+zsnf+g7HAX8BfzvxLn9uIiEjml3VdB/zLWjsR+Ffv8750WWun934prD9AumPdfP4fn+dnb/zsP96LxqNc+/y1PLz6YQCqCqrwuTXOLCLSl0zHsM8FTu59fA/wAvCdDGvKASTgCfDLj/yScfnjAFjbvJbtXds5YdQJeN1e4jZOwib28VGKiOz/MhrDNsa0WGsLex8boHnn8z22iwFLgBhwo7X2iX7qzQHmAIwdO/aoTZs2DfnYZP904xs38tiax3jpwpcIeALEE/EBr7UWEfkgGWgMe9DANsY8B1T08db3gXt2D2hjTLO1tqiPGqOttVuMMROA54GZ1tp1A+1Xk84OTE1dTbT0tDA2f+z7Z5CLiEhmk86stacOULjeGDPSWltnjBkJbO+nxpbe7+uNMS8ARwIDBrYcmEqCJZQES/b1YYiIZJ1MJ509BXyu9/HngCf33MAYU2SM8fc+LgVOBFZkuF8REZEPlEwD+0bgNGPMGuDU3ucYY6qNMXf0bnMosMgYsxSYT3IMW4EtIiKShoxmiVtrm4CZfby+CLii9/ECYEom+xEREfmg0+01RUREsoACW0REJAsosEVERLKAAltERCQLKLBFRESygAJbREQkCyiwRUREsoACW0REJAsosEVERLKAAltERCQLKLBFRESygAJbREQkCyiwRUREsoACW0REJAsosEVERLKAAltERCQLKLBFRESygAJbREQkCyiwRUREsoACW0REJAsosEVERLKAAltERCQLKLBFRESygAJbREQkCyiwRUREsoACW0REJAsosEVERLKAAltERCQLKLBFRESygAJbREQkCyiwRUREsoACW0REJAsosEVERLKAAltERCQLKLBFRESygAJbREQkCyiwRUREsoACW0REJAsosEVERLKAAltERCQLKLBFRESygAJbREQkCyiwRUREsoACW0REJAsosEVERLJARoFtjPm0MWa5MSZhjKkeYLtZxphVxpi1xpjrMtmniIjIB1GmLex3gE8CL/W3gTHGDfweOBM4DLjYGHNYhvsVERH5QPFk8sPW2ncBjDEDbXYMsNZau7532weAc4EVmexbRETkg2Q4xrBHAzW7Pa/tfU1ERERSNGgL2xjzHFDRx1vft9Y+6eTBGGPmAHN6n3YYY1Y5WX8/UAo07uuDcNiBeE5wYJ6Xzil7HIjndSCeEzh/XuP6e2PQwLbWnprhzrcAY3Z7Xtn7Wl/7uh24PcP97beMMYustf1OzstGB+I5wYF5Xjqn7HEgnteBeE4wvOc1HF3iC4GJxpjxxhgfcBHw1DDsV0RE5ICR6WVd5xljaoHjgb8bY/7R+/ooY8wzANbaGHAN8A/gXeAha+3yzA5bRETkgyXTWeKPA4/38fpW4Kzdnj8DPJPJvg4QB2J3/4F4TnBgnpfOKXsciOd1IJ4TDON5GWvtcO1LREREhkhLk4qIiGQBBbbDBluG1RhzuTGmwRizpPfrin1xnOkwxtxpjNlujHmnn/eNMebm3nNeZoyZMdzHOBQpnNfJxpjW3T6r64f7GNNljBljjJlvjFnRu2zwV/vYJqs+rxTPKRs/q4Ax5g1jzNLe8/pRH9v4jTEP9n5WrxtjqvbBoaYsxXPKur+BkFy10xjzljHm6T7eG57PyVqrL4e+ADewDpgA+IClwGF7bHM58Lt9faxpnteHgRnAO/28fxbwLGCA44DX9/UxO3ReJwNP7+vjTPOcRgIzeh/nAav7+B3Mqs8rxXPKxs/KALm9j73A68Bxe2xzNXBb7+OLgAf39XE7cE5Z9zew97i/AdzX1+/ZcH1OamE7a9cyrNbaCLBzGdasZq19CdgxwCbnAn+2Sa8BhcaYkcNzdEOXwnllHWttnbV2ce/jdpJXZuy5smBWfV4pnlPW6f3v39H71Nv7teekonOBe3ofPwLMNIOsBb0vpXhOWccYUwl8DLijn02G5XNSYDsr1WVYz+/tinzEGDOmj/ezzYG8/Ozxvd17zxpjDt/XB5OO3m65I0m2cnaXtZ/XAOcEWfhZ9XazLgG2A/Ostf1+VjZ5iWwrUDKsB5mmFM4Jsu9v4G+A/wIS/bw/LJ+TAnv4/Q2ostZOBebx3r/KZP+zGBhnrZ0G3AI8sW8PJ3XGmFzgUeBr1tq2fX08ThjknLLys7LWxq2100muAHmMMeaIfXxIGUvhnLLqb6Ax5mxgu7X2zX19LApsZw26DKu1tsla29P79A7gqGE6tr0p5eVns4m1tm1n955NriXgNcaU7uPDGpQxxksy2O611j7WxyZZ93kNdk7Z+lntZK1tAeYDs/Z4a9dnZYzxAAVA07Ae3BD1d05Z+DfwROAcY8xGksOcpxhj/rrHNsPyOSmwnTXoMqx7jBWeQ3I8Lts9BXy2d/bxcUCrtbZuXx9UpowxFTvHoYwxx5D8/2W//mPZe7x/At611v6qn82y6vNK5Zyy9LMqM8YU9j4OAqcBK/fY7Cngc72PPwU8b3tnNu2PUjmnbPsbaK39rrW20lpbRfJv+vPW2kv32GxYPqeMVjqT97PWxowxO5dhdQN3WmuXG2N+DCyy1j4FXGuMOQeIkZzwdPk+O+AUGWPuJzkLt9Qkl6K9geRkEqy1t5Fcxe4sYC0QBmbvmyNNTwrn9SngKmNMDOgCLtqf/1j2OhG4DHi7dxwR4HvAWMjazyuVc8rGz2okcI8xxk3yHxgPWWuf3uPvxZ+Avxhj1pL8e3HRvjvclKRyTln3N7Av++Jz0kpnIiIiWUBd4iIiIllAgS0iIpIFFNgiIiJZQIEtIiKSBRTYIiIiWUCBLSIikgUU2CIiIllAgS0iIpIF/j9OiLyoitzOUgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8,6))\n",
    "# is per-particle energy, so grab only 0 particle\n",
    "plt.plot(r, start_vals[2][:,0], label='start')\n",
    "plt.plot(r, end_vals[2][:,0], label='final')\n",
    "plt.plot(r, true_pot, label='true', linestyle=':')\n",
    "plt.legend()\n",
    "plt.ylim(-1,2)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}