Ziegler-Biersack-Littmark (ZBL) potential

zbl.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ziegler-Biersack-Littmark (ZBL) potential\n",
    "\n",
    "This notebook shows a plot of ZBL potential when inner radius cutoff is set zero.\n",
    "\n",
    "\n",
    "### References\n",
    "- https://docs.lammps.org/pair_zbl.html\n",
    "- https://docs.lammps.org/pair_gromacs.html\n",
    "- https://github.com/lammps/lammps/blob/develop/src/pair_zbl.cpp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_context('poster')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class ZBL:\n",
    "    p = 0.23;\n",
    "    # angstrom\n",
    "    a0 = 0.46850;\n",
    "    c = [0.02817, 0.28022, 0.50986, 0.18175]\n",
    "    d = [0.20162, 0.40290, 0.94229, 3.19980]\n",
    "        \n",
    "    def __init__(self, zi, zj, rc):\n",
    "        self.zi = zi\n",
    "        self.zj = zj\n",
    "        self.rc = rc\n",
    "        \n",
    "        # e * e / (4 * pi * epsilon_0) / electron_volt / angstrom\n",
    "        self.zzeij = 14.399645478425668 * self.zi * self.zj  # eV.angstrom\n",
    "        self.a = self.a0 / (zi ** self.p + zj ** self.p)\n",
    "        self.da = [dl / self.a for dl in self.d]\n",
    "        \n",
    "        # switching function\n",
    "        ec = self._e_zbl(self.rc)\n",
    "        dec = self._dedr(self.rc)\n",
    "        d2ec = self._d2edr2(self.rc)\n",
    "        # coefficients are determined such that E(rc) = 0, E'(rc) = 0, and E''(rc) = 0\n",
    "        self.A = (-3 * dec + self.rc * d2ec) / (rc ** 2)\n",
    "        self.B = (2 * dec - self.rc * d2ec) / (rc ** 3)\n",
    "        self.C = -ec + rc * dec / 2 - rc * rc * d2ec / 12\n",
    "        \n",
    "    def _phi(self, r):\n",
    "        phi = np.sum([cl * np.exp(-dal * r) for cl, dal in zip(self.c, self.da)])\n",
    "        return phi\n",
    "    \n",
    "    def _dphi(self, r):\n",
    "        dphi = np.sum([-cl * dal * np.exp(-dal * r) for cl, dal in zip(self.c, self.da)])\n",
    "        return dphi\n",
    "        \n",
    "    def _e_zbl(self, r):\n",
    "        phi = self._phi(r)\n",
    "        ret = self.zzeij / r * phi\n",
    "        return ret\n",
    "    \n",
    "    def _dedr(self, r):\n",
    "        phi = self._phi(r)\n",
    "        dphi = self._dphi(r)\n",
    "        ret = self.zzeij / r * (-phi / r + dphi)\n",
    "        return ret\n",
    "    \n",
    "    def _d2edr2(self, r):\n",
    "        phi = self._phi(r)\n",
    "        dphi = self._dphi(r)\n",
    "        d2phi = np.sum([cl * (dal ** 2) * np.exp(-dal * r) for cl, dal in zip(self.c, self.da)])\n",
    "        ret = self.zzeij / r * (d2phi - 2 / r * dphi + 2 * phi / (r ** 2))\n",
    "        return ret\n",
    "    \n",
    "    def energy(self, r):\n",
    "        if r > self.rc:\n",
    "            return 0\n",
    "        \n",
    "        e = self._e_zbl(r)\n",
    "        e += self.A / 3 * (r ** 3) + self.B / 4 * (r ** 4) + self.C\n",
    "        return e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "rmin = 0.1\n",
    "rc = 1.0\n",
    "zbl = ZBL(10, 10, rc=rc)  # Ne-Ne\n",
    "\n",
    "rs = np.linspace(rmin, rc, num=100, endpoint=True)\n",
    "energies = [zbl.energy(r) for r in rs]\n",
    "screenings = [e / zbl.zzeij * r for e, r in zip(energies, rs)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0001, 10000.0)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAekAAAHNCAYAAAA684O3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3debgcVZ3/8fcnQCCQBBKvhCUEojcBQZAtBOYqoBkc9IcLYYkLcXDEQUSJ6wAzos6gEsaNRWRwQCIgiCwyrihGAhJDCHtYYhJMCESIxEBIkLDl+/ujqk2n6e7be9ft/ryep57K7Tqn6nTdzv32OXUWRQRmZmaWPYPaXQAzMzMrzkHazMwsoxykzczMMspB2szMLKMcpM3MzDLKQdrMzCyjHKQbQNImku6WFJKObnd5zMysMzhIN8ZJwA7tLoSZmXUWB+k6SRoFnAmc1u6ymJlZZ3GQrt/XgV8Ds9pcDjMz6zAdFaQl7SppmqQrJC2QtL7S58SSPiDp95JWS1or6U5JJ0sqeY8kHQxMBj7fyPdhZmYGsGm7C9BgJwHTqs0k6QLg48A6YCbwEjAJ+A4wSdIxEfFKQZ5NgQuAsyLiMUm71Fd0MzOzjXVUTRp4gKT5eQrQC9zSXwZJR5EE6CeBvSLiiIg4EhgHPAwcCXyiSNZpwBDgG40pupmZ2cY6qiYdERfn/yypkmynp/tTI2JR3rlWSDqJ5FnzaZLOj4j16Xl7gC+TBPchkoYAw9OsW0raOiJW1/NezMzMOipIV0vSaGA/4EXgmsLjEXGLpOXAjsCBwB/SQ6OBocBlRU77A2A1sE0zymxmZt2jq4M0sE+6fzAini+RZh5JkN6HDUF6MfDWgnTbAVeR1LBnFjuRpOOB4yss237AJsCq9HpmZtYZekkqeksiYp9yCbs9SI9N94+WSbOsIC0RsZaCIVd5HccejIjbSpxrF+CQKsu4Y7qZmVlnGdtfgm4P0kPT/XNl0qxN98MacL2lVNCZLXUQMHjrrbdm7733bsClq7P6+ZdY/vTzPLvupVcdG77FZuw4YghbD9ms5eUyMxvo7r33XlavXg0b4ktJ3R6kcz3Lot4TRcTSvPOVSjMDmFHJ+STNAg7Ze++9mTVrVn2Fq9LV85Zx+vXz2TJgyxJpVgtOn7wXx07YqaVlMzMb6A499FBuueUWqOBRZqcNwarWmnQ/tEya3LE1ZdJ0jNmLV3L69fNZ38/XlvUBp153P7MXr2xNwczMulC3B+ml6X7nMmlyVcWlZdJ0jHNnLuo3QOcE8LEr7nKgNjNrkm4P0vek+z3Ssc7FTChI27EWrljDHUtWVZVnzbqX+eDFcznlqntYuKIrGhvMzFqmq4N0RDwG3A0MBo4pPC7pEJIx0U8Cc1pbutarp0b80/v+zNu/fSvHXjTHNWszswbp6iCdOivdny2pN/eipG2B76Y/Ts/NNtbJ1q57ue5z3LFkFVMvmcuP5z3WgBKZmXW3jurdLWlfNgRWgN3T/dckfS73YkQcmPfvayVdSLI4x3xJv2XDAhvDgRtIFtroeEO3aMzHIdepbMcRQ+jr7WnIOc3MulFHBWmSoDqxyOvjymWKiI9Lug04mWSykU2ABcD3gQu7oRYNNDSg5jqV/c9x+zlQm5nVqKOauyNiVkSov61E3isjoi8ihkfEVhGxX0Rc0C0BGmD8qGEcMHZkw86X61R23sxF/Sc2M7NX6aggbfWbNmkcgypaPKxy37ppIYefc6s7lJmZVclB2jbS19vDWZP3LD91Wg0WPLmG41yrNjOrioO0vcqUCWO44oSJDGtQR7KcwLVqM7NqOEhbUX29PfzPcfuhRlepSWrVHqZlZtY/B2krqa+3h+mT92z4M2rYMEzr6nnL+k9sZtalHKStrCkTxnD5RyYysYG9vnMCOPW6+Z6lzMysBAdp61dfbw9Xn3gQv/n0wbzrTds3/Px3LFnlTmVmZkU4SFvFxo8axvnv35cfulOZmVlLOEhb1Zrdqcy1ajOzhIO01aSZncpcqzYzSzhIW81yncp2225YU87vWrWZdTsHaatLX28PN37qYD592LiGz1IGrlWbWXdzkLaGmDZpPFec4Fq1mVkjOUhbw7SqVn3AV3/L5665j0tnL2HhijVNuJKZWTZ02nrSlgHTJo1n/51Hct7MRcxdsqrh5//Lmhe49q7Hufau5OcDxo5k2qRxXrfazDqOa9LWFLkJUM4+qvErahXyZChm1qkcpK2ppkwYw/SjmjNUK587mJlZJ3KQtqZr9lCtfO5gZmadxEHaWqLZncryuVZtZp3CQdpaqtlDtfK5Vm1mA52DtLWca9VmZpVxkLa2aUet+pSr7vHYajMbMDxO2toqV6s+d+ZCzrlpEdHEawXw0/v+zE/v+zM7bjOEt+22LVMP2pnxo5r/JcHMrBauSVsmtLJWDbD8mee5/PZHefu3b+XYi+a4KdzMMslB2jKjlc+q83kyFDPLKgdpy5xcrXri2JEtu2aug9n+X7mJM254wM+tzSwT/EzaMqmvt4e+3h4WrljD7MUreXD5s9y6+Cn+8uwLTb3uyrUvcvntj3L57Y96TnAzazsHacu08aOGbdSxqxUdzHJyzeCfPmw8p0wa14IrmpltzM3dNqC0uoOZx1mbWTs5SNuA044OZh5nbWbt4OZuG7By61af+fOHWPBk8wOnx1mbWau5Jm0DWruGbXmctZm1goO0dYR2DNvK8ThrM2sWN3dbx8gftnX+7xbx8/ueaEkvcNjQweyyOUt5xxu3dzO4mTWEa9LWccaPGsb579+3LTXr3DhrN4ObWSO4Jm0dK79mffmcpfzqgSdZufbFll0/1wz+rjftwCfe1uuatZlVzUHaOt74UcM48717cuZ792zpZCjgHuFmVh83d1tXaWcHM/cIN7NqOUhb1+nr7eHqEw/iN58+mKkHjqFn6OCWl8E9ws2sEg7S1rVyzeB3fuGwlo+zBq+8ZWb9c5A2o73N4O4RbmalOEibpfKbwd/1pu1bXrMGN4Ob2cbcu9usQG6c9fsmrOS8mYuYu2RVS6+fawa/4vZHOXj8a9ljh+H09fa4R7hZF3KQNiuh3eOs/7LmBa6963GuvSv5+YCxI5k2aRx9vT0tK4OZtZebu8360e4OZjluCjfrPg7SZlVoZwczcI9ws27j5m6zKhU2g9/8x6d4/OnnW1qGXI/wy29/1DOZmXUwB2mzGuWawYG2rLyVk5vJ7PLbH/Vza7MO4+ZuswZo58pb+XLPrU+56h43hZt1ANekzRqo3T3CwYt6mHUS16TNmiArPcK9qIfZwOYgbdZk7e4RnuMhXGYDj5u7zVogC83gsGEI12VzlvKON27vZnCzjHNN2qyFstIM7kU9zAYGB2mzNnEzuJn1x83dZm2U3ww+e/FKHlz+LLcufoq/PPtCS8vhRT3MsslB2iwDxo8atlFAPHfmQs65aVHLJ0YpXNTDQ7jM2svN3WYZlJWmcA/hMmsvB2mzjOrr7eHqEw/iN58+mKkHjmH0iCFtLY+fXZu1npu7zTKucI7wdi3qAR7CZdZqDtJmA0hWFvXIX4XLi3qYNY+bu80GKC/qYdb5XJM2G+CyMJuZF/Uwaw7XpM06RFZmM3OPcLPGcZA260BZGcLlHuFm9XFzt1mHykIzOLhHuFk9XJM263BZaQb3oh5m1XOQNusibgY3G1jc3G3WZYot6jHv0VU8+te/tbQcXtTDrH8O0mZdqnBRj3Y9uy5c1MOTo5ht4OZuMwOy8+zaTeFmGzhI10HSMZJukPS4pOck3SfpBEnt+vtm1hDtfnadawrf/ys3ccYND3gmM+tabu6uz2eApcBngaeAw4CLgJ2AL7WvWGb1KxzC1Y5FPTxHuHU7B+n6vCsi8seR/E7Sa4BPSfrPiFjfroKZNUpWFvXINYN/+rDxnDJpXIuvbtYebu6uQ0GAzrkHGA5s0eLimDVduxf1yDWDH/DV3/K5a+7j0tlL3BRuHa3jatKSdgUOByYA+wPjAQHHRMS1/eT9AHASsBewCbAAuBS4sIpa8VuApRHR2vEsZi3U7tnM3CPcukUn1qRPAs4BPgjsCpV1UpV0AfBDksD+e+AmkgD/HeBaSZtUcI43A1OAC2oqudkA4x7hZs3ViUH6AeDrJMGyF7ilvwySjgI+DjwJ7BURR0TEkcA44GHgSOAT/ZxjNHB1er1v1/MGzAYi9wg3a7yOa+6OiIvzf65wNNTp6f7UiPj7V/GIWCHpJGAWcJqk84s1e0vaBvgVsAo4MiJeqbH4ZgNau5vBwT3CrbN0Yk26KmkNeD/gReCawuMRcQuwHNgOOLBI/iHAz4GtgcMjYnVTC2w2ALgZ3KwxOq4mXYN90v2DEVFqEOg8YMc07R9yL0raFPgx8AbgLRGxvNyFJB0PHF9hufauMJ1Zpk2bNJ79dx7JeTMXMXfJqpZfP9cM/sv5T3DGEbu7Vm0DioM0jE33j5ZJs6wgbc53gSNIJjMZLim/pv1QRDxbkH4X4JAay2k2YBVb1OPWxU/xl2dfaFkZFjy5xuOsbcBxkIah6f65MmnWpvvC5Xnenu6/WSTPW0meZedbSgUd2VJ7kzShm3WMwkU9zp25kHNuWtSyiVG88pYNNA7SG4ZoVf13IiJ2qTL9DGBGJWklzcK1butw7WoK9zhrGyi6vuMYkBunMbRMmtwxj+kwa7C+3h6uPvEgfvPpg5l64Bh6hg5ueRncwcyyykE6aYIG2LlMmp0K0ppZg7W7R3iuKfzwc25l9uJiM/6atZ6DdDLXNsAe6XCqYiYUpDWzJmrnxCi5DmauVVsWdP0z6Yh4TNLdwL7AMcBl+cclHQKMJpmNbE7rS2jWndrZIzxXq75szlLe8cbtmXrQzu5cZm3R9UE6dRbJRCZnS/pDRCwGkLQtyTArgOleetKs9drZI9yzl1m7dVxzt6R9Jd2e20hqyABfK3j979LVsS4kmVVsvqSfSboeWATsDtxAstCGmbVZu5rCc53LTrnqHs8Lbi3TiTXp4cDEIq+Xnb0gIj4u6TbgZJKhT7mlKr9PdUtVmlmTtWuO8AB+et+f+el9f3bN2lqi42rSETErItTfViLvlRHRFxHDI2KriNgvIi5wgDbLpnb2CPewLWuFjgvSZtadcs3gu23Xug5eHrZlzeYgbWYdo6+3hxs/dXDLa9UetmXN4iBtZh2nHZ3LcrXq/b9yE2fc8IA7l1lDdGLHMTOzto2z9rAtayQHaTPraO0cZ53rXOblMa1Wbu42s67S6g5m7lxm9XCQNrOu044OZrnOZZ4MxarhIG1mXavVHcxyk6G8/du3cuxFc1yztn45SJtZVytcz3r0iFKL4TWWJ0OxSrjjmJkZG2YvA1i4Yg3n/24RP7/viaZ2MMs9r/7l/Cc444jd3QvcXsU1aTOzAuNHDeP89+/bsg5mngzFSnGQNjMroZUdzNwL3IpxkDYz60crh225Vm35qn4mLWkQsDfJcpDbAz3AEOCvwEqS5R1viwh/FTSzjpGrVbdiMhQ/q7acioJ0GpiPAD4MvA0YWpgENv7MSvoj8CNgRkQsq7+oZmbtN23SePbfeSTnzVzE3CWrmnqtBU+uYeolc5k+eS+OnbBTU69l2VQ2SEsaAnwK+CQwCv7+WOZF4I8kNeengXXAiHTbBdgO2A34EvBFSb8BvhQR8xr/FszMWit/XvDL5yzlVw88ycq1LzblWusDTr3ufoJgyoQxTbmGZVfJZ9KSTgIeAb5KEnRvBz4NHAAMi4g3RcSkiDg6Io6LiP8XEf8QETsAOwGTgUuBZ4DDgdslXSfp9U1+T2ZmLZEbtnXnFw5raueyAE69br4nQOlC5TqOXQBsAXwF2CUi+iLi3Ii4MyJeKnfSiFgeETdExAkkAf49wG3AkcAHG1R2M7PMaMXsZZ4ApfuUa+7+AnB+RNQ1yWxEvAz8DPiZpD5gm3rOZ2aWVfnN4M2aDMWdyrpLyZp0RHyt3gBd5JyzI+IXjTynmVnWtGIyFA/V6g5lx0lLGt6qgpiZdZpmT4biCVA6X3+TmTwp6UpJh0tq1YpuZmYdpdmTobhW3bn6C9JbAFOAXwCPSzpb0h7NL5aZWWdxrdpq0V+Q/iQwj2R89PbA54D7Jd0p6ROS3GPBzKwKrahVT71kLj+e91hTzm+tVTZIR8QFEXEgycQkZwGPkQTsfYFzgeWSfiLpSEle9tLMrALNrlXnJkBxjXrgq2iBjYhYGBH/ERG7kEwL+gNgLbAZ8G7gWuAJSedJ2r9ZhTUz6yTNHFsdwMeuuMuBeoCrehWsiJgVER8mmaTkOOAmYD3wGuBkYK6kByV9XtIODS2tmVmH6evt4eoTD+Lso/ZseK16zbqX+eDFcznlqntYuKKhI2qtRWpeqjIino+IKyPicJJpQE8FHiBpDn8DMB14VNKNDSmpmVkHmzJhDNOP2pNBTWj//ul9f+bt377V04oOQA1ZTzoinoyIr0fEm0ieV19I0tqyCXBYI65hZtbppkwYw+UfaV6nsjuWrHKnsgGmIUE6R9KBwL+SDNvyuGozsyq5U5nlqztIS9pF0hmSFgKzgROBkcDLwA0kq2GZmVkVmjlUy53KBo6agrSk4ZJOkHQryXKWXwZ6SWrPdwOnADtExOSI+L9GFdbMrJs0s1a9Zt3LHOem78yreGyzpEHAO4APAe8CNmdDk/YTwBXADyLioUYX0sysm02bNJ79dx7JmT9/iAVPNq6XdqRN3zuOGOLVtDKq35q0pH0lnQP8GfgpcDTJdKEvAFcD7wR2iohTHaDNzJojv1bdSG76zrb+VsF6gGRa0E8C25LUnOeQPHfeLiLeHxE3RsT6ppfUzMyYNmk8PzxhIsO2aNwkj7nx1F6gI3v6q0nvThKYlwFfAcZFxJsj4n8j4tmml87MzF6lr7eH/zluPxq9NqEX6Mie/oL0ZcDbImJsRHwxIh5pRaHMzKy8vt4epk9u/OQnXqAjW8q2l0TE8S0qh5mZVWnKhDGMHrEl581cxNwlqxp23txY6iCYMmFMw85r1avpoYYkAUeSzCa2EzAkIiblHd8K2A+IiPh9IwpqZmav1tfbQ19vDwtXrOH83y3iZ/c90ZDzBnDqdfO57u7lTJs0zr2/26TqcdKSxgH3A9eQdCB7J3BoQbJ1wMXALEn71llGMzPrx/hRwzj//fs2vFOZpxJtr6qCtKQRwG+BPUgC9RnAqzqQRcQrwHdJOp0dVX8xzcysEs3oVOapRNun2pr0Z0mat38FTIiIrwLPl0j7s3T/jzWWzczMapDrVNbIPmUeT90e1Qbp95D8rj4XES+XS5j2BH+BZLpQMzNroSkTxnBFE8ZTeyrR1qo2SI8Fno+IhytMvxZozpprZmZWVjOaviPgtOvd9N0q1Qbp3BrR/ZI0GNiaIs+szcysNZoxnnp94NnJWqTaIL0EGJz28O7PO0mGeFVa6zYzsyaYMmEMl3+ksctezl2yiqvnLWvY+ay4aoP0L0h6bH+2XCJJrwW+QVLz9lKVZmZt1oxlL0+7br6fTzdZtUH6m8DTwEclfUvSTvkHJW0r6WPAPcDrSFbOurAhJTUzs7pNmzSeK06YyMSxI+s+VzLhiZ9PN1NVQToiVpL08H4WmAYsJVkdC0krSdaVvgDYAVgFvDcinmtgec3MrE59vT1cfeJBnH1U/cO0PDSruaqecSwibgPeBFwFvETS/C1gZLp/hWSd6f0i4q7GFdXMzBppyoQxTG9AoPbQrOapaQBdRCwDjpN0ArA/sD1JwF8B3BkRaxtXRDMza5bcIh0fu+Iu1qwrO/1FWbmhWTuOGOJ5vhuo6pp0vohYFxG3RcQ1EXF1RMxygDYzG1gaNZ7aQ7Mar64gbWZmnaFRU4l6aFZjlQzSkl7X6ItJGiTJi5OamWVQ7hl1vTw0q3HK1aQXSPqBpF3rvYikzST9K7AIOL7e85mZWXNMmTCGA+ocnpUbmuUadf3KBek7gKnAg5JulnSipNdUemIl3irpIjaMl34tcF9dJTYzs6aaNmlc3dOIJoF6PsdeNMfDs+pQsnd3RLxZ0ruBrwGHAAcDF0haBNxFsp70SpLJTV4EtgFGkCzCsT+wD7AVybCsl0jGT58ZEU817d2YmVnd+np7OGvynpx23XyiznPdsWQVUy+Zy/TJe3HshJ36z2AbKTsEKyJ+KulnwOHACcARwK7p9v4yWXPfwf4EfB+4NCKeqL+4ZmbWCo0amgVJr28Pz6pNv727I/GriDiKZDz0McA5wBySGcfWktSUVwAPAtcDnwEOiIjeiPiaA7SZ2cDTyKUuPTyrNlVNZhIRq4Dr0s3MzDpcbmhWI5q+c8OzpkzwIJ9KeZy0mZmV1aihWeDhWdVykDYzs341YmgWJL2+T7veK2dVykHazMwq0oihWeDn09VwkDYzs4rkhmY1IlB7+tDKOEibmVnFpkwYw+UfmcjEBjR9+/l0/xykzcysKn29PVx94kGcXeda1H4+3T8HaTMzq0mu13c9gdrPp8tzkDYzs5o1YnjW3CWrWLhiTYNK1FkcpM3MrC6NGJ51+ZylDSlLp6kqSEs6QdLQZhXGzMwGpnqHZ11x+zJ3Iiui2pr094AnJF0q6eBmFMjMzAae3PCsWuO0O5EVV22QfpRk+cl/Bm6WtFDSaZJ2aHzRsk/SOEk3Slor6SlJ50vast3lMjNrh3qfT7sT2atVFaQjYizwj8CVwDqgF/gq8Kikn0uaLKmqRTsGKknbADcDw4Cjgc+SLN/5/XaWy8ysnep9Pu1OZBuruuNYRPwuIo4jWbbyJOBOYBPgncA1wJ8lfVPSGxta0uw5ERgBvCciboyIy4BTgCmS9mhv0czM2qfe59Mf/cGdbvZO1dy7OyKejYiLImIisAfwLeAvQA/wKeA+SXdIOlHS8MYUN1PeCcyMiPxP0nXAC8A72lMkM7P2q/f59KOr/sbUS+a6IxkNGoIVEQ9HxOeA0cB7gbmAgP2A77Khs9k+jbheKZJ2lTRN0hWSFkhaLykkHV1B3g9I+r2k1ekz5jslnSyp1D16A/BQ/gsR8QLwCLBb/e/GzGzgmjJhDMcduHPN+deHO5JBA8dJS9qMJECfCEzIvQw8Dwwh6Wx2ZxqsN2/UdQucBJwDfBDYNb1+vyRdAPwQ2B/4PXATMB74DnCtpE2KZBsBPFPk9aeB+ie1NTMb4KYeVHuQBnckgwYEaUl7SzoX+DPwY5Km3gBuAI4AhgNvJgmC64EPAf9V73VLeAD4OjCFpFPbLf1lkHQU8HHgSWCviDgiIo4ExgEPA0cCnyiRPYqdssTrZmZdZfyoYXVPctLtHclqCtKSRkr6pKS7gbtIgthrSJp6Twd2iojJEfHLiFgfEX+IiKkkQVvA+xpU/o1ExMUR8W8R8eOIeKTCbKen+1Mj4u9f2SJiBUnNHOC0Is3eT5PUpgttkx4zM+t6jViDuptnI6t2xrF3SrqGpNZ8DrA3SUepq4C3RsT4iDg7DXCvEhG/JulctmN9xW4MSaNJnpu/SNIzfSMRcQuwHNgOOLDg8MMkz6Xzz7c58HpgQTPKa2Y20DRiDepuno2s2jHNPydpyhVwP3AxcEVEFHs2W8o6KnxW3AK5jmwPRsTzJdLMI/lSsQ/wh7zXfwmcIek1EfHX9LUjgc3TY68i6Xjg+ArLtneF6czMMm3KhDGMHrEl/379fB5d9beq8+dmI9txxBD6ensaX8AMqzZIryWpNV8cEfNquWBE7FJLviYZm+4fLZNmWUHanIuATwL/J+lMYFuSYWhXR8RDFLcLcEhtRTUzG7j6env433/en7d/+9aa8uc6kTlIl7ddRFT/NSi7couFPFcmzdp0Pyz/xYh4RtLbgPOA60l6sf8I+Lcy51pKBZ3ZUnsDW1eY1sws83Idye5Ysqqm/LlOZONHDes/cYeoKkh3WICGDc3uNfXGjoiFwOFVpJ8BzKgkraRZuNZtZh1m2qRxTL1kLutrHAMze/HKrgrS3b6edK5ff7nlN3PHuncMgJlZg9Q7G9mDy59taHmyrtre3a9UuT0v6QlJv0tXy9q2WW+kRkvTfbkR9zsVpDUzszrUMxvZdXc/3lU9vautSavKbXNgFHAoyWpZCyRlaV7re9L9HpKGlEgzoSCtmZnVqdbZyLpt3elqg/RYkik3nwGeAL4EvI1krurdgLcCXyQZR/00ydKN+wAfBR4kmejjWkmvb0Th6xURjwF3A4OBYwqPSzqEZD7yJ4E5rS2dmVnnqmc2sm6aLrTaID2UZOjRw8AbIuLMiJgVEQvT7ZaI+ArJJB8LgO8B6yLiEpJJQ24BtgA+07i3ULez0v3ZknpzL6ZN899Nf5weEetbXjIzsw5Wz2xk3TJdaLVB+gvAVsAJEVHy6X1ErCGpPQ9L8xARLwGnkjSDT6qptP2QtK+k23MbsG966GsFr+eX9VrgQpJZxeZL+pmk64FFwO4kc5B/pxnlNTPrZvV2IuuGJu9qx0kfAjwbEQ/3lzAiHpK0mryAHBF3SFpH0oTcDMOBiUVeH1cuU0R8XNJtwMkk73ETkpaA7wMXuhZtZtYcUyaMYd7Sp7n2rserztsNPb2rrUmPALYosXTjRtI0W/DqRSjW0aRVotKmd/W3lch7ZUT0RcTwiNgqIvaLiAscoM3MmmuPHYbXlK8benpXG6QfJelkdWwFaY8l6d2dm1YTSVuSdB57qsrrmplZh6p1qs9u6OldbZD+Eckz5YskTSmVSNIxJB3MgmSu75zcM2KvEmVmZoB7epdTbZCeTjJkaShwpaSlkn4o6Rvp9kNJS0mC+VDg3jRPzr+k+9/UWW4zM+sg7uldXFVBOiLWkYyLvjx9aQzJWOhPp9v709cArgDelubJ+RJJJ66L6iizmZl1GPf0Lq7a3t2kQ6/+WdKXSdZP3gfoIWkGf4pkZq4bIuJPRfJ29hN+MzOrWT09vdeue7kJJWq/qoN0TkQsIVk/2czMrCH22GE4195Vfb6hW9QczjKt2gU27pZ0l6TXNatAZmbWvWrt6b3l4H5HBg9I1XYc2x0YV6wp28zMrF619vQ+/fr5HTlmutogvRxqfq5vZmbWr1p6eq+PzhwzXW2Q/jWwpaRiU2+amZnVLdfTu5ZA3WljpqsN0l8B/gr8j6TaHrwsP68AACAASURBVByYmZn1Y8qEMZw1ec+q83XamOlqu8P1Av8BfBP4o6TLSNZZfgp4pVSmiLi15hKamVlX+tuLJcNKWbMXr2T8qGENLk17VBukZ7FhcQwBp6RbOVHDdczMrMvVOva5k8ZMVxs8l9GkFazMzMzy1Tr2uZPGTFf1TiJilyaVw8zMbCO1jpm+/u7ljB81rOb8WVJtxzEzM7OWqHXM9Pzlq5l6ydyOGDftIG1mZplV6+pYnTJuuqYgrcRkSRdK+rmkmQXHt5J0sKS3NKaYZmbWjWodMw2dMW666qfrksYB15NMEZq7bYWdydYBFwOvlzQhIu6uq5RmZta1pkwYw+gRWzL9Vw8zf/mzVeXNjZseqEOyql1gYwTwW2AP4H7gDOBVdywiXgG+SxLEj6q/mGZm1s36enuYvO/omvIO5Cbvapu7PwvsBPwKmBARXwWeL5H2Z+n+H2ssm5mZ2d9147jpaoP0e0iatj8XEWXfdUQ8ArxAMkuZmZlZXbpx3HS1QXos8HxEPFxh+rXAwHwQYGZmmVLruOeBPF662iAdQEUra0saDGxNkWfWZmZm1apl3PSwLTblqTUvNKlEzVdtkF4CDE57ePfnnSS9xyutdZuZmZVV7bjpNeteHtATm1QbpH9B0mP7s+USSXot8A2Smvf/1VY0MzOzjdUybnogT2xSbZD+JvA08FFJ35K0U/5BSdtK+hhwD/A64M/AhQ0pqZmZGcm46cs/MpFhVXQIG6gTm1QVpCNiJUkP72eBacBSYFsASSuBJ4ALgB2AVcB7I+K5BpbXzMyM1w7bnDVVDq3KTWwykFQ9LWhE3Aa8CbgKeImk+VvAyHT/CnA1sF9E3NW4opqZmSVqbboeaE3eNQ0ei4hlwHGSTgD2B7YnCfgrgDsjYm3jimhmZraxbpnYpK4R3hGxDritQWUxMzOrSLdMbOKlKs3MbMDplolNav5KIWlTkik/RwCblUsbEbfWeh0zM7NCuYlN7liyquI8E8eOHHCrYdWyVOVY4Czg3cDmFWSJWq5jZmZWzrRJ45h6yVzWFy6WXISAUyZVMg9XtlQVPCX1AnPY0JM7gL+QrB9tZmbWMrmJTU6/fn6/gTqAc9Nx0gOpybvaZ9JnAq8BlgNHA5tHxPYRMbbc1vBSm5mZsWFik4kVzOl9x5JVA26K0GqD9NtIvpC8PyKu72+5SjMzs2br6+3hlEnjUAVThQ60KUKrDdLDSJaqnN2MwpiZmdXi3JmLiAqeTcPAmiK02iC9DBgkVfJ9xczMrPkWrlhTVS9vGDhThFYbpH9E0qN7UhPKYmZmVrVOniK02iA9HbgPuCgdimVmZtZWnTxFaLXjl48FLgX+E5gv6VpgHlC2zSAiLquteGZmZuV18hSh1ZZwBknv7twz6anp1h8HaTMza4pOniK02iB9K0mQNjMzy4ROniK0qiAdEYc2qRxmZmY1q2aK0EEaOFOEehUsMzMb8HJThA7qZ4DwIMH0yXsNiKZu8MIXZmbWIaZMGMPoEVty3sxFzC3S9D16myG8dbdt2XvMNm0oXW3KBmlJpwDPRcQlRY4NBQZFxLNl8n8bGB4RH6m7pGZmZv3o6+2hr7eHhSvWcPmcpdy84Ckef+Z5AB5/5nkuv/1RLr/9UQ4YO5Jpk8ZlvkbdX3P3OcB/lTi2COjvKf37gOOrLJOZmVld7ln2ND+cu+zvAbrQQFlso5Jn0uVa+D09qJmZZcrsxSsrWr5yICy24Y5jZmbWUc6duaiiXt6Q/cU2HKTNzKxjdNpiGw7SZmbWMTptsQ0HaTMz6xidttiGg7SZmXWMTltsw0HazMw6RqcttlHJV4eRkn5X7HWAEsc2SmNmZtYKnbbYRiVBejBwaJnj5Y6BV80yM7MW6qTFNvoL0j9oSSnMzMwaJLfYRn8TmgyExTbKBumI+HCrCmJmZtYo/S22MXHsSE4ZAHN3Z7M7m5mZWZ3yF9uYvXgla9e9zNAtNqWvtyezz6ALOUibmVlHGz9q2IAJyoU8BMvMzCyjHKTNzMwyys3dZmbWVQbSM2oHaTMz6wqzF6/k3JmLik50csDYkUzLYG9vN3ebmVnHu3reMqZeMrfkTGR3LFnF1Evm8uN5j7W4ZOU5SJuZWUebvXhlvxObAKwPOO36+zO1bKWDtJmZdbRzZy6qaIpQSAL1eTMXNbdAVXCQNjOzjrVwxZqqFtsAmLtkFQtXrGlSiarjIF0HScdIukHS45Kek3SfpBMkqd1lMzMzam66zkqTt3t31+czwFLgs8BTwGHARcBOwJfaVywzMwNYu+7lluZrNAfp+rwrIvK/bv1O0muAT0n6z4hY366CmZkZDN2itjBXa75Gc3N3HQoCdM49wHBgixYXx8zMCtQ67jkr46UzHaQl7SppmqQrJC2QtF5SSDq6grwfkPR7SaslrZV0p6STJTX7Pb8FWBoRf2vydczMrB/jRw3jgLEjq8ozcezIzMxAlukgDZwEnAN8ENgVqKhDlqQLgB8C+wO/B24CxgPfAa6VtEkzCivpzcAU4IJmnN/MzKo3bdI4BlXYnXeQ4JRJ45pboCpkPUg/AHydJPD1Arf0l0HSUcDHgSeBvSLiiIg4EhgHPAwcCXyiSL6tJe1WwbZlieuOBq5Oy/jt2t6umZk1Wl9vD2dN3rPfQD1IMH3yXplp6oaMdxyLiIvzf65wZNPp6f7UiPj7iPSIWCHpJGAWcJqk8ws6dh0JXFrB+Q8DfltQrm2AXwGrgCMj4pVKCmpmZq0xZcIYRo/YkvNmLmJukXHTE8eO5JQMzt2d6SBdrbQ2ux/wInBN4fGIuEXScmBH4EDgD3nHZgAzarjmEODnwNbAQRGxupaym5lZc/X19tDX2+NVsNpon3T/YEQ8XyLNPJIgvQ95QboWkjYFfgy8AXhLRCyv53xmZtZ840cNy2xQLtRpQXpsun+0TJplBWnr8V3gCJLJTIZLOjDv2EMR8Wx+YknHA8dXeO69G1A+MzMbwDotSA9N98+VSbM23Tfia9Tb0/03ixx7K8nz73y7AIc04LpmZtYFOi1I53qWVbjeSX0iYpcqsyylgh7qqb1JnnObmVmX6rQgnVu2ZGiZNLljLV/ipJrOaZJm4Vq3mVlXy/o46WotTfc7l0mzU0FaMzOzTOq0IH1Put8jHRpVzISCtGZmZpnUUUE6Ih4D7gYGA8cUHpd0CDCaZDayOa0tnZmZWXU6Kkinzkr3Z0vqzb0oaVuSIVMA072MpJmZZV2mO45J2pcNgRVg93T/NUmfy70YEQfm/ftaSReSLM4xX9JvgZeASSRLSN5AstCGmZl1uazPPpbpIE0SVCcWeb3sEiUR8XFJtwEnk/SQ3gRYAHwfuNC1aDOz7jZ78UrOnbmIO4rM433A2JFMy8g83pkO0hExiwqXpyyS90rgyoYWyMzMBryr5y3j9Ovns77EjBp3LFnF1EvmMn3yXhw7YafiiVqkE59Jm5mZFTV78cqyATpnfcBp19/P7MUrW1OwEhykzcysa5w7c1G/ATpnfcB5Mxf1n7CJHKTNzKwrLFyxpugz6HLmLlnFwhUtn6Dy7xykzcysK9TadN3OJm8HaTMz6wpr173c0nyN4CBtZmZdYegWtQ1oqjVfIzhIm5lZV6h13HM7x0s7SJuZWVcYP2oYB4wdWVWeiWNHtnUGMgdpMzPrGtMmjWNQhVNkDRKcMqnsBJdN5yBtZmZdo6+3h7Mm79lvoB4kmD55r7ZPDZrpaUHNzMwabcqEMYwesSXnzVzE3CLjpieOHckpnrvbzMysPfp6e+jr7fEqWGZmZlk1ftSwTAXlQn4mbWZmllEO0mZmZhnlIG1mZpZRDtJmZmYZ5SBtZmaWUQ7SZmZmGeUgbWZmllEO0mZmZhnlIG1mZpZRDtJmZmYZ5SBtZmaWUQ7SZmZmGeUgbWZmllEO0mZmZhnlIG1mZpZRDtJmZmYZ5SBtZmaWUQ7SZmZmGeUgbWZmllEO0mZmZhnlIG1mZpZRDtJmZmYZ5SBtZmaWUQ7SZmZmGeUgbWZmllGbtrsAZmZmWbBwxRpmL17J2nUvM3SLTenr7WH8qGFtLZODtJmZdbXZi1dy7sxF3LFk1auOHTB2JNMmjaOvt6cNJXNzt5mZdbGr5y1j6iVziwZogDuWrGLqJXP58bzHWlyyhIO0mZl1pdmLV3L69fNZH+XTrQ847fr7mb14ZWsKlsdB2szMutK5Mxf1G6Bz1gecN3NRcwtUhIO0mZl1nYUr1pRs4i5l7pJVLFyxpkklKs5B2szMuk6tTdetbvJ2kDYzs66zdt3LLc1XKwdpMzPrOkO3qG0Ecq35auUgbWZmXafWcc+tHi/tIG1mZl1n/KhhHDB2ZFV5Jo4d2fIZyBykzcysK02bNI5BqiztIMEpk8Y1t0DFrtvyK5qZmWVAX28PZ03es99APUgwffJebZka1HN3m5lZ15oyYQyjR2zJeTMXMbfIuOmJY0dyShvn7naQNjOzrtbX20Nfb49XwTIzM8uq8aOGtT0oF/IzaTMzs4xykDYzM8soB2kzM7OMcpA2MzPLKAdpMzOzjHKQNjMzyygHaTMzs4xykDYzM8soB2kzM7OMcpA2MzPLKAdpMzOzjHKQNjMzyygHaTMzs4xykDYzM8soB2kzM7OMcpA2MzPLKAdpMzOzjHKQbgBJm0i6W1JIOrrd5TEzs87gIN0YJwE7tLsQZmbWWRyk6yRpFHAmcFq7y2JmZp3FQbp+Xwd+DcxqcznMzKzDZDZIS9pV0jRJV0haIGl9pc98JX1A0u8lrZa0VtKdkk6W1ND3K+lgYDLw+Uae18zMDGDTdhegjJOAadVmknQB8HFgHTATeAmYBHwHmCTpmIh4pd7CSdoUuAA4KyIek7RLvec0MzPLl+Ug/QBJU/KdwF3AJcAh5TJIOookQD8JHBwRi9LXRwE3A0cCnwDOLci3NbB9BWVaFhF/S/89DRgCfKPC92NmZlaVzAbpiLg4/2dJlWQ7Pd2fmgvQ6blWSDqJ5LnxaZLOj4j1efmOBC6t4PyHAb+V1AN8meQLwRBJQ4DhaZotJW0dEasrKbCZmVkpmQ3S1ZI0GtgPeBG4pvB4RNwiaTmwI3Ag8Ie8YzOAGVVcbjQwFLisyLEfAKuBbao4n5mZ2at0TJAG9kn3D0bE8yXSzCMJ0vuQF6RrsBh4a8Fr2wFXkdSwZxbLJOl44PgKr3EQwL333suhhx5aQxHNzCyL7r333tw/e/tL20lBemy6f7RMmmUFaWsSEWspGHKV13HswYi4rUTWXejnuXqh1atXc8stt1RXQDMzGwiG9pegk4J07s0+VybN2nQ/rMllKWUpUGnEfTOwCUnz/ZxmFaiL7Q1sTfJo4t5+0lptfI+bz/e4uZp1f3tJYtaS/hJ2UpDO9SyLdlw8IpbmlaFUmhlU+Oxb0iySWveciDi0rsLZq+Td33t9f5vD97j5fI+bKwv3N7OTmdRgTbov13yQO7amTBozM7NM6KQgvTTd71wmzU4Fac3MzDKrk4L0Pel+j3TccjETCtKamZllVscE6Yh4DLgbGAwcU3hc0iEk45ufxB2xzMxsAOiYIJ06K92fLenv488kbQt8N/1xesFsY2ZmZpmU2d7dkvZlQ2AF2D3df03S53IvRsSBef++VtKFJItzzJf0WzYssDEcuIFkoQ0zM7PMy2yQJgmqE4u8Pq5cpoj4uKTbgJNJus5vAiwAvg9c6Fq0mZkNFJkN0hExi37GHZfJeyVwZUMLZGZm1mKd9kzazMysY2S2Jm3MIJkffGlbS9G5ZuD722wz8D1uthn4HjfTDNp8fxXRllk0zczMrB9u7jYzM8soB2kzM7OMcpBuIEkfkPR7SaslrZV0p6STJdV0nyUdLuk3klZJ+pukByT9h6TN+8k3UdJPJP1F0jpJiyT9t6Sta3tn2dCI+ytpkKR/kPSV9FyPS3pR0gpJv5T03jJ5vywpymzrGvNO26dRn+F671Wj/y9lSYM+x7v0c3/zt4ML8nbk51jSrpKmSbpC0gJJ69P3c3Sd563p91Xr3+9C7jjWIJIuAD4OrANmsmESle8AkyQdExGvVHG+fwPOBl4h6bjwNMm4768AR0iaFBF/K5Lv/cDlJOPDZwPLgQOBzwNHSuqLiL/U+j7bpYH393Uk9wVgFXAncFv6+juAd0iaAfxLlO6wcR/F15Z9qbJ3k02N/gynqr5XTSpHJjTwva0FflDm+O4kaxWsAe4qkabTPscnAdMaecJaf1+1/v0uKiK81bkBR5GsY/0EMC7v9VHAQ+mxaVWcb39gPfAcMDHv9aHALen5vl0k32jgb+kH4z15r28K/CjN95N236923l/g9el/tsOBTQqOHULyxy+ADxfJ++X02JfbfU+yfI/ruVeNLkeWtla+N+CX6fm+16jfTdY34ATgv4Fj0//ns9L3eXQrf1+1/v0uWY5239hO2EhqYwF8qMixQ/J+0YMqPN+1aZ4vFjn2ujQIvwBsU3DsG2m+7xfJNxxYnR7fvd33rJ33t59rfSE938wixzryj1sz7nEdQbplv+uBfo/LXGdH4OX0fBOLHO/Yz3HB+6w3SNf0+6r173epbcA/32k3SaOB/YAXgWsKj0fELSRNztuRNDv3d77BJM2uAD8scr4/kaziNRh4Z8Hh3PPUYvmeBX5WkC7zGn1/K5BbxnR0A841ILThHme6HM3Q4vd2PMnjrgcjYm6d5+pKtf6+6vz7XZSDdP32SfcPRsTzJdLMK0hbzq7AlsCqiHik0vNJGk7SxJN/vJ5yZEWj729/cnPDP1Emzb6Szpb0PUnTJR2Z/uccqJp5j6u5V63+XbdSK9/b8en+kn7SddrnuJFq/X3V9Pe7HHccq9/YdP9omTTLCtJWcr5lZdIUO98u6f6ZtNZcbzmyotH3tyRJWwKnpD9eVybpu9It3+OSjku/YQ80zbzH1dyrlv2u26Al703SIUAvSQ3win6Sd9rnuJFq/X3V+ve7JNek6zc03T9XJs3adD+siedrdDmyopXv67sk/3EeAr5X5PgjwOnA3sDWwGuBt5F0BhkN/FLSm+osQzs04x7Xcq869TMMrXtv/5LufxoRT5VI06mf40bKzN9h16Trl1upq1Hzq9Z6vkaXIyta8r4knQH8M0nnumMj4oXCNBFxeZGsNwM3S7qWpDfoV4EjmlnWJmj4Pa7xXnXqZxha8N7SR165McHfL5Wugz/HjZSZv8OuSddvTbofWiZN7tiaMmnqPV+jy5EVTX9fkj4D/BfJN9x3RMSDNZzmv9L9YZI2q6UcbdTqz06pe9Wpn2FozXt7H8nz0MeBX9d4joH8OW6kzPwddpCu39J0v3OZNDsVpK3kfGOqPF/u39uk36jrLUdWLE33jbq/G5H0SeCbwPPAERExp9pzpBak+8FAT43naJel6b4p97iIUveq1eVopaXpvpnvLdfUPSMi1td4joH8OW6kpem+2t9X7t/V/v0uyUG6frkhO3tIGlIizYSCtOUsIAkYIyW9vkSaAwrPl3YWy/UmnPCqHCXyDQCNvr9/J+lk4DyS2YTeXWdnmdfk/XttyVTZ1LR7XEKpe9XqcrRSU9+bpN2BiSTNrJdWX7y/G8if40aq9fdV09/vchyk6xQRjwF3k3zzPKbweNrbcjTwJMn4uP7O9yLwq/THDxY53+uAg0h6b/6i4PD/lck3nA09OX/SXzmyotH3Ny/fx0im9nsBeG9E/LbOoh6b7v8YEQOqKbZZ97iMoveqDeVomRa8t4+k+5vTsbi1GrCf40aq9fdV59/vkoXxVv/MNkezYfaZ3rzXtwUepMj0ccAnSL51XVbkfBPYMK3cAXmvD2XDLDrFpgXdiQ3Tgr477/VNgasYuNOCNvr+fjS9v+uAd1ZYhjHAB4DNC14XMDW97wGc2O771e57XM+9qqUcA2Vr9Oc4L81mwIo0/we6+XNc8J5yfytLzjgGnJXe37Ma9Vmkxr/fJcvY7hvZKRvJ8J0gaer4GXA9G6bh/Amvnif6y+mxWSXO92/p8ZeB3wA/zvuPeDuwZYl870/zrAduJZmze2mabxGwbbvvVTvvL8mwk/XpsYeBGSW2bxTJF8CzJJMRXJeW40/p6wGc3+77lKF7XPO9qrYcA2lr9N+JNM2RaZqngS36uX7Hfo6BfdO/jbnt2fT9LMx/vSDPjDTNjEZ+Fqnx73fRc7X7xnbSRvINdXb64XiOZPWZkykyF2+F//kOB25K//M9T/Lt7T8o+BZcJN9E4AbgKZLm3MUkE89v3e571O77Cxya98eo3La0IN9r0nt4M/AYSY1jHckXoB8Bb2v3/cnQPa77XlVTjoG2NeHvxM/SNBdUcO2O/RxX+n+7IM8MygTpej6L1Pj3u3BTejIzMzPLGHccMzMzyygHaTMzs4xykDYzM8soB2kzM7OMcpA2MzPLKAdpMzOzjHKQNjMzyygHaTMzs4xykDYzM8soB2kzsyZT4jOSHpG0StKPJI1qd7ks+zwtqJlZk0n6DPBNkhXqXgC2JFkKcWJEvNzOslm2uSZtZtZ804A7gO2AYcBXSVZtenM7C2XZ5yBtZtZEkgSMBq6NiJURsZ5kCUSA7dtXMhsIHKTNBiBJsySFpOPbXZZuJek2SS9L6i2XLpJnivOBD0naVdLWwBdI1jW/vch5B0laIGmtn1ubg7RZC0makQbX/O0lSX+VtFjSDZL+XdLYFpbpeElflrR3q6450El6N9AH/CgiFleQ5VPA64EFwDPAR4FpEbGkMGFa0z4L2Ao4o2GFtgHJQdqsPV4CVqTbKpKORK8H3kPyvPIRSddIem2J/MuAPwKrG1CW44EvAQ7SFZA0CPgaECS/q0o8SPI7z3koIr5TJv0PgSXAv7byC5tlj4O0WXv8ISK2S7dRETEEGAG8A7iaJAAcDdwjacfCzBHxoYjYLSJ+0tpiG/BPwB7AbRHxcIV5/hMYDvxv+vNukjYrlTjt8f0DYDPgE3WU1QY4B2mzjIiIZyLixoh4H/D/gHXAjsB17S2ZFTgh3f+oksSS3kDSvP0gcGr68mBg136yXpXup5YL6NbZHKTNMigibgQ+l/44UdK78o+X6zgmabCkaZL+IOmZ9Jn3Ckn3SbpA0kFpuuMlBXBImvXSgmflSwvOO1LSP0u6Lu3YtEbSc5IekvQtSTsUey+SlqbnOzQ9x7ckLZH0gqTlkv5XUtlezpLeIOl/JC1Mr/mMpPmSzpO0X5l8b5T0/fR669J8syV9rJbAJ+k1wLtIWjquqTDbN4BNgdMi4mngsfT1vcplioiFwH3Aa4Ejqi2rdYiI8ObNW4s2YAbJH/hZFaQdDDyZpr+q4Nis9PXjC17fNO9YkPQgfhp4Oe+1H6Vpp6TnfzF9fXX6c26bV3Dub+SdI5c+/7x/AfYq8j6WpsePy/v3cyQtBbm8S4ARJe7DJwuusxb4W97PRe8lSTPxKwX58s9zM7Bllb+/o9O8f6ww/T+m6W/Je+3G9LXpFeS/IE17Ybs/u97as7kmbZZREfEi8Lv0x7dUmO0DJDXjvwFTSYLQCGBzYGeSwHVfev6rI2I74A9p3mmx4Tn5dhExoeDcy4HpJJNwDIuIrdPz7g/8mqTGd2U6LriY80m+MPxDRGwFDCXpKPcMsAtwemEGSccA5wGbANcCu0fEUJKezzuQBP67iuR7T3q954F/B0al+YYAbyfpdHco8O0SZS2lL92/6ppFyjCIZJYxgH/LO/Rgut+zguvdme4r/f1bp2n3twRv3rppo4qadJr+dDbU/DbLe30WxWvS36XKmlepc1X5vjYnCT4BHFJwbGn6+pPAa4rk/Wx6/E8Fr29G0jQcwJVVlGWTvGseWSLNWJKa9UvA9lWce3Z63tMqSHtCmvbagtf/JX39sQrOsS8bWkSGtfvz6631m2vSZtn2dN6/R1aQ/tl039KZrCLiBeCm9Me+Esm+FxF/LfL6Del+rKSt8l6fRDJT1yvA56sozqEkrQZLo0Tv90jGJ99O8njg0CrOnbuvK8slkjQUOJOkef3fCw7natKjJY3o53q56wjwxCZdyEHabOCoZDWcX6X790j6qaTJaWenhpC0m6TvSLpf0rOS1uc6mpHMTw1JM3Qx80q8vjzv39vk/fvAdH9fROSn6c8/5Moh6clSGxu+TOxUxbl70v3TZVMlvbi3Ay6OpANYvofy/l2281jBdXpKprKOtWm7C2BmZeXXtPoLDETELZK+CHyRpBfyuwAkLQB+AVwUEYtqKYik9wGXkTRDQ9IEu5pkVSdInjFvlW7FrClR5nV5j7Hze1znao7LqixqrrY7mMpqn1tWce7N0/2LpRJIGk3ShP8cyfjojUTEGknLgDEkz6VvKXO9dXn/HlJFOa1DuCZtlm25zkWPR8RLZVOmIuJMYDzJ8+xfkzSB70YSOB6S9KFqC5HOfPa/JEH0apLOYltExIhIO5qxoRNWqY5jVV+2xny5v2s/iQhVsH25inOvSvfblElzFklA3Qp4osg0sEESoKH/mnT+l7RijwqswzlIm2WUpMEkz2UBfl9N3ohYEhHTI+JwkmfZbwVuJWk9+66kbasszjtIasoPAR+IiLuKfGlo9DPTJ9P9zlXmW5Hud29gWXJyz4iLPktOx2x/sIrzVROkyz4Ht87kIG2WXR8FcsH0h7WeJCJeiYhZJBNivERSw9s/L8n6dF+u5jo63d8fyQIQG0mHXb2t1jKWkFshaq9iU6OWMSfd7yppjwaX6Y/pvtR82t8iuY9nkATYUttJafo3lhmyBsnQNNgwht26jIO0WQZJ+ifg6+mPcyLiFxXmG1zm8IskPaVhw7NV2NAjvFwTbm4hj1JB5aMkC4Q00kySTmWbsOFeVJov9xz725I2KZWwgt7VhWan+/0LD0iaDBxMUuM9J5JpXotuwP1ptq2A15W5Xm6s+uxiX46s8zlIm2WEpK0l/ZOkq4BfkjzXfIxklqtKXSbp0vQ8gZp7fQAAAn9JREFUw/LOvQvJgg1bkEzwkd98nhsSNFnJWsfF/Jakd/kbgfMkbZOed7ikz5PMjNXQZ6Zpc/pn0x/fL+nHknbLHZe0vaSPSjqvSL5PpuU9DPiNpIm5LxeSNpW0n6TpwJ+qLNZt6X6f/OCfTjF6dvrjNyJibT/nye/xXa7JOxekq3rcYR2k3QO1vXnrpo0Nk5m8yMZTcD7HxlNurifpoNVT4jyzKD6ZyQ0F53i64NwvA1ML8uxG0kM7SJrDl5NMBnJbQbpvFZRxFRum2bwR+Er67xkF+Zamrx9a5r7kzrlLkWOfYePpPddQ2bSgH857X0Hy5WQlG08NGlX+/gQ8kuadlPf6p9PXngKGVniuv6Z5vlTi+BYkrRzrgde3+7PrrT2ba9Jm7bEZSUerUSTjX18gqdX9FPgPkj/KUyKi2s5Cp5FMQXljer7BJM3FjwCXAvtGxOX5GSJiAUmN80aSZu3tSDprjS5I9xngX4F70vJuCtwLfIpk1a6XqyxrRSLiW8A+afmXkty7dSRNxueSBMhi+S4lWWnqHJLWgpeBrUmC480kC5jsUmVZAvh++uP7IFl4hOQZNFRWi87J1aZL1aSPAIaRfAl5pJpyWudQ8pkzM7NKpKt9LSWp0e8QyWxrzbjOdcBkkt70V/WX3jqTa9JmZlWIiD8DF5EMbftwM64hqZdk8ZGHSB57WJdyTdrMrErpOPNHSJ5xj4uIhjb1S7qEZCGOIyPihv7SW+fytKBmZlWKiL+kM7e9ieTZ/dJGnTtd4vIR4PMO0OaatJmZWUb5mbSZmVlGOUibmZlllIO0mZlZRjlIm5mZZZSDtJmZWUY5SJuZmWWUg7SZmVlG/X8J1mrTMv3kbQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 486x486 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(6.75, 6.75))\n",
    "ax.scatter(rs, energies, c='C0')\n",
    "ax.set_xlim(0, None)\n",
    "ax.set_xlabel(r'Distance ($\\AA$)')\n",
    "ax.set_ylabel('Energy (eV)')\n",
    "ax.set_yscale('log')\n",
    "ax.set_ylim(1e-4, 1e4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.9"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}