{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy import spatial\n",
    "import matplotlib.pyplot as plt\n",
    "from ipywidgets import interact\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lennard-Jones\n",
    "\n",
    "Use proper vectors (r is a vector as function arguments and force should have direction)\n",
    "\n",
    "$$ V(r) = 4 \\varepsilon \\left[\\left( \\frac{\\sigma}{r} \\right)^{12} - \\left( \\frac{\\sigma}{r} \\right)^{6} \\right] $$\n",
    "\n",
    "$r^{-12}$ term is repulsion due to Pauli Exclusion Principle and $r^{-6}$ term is attraction due to dipole-dipole interaction (van der Waals) between non-polar atoms.\n",
    "\n",
    "$V(\\sigma)=0$. and $V(2^{1/6} \\sigma) = V_\\text{min} = -\\varepsilon$\n",
    "\n",
    "$$ F(r) = -\\nabla V = \\frac{24 \\varepsilon}{\\sigma} \\left[2 \\left(\\frac{\\sigma}{r} \\right)^{13} - \\left( \\frac{\\sigma}{r} \\right)^{7} \\right] $$\n",
    "\n",
    "For Argon the values are\n",
    "\n",
    "- $m = 6.69 \\times 10^{-26}$ kg\n",
    "- $\\varepsilon = 1.65 \\times 10^{-21}$ J\n",
    "- $\\sigma = 3.4 \\times 10^{-10}$ m\n",
    "\n",
    "We'll choose all units as 1.0. $m = 1, \\sigma=1, \\varepsilon=1$. Which gives a time unit\n",
    "\n",
    "$$\\tau = \\sqrt{\\frac{m \\sigma^2}{\\varepsilon}} = 2.7 \\times 10^{-12} \\,\\text{s}$$\n",
    "\n",
    "This means a time duration $\\Delta t = 1.0$ correspondes to 2.7 picoseconds in a simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlcAAAFwCAYAAACVel6XAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl8FdX5x/HPZN83CBBICAEUQVZFEQUFF0Cx4oqKAlFr\nrVWL/upWtRbU2rrUWlvbqlVAqC2IC4qKiBACIpuA7LIEwhIChGxk3+b3x5CEsN7kLnNv8n2/XvOa\nOffOnfsEj/DkzDPnGKZpIiIiIiKu4Wd3ACIiIiLNiZIrERERERdSciUiIiLiQkquRERERFxIyZWI\niIiICym5EhEREXEhlyVXhmH4GYax2jCMz1x1TRERERFf48qRqwnAJhdeT0RERMTnuCS5MgwjEbgG\n+LcrriciIiLiq1w1cvUX4DFA072LiIhIi+Z0cmUYxkjggGmaawHj6CYiIiLSIhnOri1oGMaLwJ1A\nFRAKRAIfm6Y57rjzNKolIiIiPsM0zSYNGDmdXDW4mGFcBvzGNM3rTvKeqUWim5e2beHgQdi/H9q1\nc911J06cyMSJE113QWBv4V6S/pJE+8j27Pu/fS69ttjLHf1Fmif1FWkMwzCanFxpnitpsuhoa5+f\n79rrfvrpLtdeEIgIigDgSPmRxn1QfxF7vV27dtkdgvgI9RXxlABXXsw0zUXAIldeU7xXTIy1d3Vy\n5Q61yVVRRRGmaWIYKg0UERH30MiVNFltclVQ4Nrrvv56qmsvCAT4BRASEIKJSUllieMf1MiV10tN\nTbU7BPER6iviKUqupMliY619To5rrztkyBDXXvCoyKBIwBq9kubDXf1Fmh/1FfEUJVfSZG3aWPtD\nh1x73dTUNNde8Ki6uquKRtRdaeTK66WlpdkdgvgId/aVTp06YRiGNh/cOnXq5PL+4NKaK2lZapOr\ngwftjcNRkcEauRIR98jMzERPxPsmd9TgauRKmsxdI1dTpgxx7QWPOrao3WEaufJ6utUjjlJfEU9R\nciVN5nMjV0drrho9HYOIiEgjKLmSJnNXcuXumiuNXDUvqrkSR6mviKcouZIm87WRqyYVtIuIiNf5\n4IMPGDFihEPnTpo0ibFjx7o5ooaUXEmTuSu5clfNVZOmYtDIlddTHY04qqX2lZSUFBYsWHDC65mZ\nmfj5+VFTU+Ox72yKk8U5ZswY5s6d6/A1PD1xtJIrabKoKAgKgqIiKGnEvJx2afISOCIizZQvrFZR\nu6qGLz2NqeRKmsww3PPEoLtqrpo0FYNGrrye6mjEUeorjTN06FCeeuopBgwYQHR0NDfccAP5x6x3\n9tlnn9GzZ0/i4uK4/PLL+emnnwAYN24cu3fv5mc/+xlRUVG8+uqrACxbtoxLLrmE2NhY+vXrx6JF\nixp817PPPsugQYOIiopixIgR5ObmAnDZZZcBEBMTQ1RUFMuXL2fq1KkMHjy47vMPP/wwHTt2JDo6\nmgsuuIAlS5a4/c/ndJRciVN8qe6qSQXtIiIt2LRp05gyZQrZ2dn4+/vz0EMPAbB161bGjBnDG2+8\nwaFDh7j66qu59tprqaqq4v3336djx47MmTOHwsJCHn30UbKysrj22mt59tlnycvL49VXX+Wmm27i\n8OHDdd/13//+l6lTp3Lo0CHKy8vrkrL09HQACgsLKSwsZMCAAUDDUbcLL7yQdevWkZeXx5gxY7jl\nlluoqKjw1B/TCZRciVPckVy5u+ZKM7Q3Ly21jkYaz66+Yhiu2zxt7NixdO/endDQUJ5//nk+/PBD\nTNNk5syZXHvttVx++eX4+/vz6KOPUlpaytKlS+s+e+xtvOnTpzNy5EiGDx8OwBVXXEH//v358ssv\n686566676NKlC8HBwYwePZq1a9c2iOV0twXHjBlDTEwMfn5+PPLII5SXl9eNpNlByZU4RSNXIiLN\nV1JSUt1xcnIylZWV5OTkkJWVRXJyct17hmGQlJTEvn37TnqdzMxMZs6cSVxcHHFxccTGxvLdd9+R\nnZ1dd067du3qjsPCwigqcvzv6ldffZUePXoQGxtLbGwshYWF5Lh64dtGUHIlTnFHcuXumiuNXDUv\nqqMRR9nVV0zTdZun7dmzp+44MzOTwMBAWrduTfv27cnMzDzh3MTERODEQvmkpCTGjRtHbm4uubm5\n5OXlceTIER577LEzxnCmovvFixfzyiuvMGvWLPLy8sjLyyMqKsrWAnglV+IUjVyJiHi/iooKysvL\n67bq6mqHPjd9+nS2bNlCSUkJv//977nlllswDIPRo0fzxRdfsHDhQqqqqnj11VcJCQlh4MCBgDUK\nlZGRUXedO++8k88//5x58+ZRU1NDWVkZixYtIisr64wxxMfH4+fnx44dO076flFREYGBgbRq1YqK\nigqee+45jhyx96lwJVfilPh4a+9TNVeNmYpBI1deTzVX4qiW3FdGjhxJWFgYoaGhhIWFMWnSpLr3\nTjcyNHbsWMaPH0/79u2pqKjgr3/9KwBnn30206dP58EHHyQ+Pp4vvviCzz//nICAAACefPJJnn/+\neeLi4njttddITExk9uzZvPjii8THx5OcnMyrr75aN3fV6WIIDQ3l6aef5pJLLiEuLo4VK1Y0eH/4\n8OEMHz6cs88+m5SUFMLCwhrczrSD4alhM8MwTF+ao0Ic8+WXMHIkDBsGX39tdzSntyN3B13/1pWU\nmBQyJmSc+QMiIg7ytXmYANavX89ll11WN+XB8YYOHcrYsWO5++67PRyZZ53qv93R15v0GIFGrsQp\nvjTPVZOWv9HIlddTzZU4Sn2lXk1NDTNmzKB///52h9IsBdgdgPg2X6q5atIkoiIizVBycjIJCQlM\nnjz5lOf4wuzt3kq3BcUppaUQFgaBgVBebs88LI4yTZOA5wOoMWuo/F0lAX763UJEXMMXbwuKRbcF\nxeuEhkJkJFRWQkGB3dGcnmEYemJQRETcTsmVOM3VtwbdVXMFTZiOQTVXXk91NOIo9RXxFCVX4jSf\nqrtqynQMIiIijaDkSpzm6uTKXfNcgUaumqOWPHeRNI76iniKkitxmk+NXDVlCRwREZFGUHIlTlPN\nldhJdTTiKPWVxjl48CCXXnop0dHRDq0BKPWUXInTfGrkKkhzXYlIy5KSksKCBQtOeD0zMxM/P7+6\nJWiO9/bbb9OmTRsKCgp45ZVX3B1ms6LkSpzmizVXheWFjn1AI1deT3U04ij1lROdbqLQzMxMevTo\n0aTrOrowdHOl5Eqc5ksjVzEhMQDkl+XbHImIiPe66667mDp1Ki+99BJRUVEsWLCAiooKHn74YTp0\n6EBiYiKPPPIIlZWVACxatIikpCRefvllEhIS6tYjnD17Nv369SM6OpqzzjqLefPmAVBYWMjPf/5z\n2rdvT1JSEr/73e+a1SSsSq7Eab5UcxUXGgdAXmmeYx/QyJXXUx2NOEp9xXGTJ0/mjjvu4IknnqCw\nsJDLL7+cF154gRUrVrBu3Tp+/PFHVqxYwQsvvFD3mezsbPLz89m9ezdvv/02K1asYPz48fz5z3+m\noKCA9PR0OnXqBMD48eMJCgoiIyODNWvW8M033/Dvf//bpp/W9bT+hzjNl0auYkNiAcgtPfkq8CIi\nrmZMct26YObv7Rvd+eCDD3jzzTdp1aoVAL///e/55S9/yaRJkwDw9/dn0qRJBAYGAvDee+9xzz33\ncPnllwOQkJBAQkICBw8e5KuvvqKgoIDg4GBCQkJ4+OGHefvtt7n33nvt+eFcTMmVOK1VK2tNwdxc\nqKqCACd7lTtrrmJDreQqr0wjV82F6mjEUeorzsnKyqJjx4517eTkZLKysura8fHxdYkVwJ49exg5\ncuQJ18nMzKSyspKEhATAWvfVNM0G1/Z1TidXhmEEA+lA0NHrzTJNc5Kz1xXf4e9vJVg5OXD4MLRt\na3dEp1Y7cuVwciUi4iQ7R5tcqUOHDmRmZtK9e3fASpLat29f9/7xxfFJSUns2LHjhOskJSUREhLC\n4cOHT1tQ78ucrrkyTbMcGGqaZj+gL3C1YRgXOh2Z+BRX3hp0Z81V3ciVaq6aDdXRiKNacl+pqKig\nvLy8bmvK03y33XYbL7zwAjk5OeTk5PD8888zduzYU55/zz33MHnyZBYuXIhpmmRlZfHTTz/Rrl07\nhg0bxiOPPMKRI0cwTZOMjAzS09Od+RG9iksK2k3TLDl6GIw1etU80nRxmK/UXWnkSkRaopEjRxIW\nFkZoaChhYWF1dVJw6ukYjn/9mWeeoX///vTu3Zs+ffrQv39/nn766VN+5wUXXMDkyZN5+OGHiY6O\nZsiQIezevRuA999/n4qKCnr06EFcXBy33HIL2dnZLvhJvYPhikcfDcPwA34AugBvmqb525OcYzan\nxyyloVtvhZkz4YMP4Pbb7Y7m1HJLc2n1ciuig6PJf1LTMYiIaxiG4XNTCaxfv57LLruM3NyW/YDP\nqf7bHX29SfctXTVyVXP0tmAiMMAwjKbNOiY+y1dGrqKDowEoKC+guqZlT3InIi1XTU0NM2bMoH//\n/naH0iy59GlB0zQLDcNYCIwANh3/fmpqat0cFzExMfTt27fu6Y3ae+Fq+2a7uNhqHzzo/PVSU9NI\nTcVt8caExJC/OZ858+YwasSo05+flgYTJ9r+56v2qdu1x94Sj9re2659zZ3X9xXJyckkJCQwefJk\nu0PxCrX/DdPS0ti1a5fT13P6tqBhGK2BStM0CwzDCAW+Bv5kmuaXx52n24LN2FtvwS9/CT//Obzz\njnPXSk1Nc+t0DJ3/2pmd+TvZ9tA2usZ1Pf3JEyeqqN3LpaWl1f0jJ3I67uwrvnhbUCzeelswAVho\nGMZaYDnw9fGJlTR/rrwt6M7EChr5xKASK6+nxEocpb4inuL0bUHTNNcD57kgFvFhvlJzBXpiUERE\n3EtrC4pL+Mo8V6CRq+bGF+tdxB7qK+IpSq7EJTRyJSIiYlFyJS4RFQVBQVBUBCUlZz7/dNxecxWi\nkavmRHU04ij1FfEUJVfiEoZRP3p16JC9sZxJXGgcYE0oKiIi4mpKrsRlXJVceazmypHbghq58nqq\noxFHtcS+cvXVVzPxJH+PzZ49m4SEBGpqalz2XXl5edxwww1ERESQkpLCf//7X5dd29couRKX8ZW6\nK9VciUhLMX78eKZPn37C69OnT2fs2LH4+bkuDfjVr35FSEgIhw4dYvr06dx///1s3rzZZdf3JUqu\nxGVclVxpnitpDNXRiKNaYl+5/vrrOXz4MEuWLKl7LT8/nzlz5jBu3DiXfU9JSQkff/wxL7zwAqGh\noVxyySWMGjWKadOmuew7fImSK3GZ+Hhrr5ErERHvEBISwi233ML7779f99qMGTPo3r07PXv2POln\nHnjgAWJjY4mLi6vb1x737dv3pJ/ZunUrgYGBdOnSpe61Pn36sHHjRtf+QD5CyZW4jKtGrjTPlTRG\nS6yjkaZpqX1l/PjxfPjhh1RUVAAwbdo0xo8ff8rz33zzTfLy8sjNza3b1x6vXbv2pJ8pKioiKiqq\nwWtRUVEcOXLEdT+ID1FyJS6jmisRkZM7fplSZ9uNcckllxAfH8+nn35KRkYGK1euZMyYMU272ClE\nRERQWFjY4LWCggIiIyNd+j2+wumFmx3+Ii3c3Ox9+SWMHAnDh8PcuXZHc2o1Zg0BzwVgYlL5u0oC\n/JxeBUpEWjhvX7j5+eefZ9myZQwYMIBVq1bx2WefnfLc+++/n+nTp2MYDdcsNk2TTp06sX79+hM+\nU1JSQlxcHBs3bqy7NThu3DgSExN58cUXXfvDuJg7Fm5WciUus2oVXHAB9OsHq1fbHc3pxb0UR15Z\nHoceO0TrsNZ2hyMiPs7bk6vMzEzOPvts2rZty1/+8hduuukml3/HmDFjMAyDd955h9WrV/Ozn/2M\npUuX0r17d5d/lyu5I7nSbUFxGV+puYJG1F2p5srrtdQ6Gmm8ltxXkpOTufjiiykpKeG6665zy3e8\n+eablJSU0KZNG+68807+9a9/eX1i5S66HyIuc+zTgqZpzdrurVR3JSItzcKFC916/djYWD755BO3\nfoev0G1BcamoKDhyBPLyICbG7mhO7appVzE/Yz5z75jL8K7D7Q5HRHyct98WlFPTbUHxenpiUERE\nWjolV+JSrkiuPFJzFaKaq+aiJdfRSOOor4inKLkSl/KZkaujBe25pbk2RyIiIs2Nkitxqdrk6tCh\npl/D3WsLAsSFxgEO3BbUyJXXa4nrxUnTqK+Ipyi5EpfymZErR28LioiINJKSK3Epn6m5CnWwoF0j\nV15PdTTiKPUV8RQlV+JSPjdypacFRUTExZRciUu5IrnyRM2VZmhvPlRHI45SXxFPUXIlLnXsLO3e\nTCNXIiKW/fv3M3fuXLvDaFaUXIlL1Y5cHTjQ9Gt4subqjFMxaOTK66mORhylvnKidevWMWfOHJYv\nX878+fPtDqfZUHIlLtW6NQQGwuHDUFpqdzSnFhUchYFBUUURldWVdocjImIL0zTrln7R8j2uo7UF\nxeW6dIGMDNiyBbp1szuaU2v1cityS3M5+OhB4sPj7Q5HRHyYN68t+MMPP7Bt2zb+9a9/ccMNN/D+\n++/z6aefkpSUBEB2djY//vgjw4e3zHVWtbag+ITkZGufmWlvHGeiuisRaQmCgoLo0aMHgYGBTJgw\ngaVLl9YlVgDt2rVrsYmVuyi5EpdzNrnyRM0VOPjEoGquvJ7qaMRRtvaViRMb/n3ibLsRevXqxbx5\n87j55psBCA4ObtJ1xHEBdgcgzY9GrkREjnN8YuRsu5HmzZvHW2+95dQ1xHEauRKXcza58sQ8V6CR\nq+ZCcxeJo1pyXykrKyMlJcXuMFoMJVficr42cnXG6RhERHxcenq63SG0KEquxOV8pubKkduCGrny\neqq5Ekepr4inOJ1cGYaRaBjGAsMwNhqGsd4wjF+7IjDxXUlJYBiwbx9UVdkdzanFhcYBDiyBIyIi\n0giuGLmqAv7PNM1zgYHAA4ZhnOOC64qPCgqChASorrYSrMbyeM2VRq58Wkuuo5HGUV8RT3E6uTJN\nM9s0zbVHj4uAzUAHZ68rvq321uCuXbaGcVp6WlBERNzBpTVXhmF0AvoCy115XfE9ztRdaZ4raQzV\n0Yij1FfEU1yWXBmGEQHMAiYcHcGSFswXnhjUyJWIiLiDSyYRNQwjACuxmmaa5uxTnZeamkqnTp0A\niImJoW/fvnX3wGt/o1C7ebTLy612ZmbjPz9lyhCPxJtdlA1YUzGc8vyjI1d2/3mqfer2kCGe6S9q\nq32mtviu2v+GaWlp7HJBPYtLFm42DON9IMc0zf87zTlauLkF+eoruOYauPJK+OYbu6M5uYKyAmJe\niiE8MJyipzTYKiJN580LN8vpeeXCzYZhXALcAVxuGMYawzBWG4Yxwtnrim/zhZqryOBI/Aw/iiuL\nqayuPPlJqrnyeho1EEepr4inOH1b0DTN7wB/F8QizUhtcrV7N9TUgJ8XTlfrZ/gRGxLL4dLD5JXl\n0Sa8jd0hiYhIM+CF/+RJcxAeDq1aQXk5HDzYuM96ap4rcOCJQY1ceb3auheRM2mpfWXv3r2MHz+e\ndu3aMW3atLrXt2/fTs+ePXn88cfJz8+3McLmR8mVuI2eGBQRsV9iYiITJkwgIiKCsWPH1r0eGhrK\ns88+y8svv0xMTIyNETY/Sq7EbZqaXHmq5go0ctUcqI5GHNWS+0pKSgp79uyhpqam7rWZM2cyevRo\nG6NqvpRcidv40shVbmmuzZGIiLhPbGws4eHh7N69G4CPPvqIm266yeaomi8lV+I2TU2uPFpzdabb\nghq58nottY5GGq+l95WUlBR27tzJ/v37qaqqomPHjqc8Nz09nffee8+D0TUvLplEVORkfGLkypEl\ncEREnLVuorXvPdE17Sbo3LkzO3bsYO3atTzyyCOnPffSSy/l0ksvbfJ3tXQauRK38YWaq7jQOEAj\nV76sJdfRSOPY2ld6T2yYGDnbboKUlBTeeecdrr/+egDeffddHnjgAZYsWcJTTz3F0qVL+f3vf09u\nbi7/+c9/ePnllxuct3jxYp555hmnYmgplFyJ2xybXHnrxMV6WlBEWoqzzjqL8ePHk5KSAsCVV15J\nTEwMgwYNYufOnVx88cWkpKSwadMmBg0aRFlZWYPzBg8eTFZWlp0/gs/QbUFxm7g4a76rwkLIz4fY\nWMc+p3mupDFaeh2NOK6l95V77733hNeCgoIACAkJAawlX2qfKDx2SZja87TEj2M0ciVuYxhwdJ1u\nr6270siViLRUixYtYs2aNfz4449s2LCBNWvWsGTJEhYtWsSiRYv48ccfqaysPOG8DRs22B2613PJ\nws0OfZEWbm6RRo6EL7+ETz+FUaMc+0xqaprHRq/WZq+l31v96NWmF+vuX3fiCRMnavTKy6WlpbX4\nEQlxjDv7ihZu9l1euXCzyOl4+xODmudKRERcTcmVuFVTkitbaq70tKDP0qiVOEp9RTxFyZW4lbeP\nXEUGReJv+FNSWUJFdYXd4YiISDOg5ErcqinJlSfnuTIM4/RPDGrkyutpnitxlPqKeIqSK3Erbx+5\nAj0xKCIirqWnBcWtamogNBQqKqC4GMLC7I7oRAP+PYAV+1aw9O6lDEwaaHc4IuKD9LSg79LTguJz\n/PwgKck6ProYu9fRyJWIiLiSkitxu8beGvRkzRXUPzF40ukYVHPl9VRHI45yZ19JTk7GMAxtPrgl\n1/4j5UJa/kbcztvrrupGrk61BI6IyBns2rXL7hDEi2jkStyuscmVJ+e5gjPcFtTIldfT3EXiKPUV\n8RQlV+J23j5yFRcaB2jkSkREXEPJlbhd7eLNjo6a21VzpZEr36SaK3GU+op4ipIrcTtvH7nS04Ii\nIuJKmudK3K6yEkJCrOOyMggMtDee46XtSmPo1KEM7jiY9LvS7Q5HRES8gOa5Eq8WGAjt21sTiu7d\na3c0J6oduTrpVAwiIiKNpORKPKIxtwZVcyWNoToacZT6iniKkivxCG+uu9I8VyIi4kpKrsQjGpNc\neXqeq4igCAL8AiitKqWsqqzhmxq58nqau0gcpb4inqLkSjwiJcXab99ubxwnYxgG8WHxABwqPmRz\nNCIi4uuUXIlHnHuutd+w4cznerrmCqBdRDsAsouyG76hkSuvpzoacZT6iniKkivxiJ49rf2mTVBd\nbW8sJ3PK5EpERKSRlFyJR0RFQceOUF4OO3ac/lxP11yBRq58mepoxFHqK+IpLkmuDMN41zCMA4Zh\nrHPF9aR5qh29Wr/e3jhORiNXIiLiKq4auZoMDHfRtaSZ6tXL2p+p7ko1V9IYqqMRR6mviKe4JLky\nTXMJoEmC5LR8YeTqQPEBmyMRERFfp5or8RhHR67sqLlqG94W0MiVL1IdjThKfUU8RcmVeEy3buDv\nD9u2QWmp3dE0pJorERFxlQBPfllqaiqdOnUCICYmhr59+9b9JlF7L1zt5t0+66whbNkC06encdZZ\nJz8/NTWN1FQ8Gt+21dtgJ2QHZjd8Py0NJk70mj8/tU9s1x57Szxqe2+79jVviUdt72rXHu/atQtn\nGaZpOn0RAMMwOgGfm6bZ6xTvm676LvFdo0fDhx/C1KkwbtzJz0lNTfP4rUHTNAl7MYyyqjKO/PYI\nEUER1hsTJ+rWoJdLS0ur+0tS5HTUV6QxDMPANE2jKZ/1c1EAHwBLgbMNw9htGMZdrriuND+1Re2n\nq7uyo+bKMIyT3xpUYuX19I+lOEp9RTzFJbcFTdMc44rrSPNXW9TurU8M7srfxYGiA3SN62p3OCIi\n4qNcMnIl4ihHRq7smOcKTvHEoEauvN6x9RIip6O+Ip6i5Eo8qnNnCA2FvXshP9/uaBrSE4MiIuIK\nSq7Eo/z9oUcP6/hUo1d21FzBKZIrjVx5PdXRiKPUV8RTlFyJxzlya9AOGrkSERFXUHIlHnemona7\naq7qkqtijVz5EtXRiKPUV8RTlFyJx3n7yNWBIq0vKCIiTeeySUTP+EWaRFSOysqCDh0gNhYOHwaj\nSVO0ud7OvJ10fqMzSVFJ7H5kt93hiIiIjWyfRFSkMRISrMQqLw/277c7mnptI+qnYtAvAiIi0lRK\nrsTjDKP+1uDJ6q7sqrkKCwwjKjiKyppK8sryrBdVc+X1VEcjjlJfEU9RciW2qC1q99a6Kz0xKCIi\nTaXkSmxxupEru+a5gpMkVxq58nqau0gcpb4inqLkSmzh7SNXemJQRESaSsmV2OLcc639pk1QXd3w\nPbtqrgDahWvkyteojkYcpb4inqLkSmwRG2tNx1BaChkZdkdT79gnBkVERJpCyZXY5lS3Br2i5qpY\nI1e+QnU04ij1FfEUJVdim9MVtdtFTwuKiIizlFyJbU41cmVrzZWeFvQ5qqMRR6mviKcouRLbePPI\nlZ4WFBGRptLagmKb0lKIiLBmbC8qgpAQuyOCyupKgl8IxjAMKp6pwN/P3+6QRETEBlpbUHxSaCh0\n7WpNxbBpk93RWAL9A2kV1ooas4ZDJYfsDkdERHyQkiux1YUXWvslS+pfs7PmCo6ru1LNlddTHY04\nSn1FPEXJldiq9slob/o7T08MioiIM1RzJbbKyIAuXSAuDg4dAj8vSPfHfjKW6eumM3nUZFL7ptod\njoiI2EA1V+KzUlIgKQlyc73nqcHaJXD0xKCIiDSFkiuxlWGceGtQNVfSGKqjEUepr4inKLkS2w0d\nau0XLrQ3jlp16wsWq+ZKREQaTzVXYrudO6FzZ2sx55wc++uu5mfM56ppVzGk0xAWjveSjE9ERDxK\nNVfi0zp1go4dIS8P1q2zOxo9LSgiIs5RciW2M4z6W4Npaaq5ksZRHY04Sn1FPEXJlXiF2qJ2b6i7\niguNI8AvgPyyfKpqquwOR0REfIxqrsQr7NplTcsQE2PVXfnbvKRf4muJ7Duyj8yHM+kY3dHeYERE\nxONUcyU+r1Mna8vPV92ViIj4NiVX4jVqbw3+4hdpNkZhqZ2OIfqPf7E5EjkT1dGIo9RXxFNcklwZ\nhjHCMIwthmFsNQzjCVdcU1qe2uQq2wsGi2pnaS+qKLI5EhER8TVOJ1eGYfgBfweGA+cCtxuGcY6z\n15WWpzYQSSKbAAAgAElEQVS5OnJkCNXVtoZSd1vwizsusDcQOaMhtR1H5AzUV8RTXDFydSGwzTTN\nTNM0K4H/AaNccF1pYZKTraL2ggJYu9beWGqTK60vKCIijeWK5KoDsOeY9t6jr4k0mjXfVRp2l0bU\nJleXTllgbyByRqqjEUepr4inBHjyy1JTU+nUqRMAMTEx9O3bt26YtrbTq632e+/BrFlpnH++ffHs\nX78fdtbXXHnTn4/aaqvdtHYtb4mnMe2aGhg0aAhVVbBwYRrV1XDRRVZ7yRKr3b+/VVLx/fdWu18/\nq71ihfX5Xr2s9urV1vvnnmu1162z2uecY7U3brTaZ51ltbdssT7fubPV3rrVanfqZLUzMqx2UtIQ\nampg5840TBPat7fau3db7yckWOfv22e9Hx9vvb9/v9Vu3dpqHzhgtWNjrXZOjvX5mBirnZtrtaOi\nrHZ+vnV+RITVLiiw/vzCwqz2kSPW+6GhVru42GoHBw/BNKGkxGqbJpSXp1FVtQtnZ45yep4rwzAu\nAiaapjniaPtJwDRN86XjztM8V3JGe/ZYS+FERcHhwxDg0fS/3tbDW+n29250ju3Mjl/vsCcIEXEZ\n04SKCigvb7hVVNS/fuz+dFtl5YnHlZWOb1VVJx5XVdVvx7erqnD6H3tpiqbPc+WKf7pWAl0Nw0gG\n9gO3Abe74LrSAiUlQZcusGOHVXfVv789cRw7z5VpmhhGk/7/EpHj1NRAWRmUlDTcSkvr97Vbbbus\n7NT78nJrf+xWmzgde1xRYfdP7rzAQGuC5cBA6xdPf39rf6pjf/+Gx8e+dux2qteP3fz8Tt2uPfbz\nO/Xrx79/suPazTBO/Nyxrx1/7sleO/78Y1+r/czx++Nfa9u26f+tnE6uTNOsNgzjQWAeVg3Xu6Zp\nbnb2utJyGUYaMIS0NPuSq8igSEIDQnlsXglFvykiMjjSnkDkjNLS0upun4jrmKaVnBw5cuJWVFS/\nHd8uLj75dmwiZZeAgDRCQ4cQFATBwdTtjz0OCjrxOCjISmhO1j5235gtIODE49rkqLZ9bCLl52ff\nn5uvq6iuoKSypMFWWllq7StK69qlVQ2PneGSmy6mac4FurniWiLt2sH27dY6g48+ak8MhmEcHb3a\nyYHiA0quxOeUlUFenrXqQX5+w+OCghO3wsITtyo3La0ZEgJhYQ230FBrO/b4+C0kpH5/7HFwcP1r\nJzuuTZQWLaqf8kW8i2malFSWUFRRxJGKIxRVFNVtR8rr28WVxRRXFNcfVx49riiue6+ksoTiyuK6\nRMqONWK1tqB4naws6NDB+ksxOxuio+2J4+J3L+b7vd+z+K7FDOo4yJ4gpMWrrLTqD3NyrO3w4RO3\n3FwreTp2X1bm/HcHBVn1j5GRDbeIiPp97RYeXr+v3SIirGSptl2bSGkUpvmoMWsoqiiioKyA/LJ8\nCsoLKCwvpKDs6P5ou3Y7UnGk/rj8SN1rRRVF1Jg1bonR3/AnPCic0IBQwoPCCQsMIywwjNCAUGsf\nGEpowNEtsH4/aegkW2uuRFyqfXtrSoaFC2HmTLj3Xnvi0PqC4g6maY0eZWdb24EDcPDgiduhQ1Yy\nlZ/ftO8JDITYWGuLibG22Fjrl5VTbVFR9VtkpDXiI81fdU01+WX55JbmkluaS15ZHnmleQ2O88vy\nyStruM8vy6ewvNBlSVFIQAgRQRFEBkUSGRxJRFBEwy3Q2ocHhVv7wHDCg8JPug8LDKtLpIL8g5oU\nzyQmNflnUXIlXic1NY3x44ewcCFMnWpvcvX7hZB9tZIrb+YtNVemaY0i7d1rjb5mZcH+/fXHtclU\ndnbjiqv9/KBVK4iPt/bHb3Fx9fvY2Pp9WJhVlCv1vKWvuFNVTRWHSw6TU5LDoZJD5JTk1G2HSw6T\nU2od55bmcrjkMLmlueSX5WPS9DtL4YHhxITEEB0STXRwdN0+Kjiqbn/sFhkcae2PJlGRQVYiFegf\n6MI/CXspuRKvdNNN8MAD8N13Vv1V166ej6FteFtMNHIlVuJUUAC7dzfc9uyxkqm9e2HfPuupNEdE\nRVm1hW3bWvs2bRpu8fH1W2ysbqO1dCWVJWQXZZNdlM2BogMcKD7AweKDHCg6wMGSgxwsrt/ySvOa\nlCjFhMQQFxpHXGgcsSGxxIbGEhcSR2xobF07JiSG2BBrX7tFh0QT4KdU4nj6ExGvM2XKEMBKsN5/\n39qee87zcbSLaMcvh8I9Sq68mqtGIoqKrClAdu6EXbusfe3xrl3WU3FnEhNj1Qt26GDd3m7fHhIS\n6vcJCVZCFRrqkpClkbxp1Mo0TQrLC9l3ZB9ZR7LIOpLF/iP72V90dDt6fKDoAEcqHOh8RxkYtApt\nReuw1sSHx1v7sPi611qHtaZVWCtahbaiVVirumTK38/fjT9ty6PkSrzW+PH1ydXEiZ7/7b19ZHsA\n9hbu9ewXi9sUFMC2bbB1q7XfscMaGd2xw6pzOp2wMGv9y44d67ekJGurTagiIjzzc4h3M02T3NJc\n9hTuYU/BHvYW7mVPobXfW7iXfUf2sa9wH8WVxQ5dL9g/mHYR7Wgb0Za24W1pF9GONuFtaBPehrbh\nbeuO48PjiQuN00iSF9B/AfE6qalpTJkyhCFDrH/AMjMhPd3zj1CnxKbw+4Xw31Y7PfvF0ijH19HU\n1FgjTZs3W9uWLVYytXWrVTx+KsHB1sLhnTtDp07W8bH7uDjVMPk6V9Vc1Zg1ZBdlszNvJzvzd5KZ\nn8nugt1kFmSSWWAdl1SeeUKvsMAwOkR2oH1k+7otISKBhMgE2ke2p11EO9pFtCM6OFoTGfsYJVfi\ntfz8YOxY+MMfrMJ2TydXnWI6AbArfxc1Zg1+hgpfvIm1Zhl8/721bdgAmzbBTz9Zs3efTEgInH22\ntXXtaq0GULvv0EG1TVKvqKKIjLwMduTuYEfeDnbk7mBnfn0yVV59+gK7yKBIOkZ3JDEqkaSoJJKi\nk0iMSiQxKpEOkR3oENVBSVMzpnmuxKtt3Qrdulm3W7KzrblyPKntq205WHyQPY/sITEq0bNfLnWK\ni2HdOmtJpLVrreMNG6w6qZNJSIDu3eu3bt2shCoxUQmU1CuuKGZ77na25W5j2+FtbM3dyrbD29ie\nu50DxacZ5gTiw+JJiU2hU0wnkqOTrS0mmY7RHUmOTiY6xKYJ+sRlDMPetQVF3Obss2HgQGtk4uOP\nrZEsT+oc25mDxQfZmbdTyZWHFBTA6tWwahX88IOVTG3devKFa9u2hZ49oVcvOPdc6/icc6zCchGw\n6p+yi7LZnLOZLTlbGmx7Cvec8nNB/kF0ju1Ml9gu1hbXhc6xnUmJSSE5JpmIIBXYyakpuRKvU1tz\nVWv8eCu5mjrV88nVo18XcfO5kJGXweDkwZ798hagrMxKpJYts5KpVausQvPjBQRYyVPfvtbWu7eV\nUMXHt4y5i8QxB4oOsP7gejYe3MimQ5vYeGgjGw9tJL/s6EysO4GU+vMD/QLpHNuZs1udzVlxZ3FW\nq7M4u9XZdI3rSmJUokoBpMmUXInXu/VWmDABFiyw5hVKSvLcd8eGxgJWciXOMU3r4YSlS61katky\na1SqsrLhecHB0KePtWj3+efDeedZt/Y0W7jUKqsqY+PBjaw7sI51B9ax/uB61h1Yx6GSQyc9PyYk\nhh7xPYj2j2bo0KGc0/oczml9DimxKXqyTtxCvUq8zrGjVmDd4hk1yloKZ9o0eOopz8Wy6+FU+Gwx\nO/P1xGBjVVdbtVHffQdLlljbvn0NzzEMawTqoovgggus7dxzraVbHKVRq+atoKyANdlrWLN/DWuy\n17A2ey2bDm2i2qw+4dyo4Ch6telFzzY96RHfg3Pjz6VHfA/aRbRT4bh4lJIr8Qnjx1vJ1dSp8Nvf\neu6R+JQY6x6CRq7OrKrKGolauBDS0mDx4hMn3oyLg0susZKp2oQqMtKWcMULFVUUsWb/GlZlrWLV\n/lWs3LeSbbkn3ic2MDin9Tn0bdeXXm160bttb3q37U1SVJKSKPEKSq7E6xxfcwUwbJi1TMjWrbB8\nufUPsyf0e2s2hCq5OhnTtJ7Y++YbK6FKT4fCwobndO4MgwbVb926uf5pPdVc+aYas4YtOVtYtndZ\n3bbx0MYTFgEO8g+id9ve9GvXz9oS+tGrTS/Cgxr/6LD6iniKkivxCQEBcOed8Oqr8NJL8Mknnvne\nqOAoAvwC2F+0n9LKUkIDW/a6JVlZMH++lVDNn29Nj3GsLl2s+ciGDoXLLrOmPhABa1Rq+d7lfLfn\nO5bsXsLyfcspLG+YjQf4BdCnbR/6t+9ft/Vs05Mg/yCbohZpGs1zJT4jK8uamqG42Bop8dQvoF3f\n6MqOvB1s+tUmusd398yXeonKSutJzS+/tLb16xu+n5AAV10FV1xhJVSefNhAvNvhksOkZ6aTnpnO\nkj1LWLN/zQl1UklRSVyUeBEDOgzgosSLOC/hvBb/C4x4D81zJS1C+/bwxBPw7LPwf/8HK1eCvwfW\nGu0c25kdeTvIyMtoEcnVgQP1ydS8eQ1v9YWHWyNSV11lbT16aEkYseSW5pKemc7CnQtJy0xj3YF1\nDd73N/w5P+F8BnUcxKCOgxiYOJAOUR1silbEvZRcidc5Wc1Vrd/8Bt5+G9assZ4cTE11czATJ9K5\nf2eAZvvEoGlaa/B99pm1LVvWcMLOc86Ba66xtkGDvG9KBNXR2KOsqozvdn/HNxnf8E3GN6zZvwaT\n+o4T7B/MwKSBXJZ8GYM7DmZA4gDbJ95UXxFPUXIlPiUsDP74R2sy0aeegptvtpbGcafm+MRgTY2V\nRH30EcyeDTt21L8XHAyXXw7XXgtXX20tXCximiabDm3iq+1fMW/HPBbvXkxZVVnd+7XJ1JDkIQzp\nNIQBiQMICQixMWIR+6jmSnxOTY31tODKldYtwkmT3Pt9H278kNGzRjOq2yg+ve1T936ZG1VVwaJF\n1jJCn3wC+/fXv9eqlZVMXXed9WSmuxNW8Q1FFUUs2LmAL7d9yVfbv2J3we4G7/dt15erOl/FlZ2v\nZFDHQYQFhtkUqYjrqeZKWhQ/P3jtNRg8GF55Be69171PpXWO9d3bglVVVvH/jBlWQpWbW/9ecjLc\neCPccANcfLFn6tfE++0t3MvnP33O7J9ms3DXQiqqK+reaxPehhFdRzCiywiu6HwFbcLb2BipiPdS\nciVe53Q1V7UGDYJbboEPP7RuD77/vpuCmTiRlCd+DVi3BU3T9PpJCqurrRGqGTOsUaqcnPr3unWD\nm26ytn79mkcxuuponGOaJhsObuDTLZ8y+6fZ/LD/h7r3DAwuSryIa7pew9VnXc15Cef59Hp76ivi\nKUquxGf96U9WvdC0afDQQ9Zs3+4QGxJLdHA0BeUF5JTkEB8e754vcoJpWjVUH3xgzWR/8GD9e926\nWeszjh5tLS0jYpoma7PX8uGmD5m1aVaDWdBDA0IZ3nU4o7qNYuRZI72yv4t4O9VciU974gl4+WVr\nSZX0dNfP/l3rvLfOY032Gpb/fDkXdrjQPV/SBOvXw3//a227dtW/3rVrfULVq1fzGKES55imyZrs\nNczcOJNZm2axI6/+KYbWYa0Z1W0Uo7qN4srOV2quKRFUcyUt2FNPweTJ1uLAjz4Kf/6zexKJlNgU\n1mSvISMvw/bkas8ea4Rq+nRr+ZlaHTrAbbfB7bfDeecpoRLLjtwdfLD+Az7Y8AFbcrbUvd4mvA03\ndb+Jm3vczKXJlxLgp38ORFxF/zeJ13Gk5qpWdLQ1anP11fCXv1gzhj/2mAuDmTjRmusqxipqt2s6\nhvx8mDXLSqgWLap/PTbWqj0bM8Yq8HfXyJ03Ux3NiQ4VH+J/G/7HBxs+YNneZXWvx4fFc0uPW7jl\n3FsY3HEw/n4t6ykG9RXxFCVX4vOuuMIqaL/9dnj8cWjbFsaNc+13pMRakz3tzPPcE4Pl5fDFF/Cf\n/8CcOVBx9KGtkBBryoQ77oARIyBIy64JUFVTxdztc3lvzXvM2TqHyppKAMIDw7mh+w2M6TmGKztf\nSaB/oM2RijR/qrmSZuONN2DCBGtKgc8/t0azXGXu9rlc/Z+ruTzlcr4d963rLnycmhpYssQaofrw\nQ2vECqxbfJdfbi1efeONEBXlthDEx2zJ2cJ7a95j2rppZBdZK2n7GX4M7zKcsb3Hcl236wgPCrc5\nShHfo5orEeDXv7YmxvzTn6yZ2xcsgAEDXHPt2rmu3HFb0DThxx/rC9P37Kl/r18/a4TqttusmioR\ngPKqcj7e/DH/+uFfpGem171+dquzuavvXYzrM472ke1tjFCkZdPIlXidxtRcHc804e67YcoUa9bx\nJUustfGa7GjNVXlVOaF/CMXP8KP06VKX3FrZvr0+odq8uf71jh2thOqOOzR1giNaUh1NRl4Gb//w\nNu+teY9DJYcA67bfbT1v4+5+dzMwcaDXz8Nmp5bUV8R5GrkSOcowrIWdDx2y6pX694c//AEefNC5\nGciDA4LpENWBvYV72VO4p24kq7G2bbMK0z/6CH6on6uR1q2twvTbb7emlWiJhelycqZp8k3GN7y+\n7HW+2v5V3eu92vTi/v73c0fvO4gK1n1iEW/i1MiVYRg3AxOB7sAFpmmuPs25GrkSjykpgdRUq24J\nrAlG33kH+vRp+jUvm3IZ6ZnpfDP2G67sfKVDnzFN2LTJSqY++gjWrat/LyLCWnrm9tvhyishUHXG\ncoySyhKm/TiNN1a8waZDmwBrceRbe97KL8//JRclXqRRKhE3snPkaj1wA/CWk9cRcamwMGum8s8+\ngwcesBZ5Pv98ay6sZ5+13m+slJgU0jPTz/jEYFGRVe/11Vcwd27DyT2jomDUKGv5mWHDIFRzNcpx\nso5k8bflf+Pt1W+TW2otBtk+sj0PXvAg955/L63DWtscoYiciVPJlWmaPwEY+vVJXMiZmqvjXXcd\nDB0KzzwDf/sbvPSSlXSNHw/XXmsVjJ/2FtzRmis4dVF7YSGsWgXLl8P8+bB4MVRW1r/fqpUVx803\nW9NGBAe75EeTo5pLHc1POT/xytJXmLZuWt1iyRd2uJCHBzzMzT1u1hQKLtBc+op4P9VcSbMXGQl/\n/atVIH7vvdatudqcKSEBrrnGSrQGDrQm5TzZvFGmCW2DrOTq+y0Z/HMbrFhhJVRbtljv1/Lzs641\nYoQ1HcR55zlX7yXN28p9K/nTd3/ik82fYGJiYHBT95t49OJHuSjxIrvDE5EmOGPNlWEY3wBtj30J\nMIGnTdP8/Og5C4HfqOZKvF1lJXz9tVXsPmcO7N174jmhodbM79HRVmKWlwdZWVDa+ju4ZxDsuwDe\nWVF3flAQ9O0LF15oFaNfdZU1WiVyOumZ6Ty36Dm+3WnNmxbkH8T4PuN59OJHObvV2TZHJyJurbky\nTfOqplz4ZFKv60un864HIKZ4GX3PaceQe6YAkPZuKoDaaru9fe21EHEgldH9Ia7/FL74AvJWp1JY\nAO+lTaG0FEb3t86futj6/PjBqQSElvAuENgmg2fvSqVVKzj/lin06QPfTz96/dvs//nU9u52emY6\nX749muwjB/i2FUQGRfJsdUd6xPfgmp+9bXt8aqvdUtt0SSUtLY1dqz/FWS6Z5+royNWjpmn+cJpz\nNHIlDklNdV3NVWOZpvWkYUGBNTt61GsTKX5sIu3bQ0SESdiLYZRVlVHwZIEef/cSvlJHk56ZzqRF\nk1iwcwEA0cHRPHLRI0y4aAIxITE2R9cy+EpfEe9g29OChmFcD/wNaA3MMQxjrWmaLlx0RMSzDAPC\nw62tfXsgEehW9y4pMSlsztnMzryd9GnnxLwO0mL8kPUDT377JPMz5gNKqkRaAs3QLtII135wLV9s\n+4KPR3/MDd1vsDsc8WLbDm/jmYXPMHPjTACigqN45KJHePiih5VUifgAzdAu4iHuXGNQmof9R/bz\n3KLn+Peaf1NVU0WwfzAPXfgQTw56klZhetJBpCXQIhvidVJT0+wOod7ROa5qpcSkALAz//QTiYrn\npKWl2R0CYM2oPiltEl3/1pV//fAvaswa7u57N9se2sYrw15RYuUFvKWvSPOnkSuRRtDIlRzPNE3+\nu+G/PDH/CfYWWnN7XH/O9fzh8j/QI76HzdGJiB1UcyXSCOsPrKf3v3rTrVU3tjy4xe5wxGYr9q1g\nwtwJLNu7DIB+7frx+ojXuTT5UpsjExFnqeZKxENSYq3bgrvyd1Fj1uBn6M56S3Sg6ACPz3+c9398\nH4C24W158YoXGd9nPP5+mo5fpKXTvwzidby55ioiKIL4sHjKq8vZf2S/PTFJA56so6muqeYfK/9B\nt7934/0f3yfIP4gnL3mSrQ9t5e5+dyux8nKquRJP0ciVSCN1ju3MoZJDZORl0CGqg93hiIesylrF\n/V/cz6qsVQCM6DqCv1/9d7rEdbE5MhHxNqq5EmmksZ+MZfq66bx17Vv84vxf2B2OuFl+WT5Pf/s0\n/1z1T0xMEqMS+euIv3LDOTdgGE0qxxARH+BMzZVuC4o0Ut+2fQFYs3+NzZGIu83eMpseb/bgH6v+\ngZ/hx6MDH2XzA5u5sfuNSqxE5JSUXInX8eaaK4B+Cf0AWJOt5MobuKOO5mDxQW6bdRvXz7ie/UX7\nGZg4kDX3reGVYa8QERTh8u8Tz1DNlXiKaq5EGqlvO2vkat2BdVTXVKuIuRkxTZMP1n/AhLkTOFx6\nmLDAMP54xR954IIH9N9ZRBymmiuRJkh+PZndBbvZ9KtNdI/vbnc44gJZR7K4b859zNk6B4ArO1/J\n29e+XTf9hoi0LKq5EvGwfu10a7A5+XDjh/T6Zy/mbJ1DdHA07173LvPunKfESkSaRMmVeB1vr7mC\n+luDa7PXejAYORln6mjyy/K58+M7GT1rNLmluYzoOoJND2zi7n53q2C9GVLNlXiKaq5EmkAjV77v\n24xvSZ2dyt7CvYQFhvHnYX/mvvPvU1IlIk5TzZVIE+wu2E3y68m0Cm3FoccO6R9kH1JeVc6T85/k\n9eWvAzCgwwCm3TCNs1qdZXNkIuJNtLagiIclRSURGxLL4dLD7Duyj8SoRLtDEgdsO7yNW2fdyprs\nNQT4BfDspc/y28G/JcBPfxWKiOuo5kq8ji/UXBmGUT/flSYTtZWjdTTT103nvLfPY032GjrHdua7\nu7/jd5f9TolVC6KaK/EUJVciTaS6K99QVFFE6qepjP1kLEUVRdx67q2s/sVqLuxwod2hiUgzpZor\nkSaavm46Yz8Zyw3n3MDHt35sdzhyEusOrGP0h6P56fBPhAaE8sbVb3BPv3tUIyciZ6SaKxEbaOTK\nu037cRr3zbmP0qpSzo0/lxk3z+DcNufaHZaItAC6LShexxdqrgC6te5GSEAIu/J3kVea57mYpIHj\n62jKq8r51Re/Ytyn4yitKiW1byor7l2hxEpUcyUeo+RKpIkC/ALo1aYXAD8e+NHmaARgT8EeLpty\nGf9c9U+C/IN469q3eO+69wgLDLM7NBFpQVRzJeKE+z6/j7dXv81rw17jkYGP2B1Oi/Ztxrfc9tFt\n5JTk0DG6I7NumcUFHS6wOywR8VFaW1DEJnXTMajuyjamafLa968xbPowckpyGNZlGD/84gclViJi\nGyVX4nV8peYKtMag3cqqyrj6D1fzm3m/ocas4ZnBz/DlmC9pHdba7tDEC6nmSjxFTwuKOKF32974\nGX5sOrSJsqoyQgJC7A6pxdh/ZD83zryRZduXEXZ2GFOvn8rNPW62OywREdVciTirx5s92JyzmZX3\nrqR/+/52h9MirMpaxfX/u559R/bRMbojs2+bXTeKKCLiCqq5ErGRbg161owNMxg8eTD7juzjkqRL\nWHnvSiVWIuJVlFyJ1/Glmis4ZjJRrTHoVqZp8of0P3DbR7dRVlXGPf3uYcH4BWxaucnu0MRHqOZK\nPEU1VyJO0hOD7ldRXcF9c+5jytopGBi8Nvw1JgyYoGVsRMQrqeZKxEk5JTnEvxJPeGA4BU8W4O/n\nb3dIzUpeaR43zbyJhbsWEhYYxgc3fsCoc0bZHZaINHOquRKxUeuw1iRGJVJcWcz23O12h9Os7Mzb\nycXvXczCXQtpF9GORamLlFiJiNdzKrkyDONlwzA2G4ax1jCMjwzDiHJVYNJy+VrNFWgRZ3dYsW8F\nA/49gC05W+jZpifLf778pE9jqo5GHKW+Ip7i7MjVPOBc0zT7AtuA3zofkojv0RODrvXlti8ZOnUo\nh0oOMazLMJbctYSO0R3tDktExCEuq7kyDON64CbTNMee4n3VXEmz9cnmT7hx5o0M6zKMr+/82u5w\nfNrkNZO59/N7qTarSe2bytvXvk2gf6DdYYlIC+MtNVd3A1+58HoiPqP2icHV+1dTY9bYHI1vqp1q\n4e7P7qbarObpwU/z3nXvKbESEZ9zxuTKMIxvDMNYd8y2/uj+Z8ec8zRQaZrmB26NVloEX6y5So5O\npkNkB3JKcvgx+0f3xtQMVddU8+CXD/LMwmcwMPj71X/nhctfcGiqBdXRiKPUV8RTzjjPlWmaV53u\nfcMwUoFrgMvPdK3U1FQ6deoEQExMDH379mXIkCFAfadXW22vaoND5y9atIg+pX3Yxz6+2v4VBT8V\neEf8PtAuqypj2PPDWJy5mOCuwUy/cTqtD7YmLS3NK+JTu/m0a3lLPGp7V7v2eNeuXTjLqZorwzBG\nAH8GLjVN8/AZzlXNlTRrH236iJs/vJnBHQeTfle63eH4hCPlR7h+xvUs2LmA6OBoPrv9My5NvtTu\nsEREnKq5cja52gYEAbWJ1TLTNH91inOVXEmzVlBWQKuXWwGQ83gOMSExNkfk3XJKcrjmP9ewMmsl\n7SLa8fWdX9O7bW+7wxIRAWwsaDdN8yzTNJNN0zzv6HbSxEqkMXyx5gogOiSai5MuptqsZn7GfPfF\n1AzsLdzLpZMvZWXWSlJiUlhy15ImJ1bH3/IRORX1FfEUp5IrEWno6q5XAzB3+1ybI/Fe2w5vY9B7\ng24BJj0AAA28SURBVNics5mebXqy5O4ldInrYndYIiIuo7UFRVxobfZa+r3Vjw6RHdjzyB4tLHyc\ntdlrGT59OAeLD3JR4kV8MeYL4kLj7A5LROQE3jLPlUiL16dtHxIiEth3ZB/rD663Oxyv8v2e7xky\nZQgHiw8yrMsw5o+dr8RKRJolJVfidXy15gqs33RGdB0BwFfbNKdurQU7F3DVtKsoKC/gxu438tlt\nnxEeFO6Sa6uORhylviKeouRKxMVqk6u5O1R3BfDF1i+45j/XUFxZzLg+45hx8wyCA4LtDktExG1U\ncyXiYnmlebR+pTV+hh+HHz9MVHCU3SHZZubGmdzx8R1U1VRxf//7+fs1f8fP0O90IuL9VHMl4kVi\nQ2MZmDiQqpoqvs341u5wbDN5zWRu/+h2qmqqeOzix3jzmjeVWIlIi6C/6cTr+HLNVa3aKRm+2t4y\n667eXPEmd392NzVmDc8PfZ6XrnzJbU9Oqo5GHKW+Ip6i5ErEDerqrrbPpaXdDn/5u5d58KsHAXht\n2Gs8c+kzmpJCRFoU1VyJuEGNWUPCnxM4WHyQDfdv4Nw259odktuZpsmkRZOYtGgSBgb/HPlP7ut/\nn91hiYg0iWquRLyMn+FXPyVDC7g1aJomT8x/gkmLJuFn+DH1+qlKrESkxVJyJV6nOdRcQcupu6ox\na3joq4d4ZekrBPgF8L+b/sfYPmM99v2qoxFHqa+IpwTYHYBIc3VV56vwM/xYnLmYoooiIoIi7A7J\n5aprqvnF57/gvbXvEeQfxKxbZvGzbj+zOywREVup5krEjQa+O5Ble5cx+7bZXNftOrvDcanK6krG\nfjKWGRtnEBoQyuzbZnNVl6vsDktExCVUcyXipWpvDc7aNMvmSFyrrKqMm2bexIyNM4gMiuTrO79W\nYiUicpSSK/E6zaXmCmBs77H4GX78b8P/yDqS5ZqYbFZUUcS1H1zL51s/Jy40jgXjFzA4ebBt8aiO\nRhylviKeouRKxI1SYlO4sfuNVNZU8rflf7M7HKfll+UzfPpwvt35LW3D27IodRH92/e3OywREa+i\nmisRN1u2dxkD3x1ITEgMex7Z47OF7TklOQyfPpzV+1eTFJXEt+O+5axWZ9kdloiIW6jmSsSLXZR4\nERcnXUx+WT6T10y2O5wm2VOwh8GTB7N6/2q6xnVl8V2LlViJiJyCkivxOs2p5qrWbwb+BoDXl79O\ndU21S67pKVtytnDJe5ewJWcLPdv0JD01neSYZLvDqqM6GnGU+op4ipIrEQ8Y1W0UXWK7kJGXwadb\nPrU7HIet3LeSQe8NYk/hHi5Oupj01HQSIhPsDktExKup5krEQ95c8SYPfvUgFyVexPf3fG93OGc0\nP2M+1//veoori7nmrGv48JYPCQsMszssERGPUM2ViA9I7ZtKbEgsy/YuY+mepXaHc1qzNs1i5Acj\nKa4s5s7ed/LprZ8qsRIRcZCSK/E6zbHmCiA8KJz7+98PwJ+//7PLrutKpmnyxvI3GP3haCqqK5gw\nYAJTr59KoH+g3aGdkupoxFHqK+IpSq5EPOjBCx8kyD+ITzZ/wo7cHXaH00BVTRUPfvkgE+ZOwMTk\nD5f/gb8M/wt+hv6aEBFpDNVciXjYXbPvYsraKTx4wYP87RrvmFi0sLyQW2fdytztcwnyD2LKqCnc\n3ut2u8MSEbGNMzVXSq5EPGz9gfX0/ldvwgLD2PHrHbSLaGdrPJn5mVz732vZcHADrcNa8+mtn3JJ\nx0tsjUlExG4qaJdmpbnWXNXq1bYXI88aSUllCbfOupXK6kqXf4ejVuxbwYB/D2DDwQ2c0/oclv98\nuc8lVqqjEUepr4inKLkSscE7P3uHhIgE0jPTefybxz3+/bWF64MnD+ZA8QGuSLmC7+/5ns6xnT0e\ni4hIc6PbgiI2WbpnKUOmDKGyppL/3PgfxvQa45HvzS3N5Z7P7qmbzPRX/X/F6yNe9+onAkVEPE01\nVyI+6h8r/8EDXz5AaEAoy36+jN5te7v1+77f8z23fXQbuwt2Ex0czbvXvctNPW5y63eKiPgi1VxJ\ns9Lca66OdX//+xnfZzylVaXcMOMG8krz3PI9NWYNLy15icGTB7O7YDcXdriQNfetaRaJlepoxFHq\nK+IpSq5EbGQYBv8c+U/6tetHRl4Gd3x8BzVmjUu/Y+mepQx8dyBPfvsk1WY1jw58lMV3LSYlNsWl\n3yMiIhanbgsahvEcMAqoAQ4AqaZpZp/iXN0WFDmFXfm76P92fw6XHuZX/X/Fy1e9THhQuNPXfGL+\nE8zcOBOAhIgE3vnZO4w8e6QrQhYRadZsq7kyDCPCNM2io8cPAT1M07z/FOcquRI5jfkZ8xk+fTg1\nZg0dIjvw4hUvcmfvOxs9Q3pheSF/XPxH/rLsL5RXlxMSEMJjFz/G45c8TkRQhJuiFxFpXmyruapN\nrI4KxxrBEnFKS6q5OtaVna8kbXwa5yWcx74j+xj/6XgG/HsAizMXn/GzJZUlzN4ym7tn303KX1P4\n03d/ory6nDt63cHWB7fy3NDnmm1ipToacZT6inhKgLMXMAzjBWAckA8MdToikRZscPJgVt67kmk/\nTuOpBU+xKmsVl065lOvPuZ4L219IVHBUg21X/i5m/zSbeTvmUVpVWnedi5Mu5rVhrzEgcYCNP42I\nSMt0xtuChmF8A7Q99iXABJ42TfPzY857Agg1TXPiKa6j24IijVBcUcwrS1/h5f9v735C7DyrOI5/\nf1K7UXBoLC4ybQPW2IrK0MUYEXFABJWSdCEqCKauVPAPKEUogoVuxI34ZyFCtVEQFy7SghWqSBBd\nlFI7NGiLLXXSOpK40GGoI23U4+LeJDe38+femcd33pv5fmCY+/C+3DmLw7mH5z3z3N9946rGaSuL\nhxc58dYT3HXbXdz+xttJdrWbLUmiJ+dcJbkJeKSq3rHF9Tp58iRHjhwBYG5ujoWFBZaWloAr27Wu\nXbu+er26vsr9p+5n/ZV15m6bY/3ldZ5/8nk2Lm4w/8557jx6J4cuHOLG193Yi3hdu3btehbXl16v\nrKwAcOrUqX0baL+1qp4bvv488N6q+ugW97pzpYncffcZHnxwab/DGLjvvk7nrjS9M2fOXC6S0nbM\nFU1jLztXe525+nqSowwG2c8Bn9nj+0mSJM00v/5GkiRpjF9/I0mS1BM2V+qdg3rOlXZndBhV2o65\noq7YXEmSJDXkzJUkSdIYZ64kSZJ6wuZKvePMlabhHI0mZa6oKzZXkiRJDTlzJUmSNMaZK0mSpJ6w\nuVLvOHOlaThHo0mZK+qKzZUkSVJDzlxJkiSNceZKkiSpJ2yu1DvOXGkaztFoUuaKumJzJUmS1JAz\nV5IkSWOcuZIkSeoJmyv1jjNXmoZzNJqUuaKu2FxJkiQ15MyVJEnSGGeuJEmSesLmSr3jzJWm4RyN\nJmWuqCs2V5IkSQ05cyVJkjTGmStJkqSesLlS7zhzpWk4R6NJmSvqis2VJElSQ85cSZIkjXHmSpIk\nqSdsrtQ7zlxpGs7RaFLmirpicyVJktSQM1eSJEljnLmSJEnqiSbNVZIvJ/lvkhtavJ8ONmeuNA3n\naDQpc0Vd2XNzlWQe+ABwbu/hSHD+/PJ+h6AZsrxsvmgy5oq60mLn6pvAPQ3eRwLg2LG1/Q7hCneu\nem9trUf5ol4zV9SVPTVXSY4DL1bV2UbxSJIkzbQdm6skv0zy1MjP2eHv48C9wNdGb/+/RaoD4/Tp\nlf0O4Qp3rnpvZWVlv0PQjDBX1JVdH8WQ5O3Ar4ANBk3VPLAKLFbV3za533MYJEnSzNjtUQzNzrlK\n8mfgjqr6R5M3lCRJmkEtz7kqfCwoSZIOuM5OaJckSToImp/QnuSDSZ5J8qckX9nk+vuSrCX5/fDn\nq61j0GxI8kCSC0me2uaebyd5NslykoUu41N/7JQr1hVdkmQ+ya+T/GH4D1hf2OI+a4smypfd1Jfr\nGgf5GuC7wPuBvwKPJ3moqp4Zu/U3VXW85d/WTPoh8B3gR5tdTPIh4M1V9ZYk7wK+BxzrMD71x7a5\nMmRdEcC/gS9V1XKS1wNPJHl09HPI2qIRO+bL0FT1pfXO1SLwbFWdq6qLwE+BE5vc52yWqKrfAtv9\nA8QJhh+mVfUY8IYkb+oiNvXLBLkC1hUBVXW+qpaHr18CngYOj91mbREwcb7AlPWldXN1GHhxZP0X\nNg/y3cOt2J8neVvjGHTtGM+nVTbPJwmsKxqT5AiwADw2dsnaolfZJl9gyvrS9LHghJ4Abq6qjeHW\n7Gng6D7EIenaYV3RVYaPeH4GfHG4IyFtaYd8mbq+tN65WgVuHllfOlj0sqp6qao2hq9/Abw2yQ2N\n49C1YRW4aWT9qnySwLqiqyW5jsEH5Y+r6qFNbrG26LKd8mU39aV1c/U4cGuSW5JcD3wceHj0htHn\n2kkWGRwH8ffGcWh2hK2fZT8MfBIgyTFgraoudBWYemfLXLGuaMwPgD9W1be2uG5t0aht82U39aXp\nY8Gq+k+SzwGPMmjcHqiqp5N8enC5vg98JMlngYvAv4CPtYxBsyPJT4Al4FCSFxh8T+X1DHOlqh5J\n8uEkzwH/BD61f9FqP+2UK1hXNJTkPcAngLNJnmRwwPW9wC1YWzRmknxhF/XFQ0QlSZIaan6IqCRJ\n0kFmcyVJktSQzZUkSVJDNleSJEkN2VxJkiQ1ZHMlSZLUkM2VJElSQzZXkiRJDf0PgAUXDe4dKYoA\nAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x2e0c2331da0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def lennard_jones_potential(r, epsion=1.0, sigma=1.0):\n",
    "    sigma12 = np.power(sigma, 12.0) \n",
    "    sigma6 = np.power(sigma, 6.0)\n",
    "    return 4.0 * epsion * (sigma12 * np.power(r, -12.) - sigma6 * np.power(r, -6.))\n",
    "\n",
    "def lennard_jones_force(r, epsilon=1.0, sigma=1.0):\n",
    "    sigma13 = np.power(sigma, 12.0) \n",
    "    sigma7 = np.power(sigma, 6.0)\n",
    "    return 24.0 * epsilon / sigma * (2.0 * sigma13 * np.power(r, -12.) - sigma7 * np.power(r, -6.)) \n",
    "\n",
    "def plot_lj_pot(r, epsilon=1.0, sigma=1.0):\n",
    "    pot = lennard_jones_potential(r, epsilon, sigma)\n",
    "    plt.plot(r, pot, lw=2, label='LJ potential')\n",
    "\n",
    "def plot_lj_force(r, epsilon=1.0, sigma=1.0):\n",
    "    force = lennard_jones_force(r, epsilon, sigma)\n",
    "    plt.plot(r, force, lw=2, label='LJ force')\n",
    "\n",
    "r = np.linspace(0.5, 2.5, 100)\n",
    "def plot_lj(epsilon=1.0, sigma=1.0):\n",
    "    fig = plt.figure(figsize=(10, 6))\n",
    "    plot_lj_pot(r, epsilon=epsilon, sigma=sigma)\n",
    "    plot_lj_force(r, epsilon=epsilon, sigma=sigma)\n",
    "    plt.ylim([-4.0, 4.0])\n",
    "    plt.axvline(sigma, linestyle=':', color='blue', label='V = 0')\n",
    "    plt.axvline(np.power(2.0, 1./6.) * sigma, linestyle=':', color='red', label='$r^*$')\n",
    "    plt.axhline(-epsilon, linestyle=':', color='orange', label='$V_\\mathrm{min}$')\n",
    "    plt.legend()\n",
    "    plt.title('')\n",
    "    plt.grid()    \n",
    "    \n",
    "# plot_lj_pot(r, 1.0, 1.0)\n",
    "interact(plot_lj, epsilon=(0.1, 2.0), sigma=(0.1, 2.0));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (root)",
   "language": "python",
   "name": "root"
  },
  "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}