notebooks/water_MD.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Molecular dynamics simulations of bulk water\n",
    "\n",
    "In this example, we show how to perform molecular dynamics of bulk water using the popular interatomic TIP3P potential\n",
    "([W. L. Jorgensen et. al.](https://doi.org/10.1063/1.445869)) and [LAMMPS](http://lammps.sandia.gov/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "%matplotlib inline \n",
    "import matplotlib.pylab as plt\n",
    "from pyiron.project import Project\n",
    "import ase.units as units"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "pr = Project(\"tip3p_water\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating the initial structure\n",
    "\n",
    "We will setup a cubic simulation box consisting of 27 water molecules density density is 1 g/cm$^3$. The target density is achieved by determining the required size of the simulation cell and repeating it in all three spatial dimensions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2635f5900e0d4ab29f46adec6f56e2a7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>NGLWidget</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "NGLWidget()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "density = 1.0e-24  # g/A^3\n",
    "n_mols = 27\n",
    "mol_mass_water = 18.015 # g/mol\n",
    "\n",
    "# Determining the supercell size size\n",
    "mass = mol_mass_water * n_mols / units.mol  # g\n",
    "vol_h2o = mass / density # in A^3\n",
    "a = vol_h2o ** (1./3.) # A\n",
    "\n",
    "# Constructing the unitcell\n",
    "n = int(round(n_mols ** (1. / 3.)))\n",
    "\n",
    "dx = 0.7\n",
    "r_O = [0, 0, 0]\n",
    "r_H1 = [dx, dx, 0]\n",
    "r_H2 = [-dx, dx, 0]\n",
    "unit_cell = (a / n) * np.eye(3)\n",
    "water = pr.create_atoms(elements=['H', 'H', 'O'], \n",
    "                        positions=[r_H1, r_H2, r_O], \n",
    "                        cell=unit_cell)\n",
    "water.set_repeat([n, n, n])\n",
    "water.plot3d()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'H54O27'"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "water.get_chemical_formula()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Equilibrate water structure\n",
    "\n",
    "The initial water structure is obviously a poor starting point and requires equilibration (Due to the highly artificial structure a MD simulation with a standard time step of 1fs shows poor convergence). Molecular dynamics using a time step that is two orders of magnitude smaller allows us to generate an equilibrated water structure. We use the NVT ensemble for this calculation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['H2O_tip3p', 'H2O_tip3p_slab']"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "job_name = \"water_slow\"\n",
    "ham = pr.create_job(\"Lammps\", job_name)\n",
    "ham.structure = water\n",
    "# Listing available potentials for this structure\n",
    "ham.list_potentials()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use the `H2O_tip3p` potential. The `H2O_tip3p_slab` should be used if you want to switch off the periodic boundary conditions along the z-axis  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "ham.potential = 'H2O_tip3p'\n",
    "ham.calc_md(temperature=300, \n",
    "            n_ionic_steps=1e4, \n",
    "            time_step=0.01)\n",
    "ham.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "383b3d53a6324a9eace3d6bb81c16f58",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>NGLWidget</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "NGLWidget(count=101)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "view = ham.animate_structure()\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Full equilibration\n",
    "\n",
    "At the end of this simulation, we have obtained a structure that approximately resembles water. Now we increase the time step to get a reasonably equilibrated structure "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the final structure from the previous simulation\n",
    "struct = ham.get_structure(iteration_step=-1)\n",
    "job_name = \"water_fast\"\n",
    "ham_eq = pr.create_job(\"Lammps\", job_name)\n",
    "ham_eq.structure = struct\n",
    "ham_eq.potential = 'H2O_tip3p'\n",
    "ham_eq.calc_md(temperature=300, \n",
    "               n_ionic_steps=1e4,\n",
    "               n_print=10, \n",
    "               time_step=1)\n",
    "ham_eq.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e9996981cf4047d4bca6bd6a62db4fb0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>NGLWidget</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "NGLWidget(count=1001)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "view = ham_eq.animate_structure()\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now plot the trajectories"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEKCAYAAADenhiQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzsfXfcFMX9//uzd0+hdxQFBLEgiqKiYosoGmtCLFFT/GkSTTU9JlhijEmMJUZTjearxhI10dhij70gIirSREVFBGmidJ5yd/P7Y3d2Z2dndmfv9soD8369ntdzN7c7M7s7+/nMpxNjDBYWFhYWFlnCqfcELCwsLCw2P1jmYmFhYWGROSxzsbCwsLDIHJa5WFhYWFhkDstcLCwsLCwyh2UuFhYWFhaZwzIXCwsLC4vMYZmLhYWFhUXmqAtzIaKLiGgJEc30/o7x2puJ6EYimk1ErxPRROGcZiK6jojeIqL5RHSi195CRP8iogVE9BIRjajHNVlYWFhYBMjXceyrGGO/k9rOAgDG2FgiGgzgYSLahzFWAnA+gBWMsZ2IyAHQ3zvnawA+YYztQESnArgMwClJgw8cOJCNGDEiq2uxsLCw2CLwyiuvfMQYG5R0XD2ZiwpjADwBAIyxFUS0GsB4ANMBfBXAaO+3EoCPvHMmA7jI+3wXgD8TEbGEvDYjRozAjBkzMr8ACwsLi80ZRPS+yXH1tLmcTUSziOgGIurntb0OYDIR5YloJIC9AQwjor7e778ioleJ6E4i2spr2xbABwDAGCsAWANgQA2vw8LCwsJCQtWYCxE9TkRzFH+TAVwDYBSAcQCWArjSO+0GAIsBzABwNYCpAApwJayhAF5gjO0F4EUAXKVGiuGVUgsRfZ2IZhDRjJUrV2ZzoRYWFhYWEVRNLcYYO9zkOCL6O4AHvHMKAH4o/DYVwNsAVgHYCOAe76c74dpaAJcZDQOwmIjyAPoA+Fgzp+sAXAcA48ePt+mgLSwsLKqEenmLDRG+Hg9gjtfenYh6eJ+PAFBgjM3z7Cf/BTDRO2cSgHne5/sBnO59PgnAk0n2FgsLCwuL6qJeBv3LiWgcXPXVQgDf8NoHA3iUiEoAlgA4TTjnZwBuIaKrAawE8BWv/XqvfQFcieXU6k/fwsLCwiIOdWEujLHTNO0LAeys+e19AJ9StLcB+HyW87OwsLCwqAw2Qt/CwsLCInNY5mLRZbGxo4B7Xltc72lYWFgo0GhBlBYWxvjFfXNx5yuLMaxfd4wf0T/5BAsLi5rBSi4WXRbL17UDANa3F+o8EwsLCxmWuVh0eVi/cwuLxoNlLhZdFqrUDBYWFo0By1wsuj6s6GJh0XCwzMWiy4Ks6GJh0bCwzMWiy4NZ0cXCouFgmYtFlwUXXGwmOQuLxoNlLhZdFmT1YhYWDQvLXCy6PKzkYmHReLDMxaLLwsotFhaNC8tcLLo8rOBiYdF4sMzFosvCmlwsLBoXlrlYdHnYwqMWFo0Hy1wsujBc0cWyFguLxoNlLhZdFlwtZgUXC4vGg2UuFl0W1uRiYdG4sMzFYjOAFV0sLBoNlrlYdFlYbzELi8aFZS4WXR6NbnN5bdEneHze8npPw8KipsjXewIWFuWCuojV5fi/TgUALLz02DrPxMKidrCSi0WXR4MLLhYWWyQsc7HosrCuyBYWjQvLXCy6LHzmYmUXC4uGg2UuFl0WXcXmYmGxJcIyF4suD6sWs6g2FqxYj988OM/msUsBy1wsui6s4GJRI5xx43T8/bn3sGT1pnpPpcvAMheLVNj1wkfw16cX1HsaIdi9pEW1USq5q8yW1jaHZS4WqbCho4jLH3mz3tOwsKgp7AYmPSxzseiysHtIi1rDrjlz1IW5ENFFRLSEiGZ6f8d47c1EdCMRzSai14lootfeSzh2JhF9RERXe7+1ENG/iGgBEb1ERCPqcU0WtQdXUVgjq0W1YZdYetQz/ctVjLHfSW1nAQBjbCwRDQbwMBHtwxhbB2AcP4iIXgFwt/f1awA+YYztQESnArgMwCnVn75FvWF3kRa1Ao+lsiYXczSaWmwMgCcAgDG2AsBqAOPFA4hoRwCDATznNU0GcJP3+S4Ak8ha3bYo2F2lRa1gY6vMUU/mcjYRzSKiG4ion9f2OoDJRJQnopEA9gYwTDrvCwD+xQJdyLYAPgAAxlgBwBoAA6o/fYt6w24hLGoFu4FJj6oxFyJ6nIjmKP4mA7gGwCi4qq6lAK70TrsBwGIAMwBcDWAqgILU9akAbheHUgyvXApE9HUimkFEM1auXFn2tVk0Fmz6lyjmL1uLN5aurfc0NjvYDY05qmZzYYwdbnIcEf0dwAPeOQUAPxR+mwrgbeH7HgDyjLFXhC4Ww5VuFhNRHkAfAB9r5nQdgOsAYPz48ZYidXHY91yPo652tcY2zX82sMQiPerlLTZE+Ho8gDlee3ci6uF9PgJAgTE2Tzj2CwhLLQBwP4DTvc8nAXiS1cB96JE5y3D+PbO3KE+lRr3WBp2WxWYEvsbshsYc9fIWu5yIxsHdECwE8A2vfTCAR4moBGAJgNOk804GcIzUdj2AW4hoAVyJ5dRqTZpjzcZOfPNWV3j64RE7YWDPlmoP2RDIiogzxnDu3bNx0t5DMX5E/7L7sX4bFrWD3cGkRV2YC2NMZhq8fSGAnWPO217R1gbg85lNzgAbOwMzUJPTaA531UNWr1dHsYQ7Xv4Ad7+6BG/95uiy++GsxUou6TBnyRowBowd2qfeU+lysEvNHFsOZcwQOXHHvAVtnhtOLebXc6k+3l25Ho/PW16DkaqP4/70PD7z5+frPY0uiUZ7BRoZlrmUA5GhbEGLrVEvtRZM77Arn8GZN8/ItM+XF36MEVMexJwlazLtNy1OvGYqzr7t1brOodHBl5j1TDSHZS7lgIkfo4utvVDEwo821HBCtUF2Npds+uEBbV31df+fJwm9sOCjus7jlfc/wQOzltZ1Dl0FjSS5MMYwYsqD+O1Db9R7KkpY5lIGisIKUy22Kf+ZjYm/expr2zprOKvqI6tdm3/PslIpNtALnwa1kLimvvNR46kzuyAa8Q7yx3rts+/WdyIaWOZSBkpM/Bxdds+97QZotnUWazWlmiAzySWjV5WbvlTPoNp4e/m6LkG0v/j3l3D79A+q0ve/Xl6Esb941K91sjmDP+tGutJ6rPs0sMylDIgvk+rx8mfuWFdZJbJ+J2r9ij01fwWOuOpZ3DtzSSb9VXuZLP5ko/a3SjZAF90/D+vaC2grbF6bqDg00oaicWaihmUuZaCUoBbjv29urCU7ySUb8Ptb6x3cghXrAQBzlnSN9Cr5nP41X9cmZ1cyR1POfQKdhUYnc5WDX2ED8RYruWyOKIqSi+IB+yaFzUxyyUqdlRXz5be31u9YVuq4Ws27ydHf6U0d5UsdzXmXfHQUS2X3YVE+Gpy3WOZSDkQVs+r5bq466LjF3FksYUO72S44c2+xGr9ljl+kLJv+qp3GPZeLYS4VqMWacl2LuSxb04YFK9bFHvPbh97ALdPej7Q3IiFvxDmJ0EboE9GPDM7fwBi7NsP5dAmwBLUYb2p0sTUtVFfz1vJ1+MtTC7BibTtefHeVWaLELm5zyXmSQNrnyxirizQbl0UiE+ZS6BrMZcJvnwAQn8yTe16dNmG7ULtv0G+gV7rRY27i0r+cAzc1ftzb8E0AWxxz4Tp3QENgeMBVYz/71FBJCD+4YybmpUzt7qvFKqSzvnqqxpIi1zIVU47LWH1StucFyUW+V1wtVs68eL9/fOJtEAG/P3lcwhldF77NpYEIeqMrSOKYyy2MsYvjTuYZjLc0fOufQTSz0luM/9/MuIvqanIx+vw0/aiwtq0TvVubMusvKzi+5JLuvBJjcIS9Wq3mnReekazC4t5irflc6n65RHTPa67X3ObMXDga6ZVudPoSZ3O5MuY3AABj7KcZzqVLQN75qR4w35k39qPPBmXwFiMVw5vL1mH3ix7DPa8tNugv/RzKBWNMsLmklFw07dWWZkRvsU6JuWz0JJfWpvTm16b85uWwYoJGeqcbXXKJW1GvE9H/iOirRGTTp3qQd35KrZjXttnZXBSXU44NweSuLFntxmbc89qH2mPqEUR538wPfYaadtx6rQdRupQJEre5tDaVIbnEuDhvdvBV3Q30TjfQVFSIWx3bAvgdgIMBvEVE9xLRKUTUrTZTa0yYMJdSAxr/MoHiemS1mIn9Q2dzeXPZOoyY8iBmfrAa/Xu4NXJWrmuP6ak8w3ol+MG/ZmLmB6sBAGmdpORpXv/8exnNKh5Ngs1FJI7n3Pk62r0ASO5WnK7fLYe5NOKr3OibV+3qYIwVGWOPMsa+AreM8I0APgfgPSL6Z60m2GiQPWNUBr7N11ssej05iUMUTGR1zSFPv7kCAPDQ7KW+ZeKj9Xrmwm0JtfaEbet0B0ytFqvTchAzRYiP585XFldUYbF5S2IuDajqbqS5qGC0OhhjHQDmAXgDwFoAY6o5qUaGrLNW0dJGdFvMAmq1WPh7oZRM6XW3Jeczi4CNxUlCwfHRMc+7ZzZO+OsLiXMpB5xYF1PbXOqzIMRNjuzhxnwpMj17KcdO09XRSO90o29eY1cHEQ0nonOI6FUADwDIAZjMGNuzJrNrQEQkF1WE/ubqiqxok9ViJpJLEKEfPtcn2iVm9OLw41Vj3vbSIry6aHViH+WAb9jTe4tVPvZtLy3Crhc+ksr9WuS98nqtZErdm+tVJd1dIx9v6KjZeEzxqd5odPqiZS5ENBXAcwC2AvB1xtjOjLFfMMYas3hAjRBVi0XRiD7xWUDFSOXknIVi8jXrXgrRUJ7mxSmW3LoWF9w72/ykCpDzXHDLCaIsFxs7Cli1vh0/v28ONnQUU0lN4jxlnsSZVDlqsZ6t9WMuv3nwDez1q/9hXY3LWjQSQW8o5wIF4iSXcwGMYIz9hDGWbQm+Lox2I8nFbWt0V8G0UF2OE5FcslGLmewQOfPmksut0xYlnpMFuK2nEldkUb366wffwB3T4+d+3B+fx96/ftz/niaAU5ymzBD9b2VwFzk25s1l62IzMGeJB2a5XoQb2muTkbkR6XgDTimEOIP+M4wxRkQ7EdETRDQHAIhodyK6oHZTbCzINhe1txj/LfvHf+u09zFrcXXUPUlQXY6ctkqUXG57aRFGTHkQqzd24KV3V/meX7r74ghpVYz8Arxj0kbKV4owEzQHE5bOT++aFfrtd4+9FXvuu1Jl0zRSU7zNxbibRBx59bM46LKnsuswBvwyyomzKgd8I9NIBL1L21w8/B2uFNMJAIyxWQBOreakGhkmajGOatC8C+6dg8/+uTqGahnn3j0Lfxeq3KnUfLJaTCRePAHg4k824ZTrpuHU6150++EeShJhyIsGfc6gDeZpoorLEjlNhH5bZ9F37VVBvH//fT0cv2Mi8QEBY07D2EKJVmW1WAUZquup9q3EEaGycWs6XCwaaS4qmDCX7oyx6VJb+UUgujjkOJf43UODP/0E3D79A/xGrM+t9BYLv9yyZAcA671sye+sdHffqlv2zsr1uOp/bwNwXYs58fh4Q4dfa15GQGhr64vMd8uyBLb7RY/hwEv1O3eRyEe87AwZZOBFZ3S4Ny5Tfq4U9SRu3OZUTvqhchBsdhrnnd4cJJePiGgUvHVNRCcBWFrVWTUwZM+kuOcbt7l8e/k6vxxytXH/6x9WVLeDQ+0tFv6u2lG/s3K9d6xnq1D0dOp107BsbRuAqFrsu7e/GjlehFFsTYbgDFUetqNYio3LiVOTmksu8MZOY3PRMxeTbi64dzZGTHnQeLxagDsi1Nqo3Uj0vJHmooIJc/kO3MzHo4loCYAfwM2GvEUi4sopPeC3l6/T/ibiiKuexWnXywJh9nh54cf43u2v4eIH5lXcl9Lm4siSi+C06Z1w/j1zAACtXhS42E+hWMLGjkIoEt+NcwkOGjmwp3o+wvG1xF2vLI6Ma5J2PiS5IL2XnYg03mLikTJzCbIl6CUAnaNEPXfOAZOt0XjSuI2ARpqLCom+hIyxdwEc7mVAdhhj8dV2NnPIWh/5BTviqme1v6mwqaOIbs3p8zqZgrtqLl2zSXvMJxs6sK6tgOEDusf2pZI4ZKIUR+hbvPxVop7/azfNwDNvhSW4ImMhisiZ1E1TF2Ls0D7Ya3g/r939vdaSC4+vEJ+vSV2UOJVK2mtIF+ciSi7ynNLj1UWf4MrH3sQLC1aVcXY2CJLD1pi7NBC6rFqMiI4TvzPGNsiMRT5mS0A6dUTyMR/GEP0sEFRr1B9zyBVP4VNXJHv5qL3FJMklRr1T8Diz2I3MWACXGIpEkBPeX9w/Fyf8dWowH6Q3bmeJ0P1I4d0GVJ4JOV2ci/hZI7mkGPuEv06tK2MBhOuv8aNvJJtL48xEjTjJ5QpPDRa37i6BG7m/xSCacj/4LBO5OEbUknfQXiihvbPKxmhufI45ZG2bYXliRZtsT12+pk17PleZJdFFWS1W0CQPq5fkwlEMSQTJcxAPqXTG6bzFhHlKt7Iam99iiWHpmk0Y2i9eEq4EpRrzFt8VucwBz7p5Bv43b7lZpVZDNLrkEsdclgP4fcL5b2c4ly6BSISz8IBNjbKA4PWTZgdaBhHltD8Lw6cyQl/iLm8tX4+jx6rP5552SXORI/Q7NfaIwOYSvu/T3q3NrlpkgCZ3V3zWlT6PNN5icUGU1cDvHnsT1zz9DqZOOQzb9NUnUV/X1ol7XluC0yZsl9qlOAhUTr6etW2dIQm7HuWmdR6PlaDBeYueuTDGJtZwHl0G2ghnpJNcshjbBFm+RCZqsY2deimI35+kqygxM6btSy4S83n5vY8TRgAWrdoIBobtBpRfTDUkiZhILsLnSoWt8tO/SGuUp3/JkNY+//ZHAIAn56/ArMWr8bOjRmNAz5bIcRfcOwf3zfwQo7fujX1H9k81hr+WDG7D7hc9hhahpECJRYN/k+Br4RqIoJusubbOIi64dw5+euTOGNy7tQazClC/5EBdFFFXzuC7vMOOe/Z8badhGGmz8NYCchBlWtuD9hjhs05y4UfJajEVIZPBbUyVqCnSMotSSjVaHMoNoszCoJ8EviQuuNf1EnzpvY/xzDmHRo5b8smm0PFpIKvFZn6wGs05B2O26a08XkzbVGIMORAWf7IRA3u2GBVKa7w3z2xOD89ZirteWYxCsYSrT61tvuEtL2d2heBEYdywvt734Lc0kkugFjMfuxx6lKXwrxrfkVaQke0h4bVgjElMO14HJP/ev0dT4hwyQciGku7hVLpPSKNWi5Nc4rqZsfBjrNmYPjGkvObeX7URu174SCR7wQYv9qpbGVUwOTjD/txfXsAxf3zO7BxvfR102VM4+7b4GCoZjWTQT2vnqzXqwlyI6CIiWkJEM72/Y7z2ZiK6kYhmE9HrRDRROOcLXvssInqEiAZ67f29csxve//7VXPunI6dtPdQr0WvvjF5rmmIRNxutbNYwl+eWoA2jUvshvZCIpEuB7LkYraDd//rVHZL17SFHB10MSD81snXlZc5XpUQsrlorvs9ISdYlmrSNFLsP6YuxLteIKtst9OVPyiWGE7624v4fze85LcZr1XFc93QUcSq9eEU+XytRqTfKoOxYJ0+/sYKw3MqM+hXA0YaAO+YWtuYAAPmQkQziOg7VSDaVzHGxnl/D3ltZwEAY2wsgCMAXElEDhHlAfwBwKGMsd0BzAJwtnfOFABPMMZ2BPCE971q4C9jkBk3+E0mgiYPP43kEkdQ/jntfVzx6JuhXGBA8J6/umg1vnZTZcmtTXKLhe0Q5v2ImL9sHb7/r5n+9ySbS0QdGTqmetTA5Fr/PeMD/3OWTm1p1GIr17XjpL+9qJyDVuHoXdCsJWv8NtMhTcnYxo6C12/5N6acU93cdeWN2UC8xei+6UqK1wImW7xTAWwD4GUiuoOIjqTqscExcBkEGGMrAKwGMB7ueiUAPbyxewPgmf8mA7jJ+3wT3FLMVQNflI4ieaGuyp8K5XhxsRjBg+fvkoP5xB3ps4qYkjQwidDPSlQXI951hFQX5xJXv8QUpRJLlPTCNhf1QGKd+ZcMvNiKJYZz757tSxr6+SV2FQIPptXZDGWGH2T2FtvMbqYpddjoqcUqYS7lnGuadTt8jvu/kWqopLFdypJpLZDIXBhjCxhj5wPYCcBtAG4AsIiIfklE6Vw8wjjbU3HdIEhFrwOYTER5IhoJYG8AwxhjnQC+BWA2XKYyBsD13jlbMcaWenNdCmBwBXNKBKc3qpoekbxjMf2INpf17QUsWZ0cTBknuQQpyKu3iFSj1/Nl40PHMfVyAyx/+p9Z2PH8hzE9xvNMHEc3SrPgljTl7tmJ1RPnLFmD26cvwg8FyU2FtM4dQT40tXQtd6eSLo2Zi+Gc+LOpJAi2nDNLLJ3tRLQ7NQ5rMZTafLVYVaeihJFymoh2B3AlgCsA/AfASQDWAngy5pzHiWiO4m8ygGsAjAIwDm4SzCu9024AsBjADABXA5gKoEBETXCZy55wpahZcMsApAIRfd1T881YubK8XTx/wYIkjAHkYD+TuBTGGE7861QceKn2VkbGjvtNDmosd1HFFUELj6ueh/xZdUyl6533HlejpNxdMc8fdvK1LyaOzxjDGx+uVR7TJGX21NnEOO55bQkAoDkf/2qmJci+d6IcRKmpU6K6beYmF/WTjTgWev2leUYfb+jA/GXBvS5nc8OkOKrEvupAmE1gwiD5MfW4hERXZCJ6Ba566noAUxhjPMPgS0R0oO48xtjhJhMgor/Di/JnjBUA/FD4bSrcQM1x3u/veO3/RmBbWU5EQxhjS4loCACthY4xdh2A6wBg/PjxZVEdmbmE4zHMJZfAFRl4c7k6XRvzxHd/rBiC4ksuEncpd1EpXz7luHrCvlZTgjYtPdDbBdRziEt3kiV41/95dQl+cufrymPyctroBPxj6kIAycwl7XVxwq67V0kJWdOMqVtzutPT+Jkc98fn8KGQBaJcyUWFNAbyRkCagnqNKrl8njE2iTF2m8BYAACMsRPKGdRjAhzHA+BVLrt7CTJBREcAKDDG5gFYAmAMEQ3yzjkCAC80cj+A073PpwO4r5w5mYIT+Lxf1CP4rRxX5Lid128efAOjznvIHzNuMfFjImqxMheVSu2iJjjyd+bPZ/ladfr5rN5Pnc1F3NGpdvhZqfJ4L3MEo7eMZilaTzey/Nha8vHuuekllwS1mHS8Wi1mOJZmzemeUxpG+aGUXqicZylngAjmo0Z4DPPxVm/sKNtDc317wbej6mBy7bVQl+tgEkT5OYWYuwbAK4yxeMWwHpcT0Ti4T2ohgG947YMBPEpEJbgM5TQAYIx9SES/BPAsEXUCeB/AGd45lwL4NxF9DcAiAJ8vc05G4A8r57m7io83spBMdhbC50KxFNrp3vyiW8mxo1hCq5NLsLmEJao0c1BBTbyS2/hpcfm+fGJCQK+WPNYlvERA/Isk35e4oEEgLigzJViyzUBWi2lLPEvvWHOCxJM2FRDvXmtbiWwSon2YMjSd8TgyNpc8K7G5lHGqmxhVxTzVnZl4Bcpo6yxi3MX/w2kTtsOvPrebcH6QemZjRwFvLluHPYdHHXF3+8WjaG1yMP9XR+uvw2AuPDVSPSQXE+Yy3vv7r/f9WAAvA/gmEd3JGLs87aCMsdM07QsB7Kz57W8A/qZoXwVgUto5lAvZFTmuPrlRzIdwfnshzFyacoSOotve2pSLfQk5gZUXUbogTaY1/Lq/K+Yf2Zgl70bFnw7YYQAenRufd4mx8DlPzV+Br/zjZWEO8m6caX8DotVEywXvOS6nXJS5qI+T9wQtTUlqsaTZhcG719mnIpKLoc0tdjAJus1RJZknylWLKbdJOnVZGeO9/sFqAMAL73wUGYO/oz/61+t4ZO4yvHLB4cqsEm2JSW2TZ9NR4xLgIkzUYgMA7MUY+zFj7MdwGc0gAJ9CID1sMfAN56o4lxj1jL6/4HO7VHCKMxrulmtCsKNBjeaLK2nHr+opmqvKZNyAopkSSfGwM28Ox+vIxEll0L9l2vsYMeVBfLKhI7MU/XycuEJfeVktpmUu4eNkpiQjO28xFvrvtyv6qDTORbcmXlu02qzjFH1yqJ1Q1HEuSQ4oacDf5V6t4WwRYl8zPQYkv/cmhefcvpKP4X01ZBAlgOEARP/JTgDbMcY2AdDXdN1MIbsihwz6EpFJL7mEPYmaPMLEd9pxRDGwuej7TzMX1VhKaSbynUsu+nHKUTPEZRSWBYfQdXifb3/Jraa4ZPWmDG0ubj9pUv7rmIJswE96buWqxaI2Mvd/5Dkq6JvpWtLp999buSEU68N7+/3/3jLqV4WkKameTYmlS+JS0XqVThCfv06VLdMBk3npoCtXUQuYqMVuAzCNiLih/DMAbvcM75XXzu1iiHVFltO/GLmfBB/l2i5NEcklbl7u/7iI+SSE635ET/z2P6N5mHTeR7GMkP9E5qqW0Asu/Rbniix7yzCWXXEx3nec0VYeSze2bGNpSkhhU64rss4rTGuLEVBpECWXOP1koSkfw5WPvRlpS2QuCqlSVrNy6G0u+o2NDkz6L44tj1euKtssQt/7n2V6CEOYBFH+Cm5altVwDfnfZIxd7FWm/FK1J9ho4ItL5S2WJimgmBWZ9xVVi/H2cCSzLJ2Iv+kqDZpA5I3iDuvaZ95BqcTw7soNkXN0Rtq4lzBkEzGcX1hyCf8WNehHJTCfuSB9dPYHH29Uts/9cC3aOotmjFSajwy+kejZ4u73Rg6KLwVQvlos3B7EuSRL3aZDLlgRn12gXPzpyQWRtqT1o7KH6dVi6j6Y5nMcdOv/5hcX4rG5y+LHK+OdSDqmHtVaYyUXInIAzGKM7QbgldpMqbHBN6mqOBddgJoK/JcScwlLoVRUqMXUNpc4/WlS2v81GzuxYOU67L1dNLmCLgDytw/Px7D+6qqC0Z2wgVosNKb+OFNEDPrCZ5+5COWe06rFZi3Wuxo/MGtprPeZTAC0zCXvzq97cw7r2wuJhDztTpRvSOLsU+H26A+mBGrFOjNteS0yDKskFze3mGI+BjYXY7VYKTj+jumL/PZLHpoPALj72wcIbvvyeGZjJN0+xpgfa1aPaq2xkgufq6geAAAgAElEQVRjrATgdSIaXqP5NDwiajHhmUV20AbqTsaY/+LL6hWuGvnsn19w+xfSZXxHUlFxfpMUa3P6jdNx4jUvKnWxcYxSnpuOiZjslMQIffOdoP63iCOFcPBDs5cqxjcc1EOPFn3MCSFaCVM3F3ds9/v+2w+Q+glvVkTCu04RjCqutQvunY0JlzyhnQMQbEh0m4GoWiyKrANSs+guS5tLJcGVIn7075m46nHXjsTAMOXu2ZFjbn3xfX+D8MKCwKNsXVtnYhaHD1dvwvXPv5e4jv/vufdw4wsLAdRHcjEx6A8BMJeIniCi+/lftSfWqHDddQPbRmgXniJC3z9HOEhexHw3q/r9QQXRBKIvk7ym5n64RnkcoFeLqcB/1lXmjJUOQtdcuQpAvO93zvggdM2/fXg+Vq1vF9Ri6Ylk3NFESTE94e/8JS8xhv2E6osyUxan+MbSaAYHkVjcOm0Rlq1tixyjno8kSWmYi6krer1hUhtIxjF/fB7fuEWRIVzTVUhyMRjv7leXYK6XDkh3z/I58tfGj4XsDmMveszfTALAa4s+wd2vLg6d+9V/vIxfPTAPS9cE+QjfWLoWI6Y8iPdXBaprsbRyozKXXwI4DsDFcHOA8b8tEkXGQkbzuAy8ZhmCWUhFJkKuSxK3QPhQ8i5afrn4Dla3o/M/S79HHAWE/6ME+4Bv0Ndc+/urNhhlE5YRd5Q41jl3zYq80IUSCxm00zKX9pidpEMUa9DXEfMSYyEvIX5YwWcuwXnqLAPRsRat2ogPNQlQ+VCykKVdU4pmXlkyK2RB7pJopmoddhRKeHnhJ4q+1J2FmhPGW5WQmJQjn3O04320PlArHv/XqfjRv1+Xfo+Oces0N+D6GSHzuRgr1ZDMhTH2DNwo+ibv88sA0pVv24xQYm7deE5rv3FLYIoy1WfL/QXHh09okmIk4uNc3N+SJBdOZBLVYhGmpD6WScw2yebylRtfrtgVWYZR2h1vjjzrcBrEBbMRmbmIcwSSC9T3rRS9f6rrUY35qSuewgHaBKjuWPIaLSqYGaCmoZxwdW8uv3Jk1kiSfFMFEUvfn3t7Jb76j5clySUe0YBe9XFNDpVN8FVqWJ5tu3+PZr9NLN/ccDYXACCiswDcBeBar2lbAPdWc1KNjFIprBaTfxNhovIRxWz5+cuBdPFZkd3/0VgbtQSiMkKL48sLX+fiXCq5BPb608dL7eq5buwohhwTjCWXmMN0qjkOQuCd9/P75uIvT71jNCYHr5Gzw+Ceyt/jDfrh76LdTLylssTHYp6FeJwpAlds97x+3Zu8vr32yLz1/Q/XOHekRRbxRmmJfeyx0nzOunkGnpy/wq87AyRvhqLeeGrESS5J4IzinLtm+W1cYurtBW3OWrw6xFyqmcBVBxO12HcAHAg3xT4YY2+jyjVTGhlcnSEShhFTHvR/E2HyOEtM774rZ9SNi4fSGdJ1zKVQKuH5tz8KZS6Ok1yiKcuY/98hwqRdtsKOg3smJiNkktbaOEI/heQSOZbU7tum4MTlwuPGRH4jit+BRjMYBJJCziFMP28S9t9+QMRzqMQYOgolPPvWSiUjSR1E6c/H/T953Lah+UUM+jHdD1SkKikHWZC7JAaVhqbKzyrH3xVh85BoczG8qqacU7anpMoDbn1bkJ/v8XnL8dk/v4D/vv5hcE4jSi4A2hljvpLPKzlc+5k2CLg6Q5WcL01WZA7X5sIJdRhNAkXUJdsT56UaM0JnvS6XrWnDl69/Cd+7/TXlfOW1KLs/Bwb94DcikTiq58lYQBDWbOqMLcYlwiSA1P+u2D5Wkv5ik1eOt5tCHeRQ1JNulaAzj0gufBPgqRMH927FyEE9InnKGIBfPTAP/++G6ZijcIVOK7nwTUUQK+WpyTizi8S5qPtfsc7McSAtypVi3lyWULHToN8RUx5014x0KL9HHUWziHl3vPB33XU15cyldrkf1WaGr5sSY3hbEWfUkEGUAJ4hovMAdPPS4N+JIInlFgeuzojs5FnUd95k7YjHyItNVIsVGYtdIKJL6W0vLcLuFz0aYlwc/IVZ6+10xIC3NGoxcVwSjtGlwRfHKGeZp9FPy0emHU8ei9tcWhVp8FesbY/Edez968eDsaVnyneQpVJwT8k7rlAs+c+AMYbXF7u5p9oU6UDKTv/inccDdO/3drema/esm2ZkFp8ijlHuzvq8e6JuviJMCXihFA2udRTBzUndmTLJvOOkkqrEZKudCpsLl2ZKjCltMnHJVasFE+YyBcBKuCWGvwHgIQAXVHNSjQyuzpBp7caOYmykuA4lgSlFbC75sLdH3Psn7kDPv3c21rYV3HOkNeXH1PgJ7YS5CANEmUu4n0CVB3CnNteGwtuZ8jxAHR2dhDTMReW1p5Nb7hdUBxyyJMIJX6ScAYCLH5gXW7pYpxYrsSC+yfHum1g0jjFg7SZXZcmZmujgUXYlSu+0vHQtpr0tW9tWFZfkOI+7a55+B9OEnGRpYKx2jShsg+ctJpJM6s6UScsJTZMg2vVUfXKGUypp7Kl1SDFm4i1WYoz9nTH2ecbYSd7nLVYtFrgihxfHOo+Yi0iT+wdQeIsJBKBQYgn1XHgfAuFXzCEQ9T3mIlxHOd5iLuEO78DF+cju1CWmjtsQcc6R0aoLf3vG3AgfDVzU57vi5YxFyGl4iqUSiNLVxPjS/03Dk/OXx8a5cOLleDnWxDQzDAxrPOaiupdpM6nLWZFlRikvLd3a7SwyTH2nPEIfhw3tetXTZY/Mx6nXTSur3zQOIzrPyhBzycjGo1pLcX0nZUouFgN1a5zarJYwKXN8IICLAGznHe/RELZ9dafWmDhp72GYsP0ADJKMmuvaOiOLw6QoVTgpXvg30U5QLMbbXHy1mNAmSkVyn3ynqPJYkj/LcxHHKbJAfeA4YryNQMQEusEYwxWPRhMQihjSpzXSxsv/miBCKEtMaSMDgGeFuAAO+UXuLLobijTM5YUFq/Dye5/g+4fvGJ6LwHx9tZgnuYiqoRILVJdczZHPkZuTHOXr0PlpUSks3J9uqcVJaZVgzaZODOqVjaOAiDTMRRcTJqqkxCOefnMFLv7vPLQXSnhhymHe79J91Mg6bR1RZhonjSZVtOwUHEVUarN6lHUxyYp8Pdy69q8gRCa2TIwb1hfjhvUFABw2ejCenL8CALChI5rAMK7OBwcTbBBRY7ygpkqwuai8ftzdmCzqu//l3bk8fpLaRYzL4BI+gUISjTteevWLSv2UBpEXnCG23LNYJA2IvsjFkqvC0jEoHTqKJbyzMmxcDWwuzGfKRO79KkiqD/4MOOEQMyendS2N2FzKlFyyRN4h/36s2VQdpmW6YWeIbsS4t5hOajjjxpcjbaZqsQ0K5nKnQormSJJcCn5ZjkCKERGXoqhaMLG5rGGMPcwYW8EYW8X/qj6zLoADRgX5oTZ2FCK7AxNRtCRwl7gXvFAqRaQJwM2Y3NZZVP6mFvXDLwyFjterxXRBdsUS89U1Drnj/eOF93Dcn54HEGUUJjvuSgsbpbG5ANHKlPKLXCi5zKecad396pLQd1EtJtpcmCS5qLyDRD19kTFs7EguD80he4vlJHWl/FRqsdF1HPJzrK3eGM2fBkCbccAUpkxy1uI1/kaRI6cw6PMbI6+R+2YuwcaOgnE4gkoSOVeRg4zjk43xzFc06KucI+LUjtWCCXN5ioiuIKL9iWgv/lf1mXUBiB5UmzqKEQJsUlFOLJQUTZcffC6W1LrUo69+DqN//ohwbphB6IMooxXqiiU3azIQZQKRjM+C+sunU0R48d1VuOi/QZmfcgzHuYqZi3wf40f9zYNvhL7Ltq1CseRJLpXDj0dizL9O0VuMQ5wBJxyizeX26Ysw5sJHsfCjaBkEFeRiYfJzkVELq2qpxNCvhxv0t1aRnBNATMYBM5hex6nXTcMv7p8bauO3O2zQdzuUXbK/f8dM/OK+uYr6LTrbVZQ2xEnsYr4xZX/eC/rd218L5RTjeO+jDXh5oZnbf1YwYS77wS1tfAmCvGK/q+akugrExbBRpRYz2KWvWNcuBB6GfwtJLhqby7secRG9t8Tzo3pk9z/fjb0nEKdZi1djj4sfw72vLUn2fOPMRTJMywzVibiZAbtt2ztyHSISKvwmQjXVOH5184vvAwBWe7tDVcqWtDYXHThNcTMbcFuVwuYifOYSsFit8oOP3R39/GXxzhEcgbeYRl2ZIKlWA4USQ4vnCVetgomVqPdkKR8APlztMpV1bVGpcfm6dqM0Om6f0V8q2VSJKYqWaKS9a595Byf/7cVE+01WMPEWO1Txd1gtJtfoEAnnJoUrcvoSo3qj6twP1+CNpWu1Z6psLqqYEv7CbFLofHkm1+cXfBSRVGRDrpgRQIzXkCHvkEuMYcSA+EJYupgaU6gIZZK9ZM6SNRh38f9wz2uLI0y+UGJePrny5iXeAtEVOee7cPP4BEEtJpzPHUPkXHOAuiTuBfdG1Su+t5jGrbrWajHZ9nPD8+9hhWFmZxlxjDAT5iK8xxfcOwf3zVyi7Nf1+pMnp+5bReATio/ivHtmK70bTfH4GyswfeHHWLq6OoGwMkxyi21FRNcT0cPe9zFE9LXqT63xIe40NnZECzyZeIuJiJNcvnnrq8pKfByqSGumcAJY7r3A3M1VhBiIKTPKx98Ii9phmwuXXKLEL2JzYcmqikoN+vfODMeuxLkic8zzGPfzb69SliamCtRi4n0pCMwlYMoEhjDBCTtXeN5iCuqjcsy4dVo0Macc56JyRf5ofTuueHQ+Oosl/NPLslsNMBa41fMUR/OWrsW+lzyB1xZFsxUn9xdte3TuMpx18wz/Pl5/+nj88rO7puqX3yI5K/bzb3+kHDNHZMyU5Zott0x7PzZBKuAGR//kztdjjzGBypusGjDxFvsHgBsBnO99fwvAv+B6kW3RENU3GxXlbkVi8dbydViwYj2OGTsk0g9fqGltBeFjw33xNplhcWIkpvXmuOnFgKB0SkQrmv4lkFxyMcxFpdtPuq5KJRe5zG6JsUTmwgMWHYrOr1AqwXEqkVwInB2fd89szFu6BsWS4MLtxbkUhd28OAXOkJSSS0JhKbEPngOPjyGCMYbz75mNR+cuxx3TPzBOHV8ORE+4Zuma5i9bhz2H90vVn2o98WzlZx3sRky05HM4YsxWEbuKDl+/eQbe8cp6t0tSxtq2TiVzUSVi1a30xyS7yM8zLmcQBxMv1ixgot0eyBj7N4ASADDGCrAuyQCkHWmRhYzzQFiH/umrnsW3/6muVOCb4mVbQYo1oItz0RHyVYqaEBwbOwo48+ZwMSVdVIS4A1dt7WWbi8rlU8T3J+0YtdNUiFIp2Y34155Rn1TMhce5lDm+LHDcOm0RxAqkPEI/YCJOSNXjG/QVxiiV5KLCIiFAU5XVmyFI0FlNxgKEVYDR5KzpCZ94SnuhGFL5nnztiwBcBp5GIhaJ/7XPvBv6be2mqFcYH6MeLt1p0TA2FwAbiGgAPHpCRBMA6IuKb0EQF6vKQyPtQ4zzFjM9V+xiySebtIT8xZh0GiqVmS7lfrHEfBdZ1bsbtbnoX7hvTxyFHx6xU0jd+OUJlVfYNpFcOBxFluMCj3Mpk7uoJLGQtxhxbzFOcMPqFT6fbk3R3Ga/ljzdTOcTmVKFNHCHwT2NiXdJULvKZSUYc7NB3/D8e8Y2S1EVfMwfnsMuFz4SOYaIKla3crQViso1nHOikkutpIQ0qFXhMBO12I8A3A9gFBG9AGAQgJOqOqsuAnGxzl4S5bemzEVbx7wctZjwok3+ywuYuPMg4z44VC+hvPs+6upnceFnxkgGfZXNJXzi7tv20TJN3iyeEidhmSLN5lGl2ijyOJcyZRcVc1m9sVPIJu1KLoFtJawW4wblPl4NFlNs27eb0nPI0UgulWyyv7jvcCz+ZBNueOG9xGMZmB/oJ6v6iiWGvz/3Lq549M1Qbr3Y/oR5c1VWU45CNk+HKndx58gJOfREONJzAxpTcqlVKphE5sIYe5WIDgGwM1zFx5uMMbVT+haGJNuA6a5FxRjc9jTMxe8kBNn+UC5ke8OqDR34/h0zsd2A7oHNRUEL+MZ0WP9u6N+9GT1akis2iERgwvYD8PCcZRXNPU3KdKKoxFgolSqSXHTniS7cgBuJ35QjT00WTIJLxX27pWMuuoqRKsmlUtdjBvOaOaLDSKSUNwsk5w3tZkGiqvekKeegU3jujkOZqVvdflRqMYWEWiMpIQ3SOhqVC6OtAWOswBibyxibYxlLgCQxO7W3mLShSFqXqoh6+ZQezSbCaRiqF0LHSIslMRhQL7k4nlrCtQOpx+WXIxKBL+03HL1a01+DiBOveREvLPgo0r7LkN74zB7bhNocCq7/sNFuTbxK41x0945fJh/vnRXrvYzbaq+jvikll7jNiTynDR3Fsjy1RJjen//NW+57MTblJbVphTYXDlkd61By4KgpdJJLTmGvq5V9Iw0ayaBvoUGS5NJZLEVcDn0onm9abzFxgXcWuGotfI6qwFUSVMxFd6UlwVtMdTtEN+W846CzWNJeV2tTwIg4cg5hv5H9U8xeM0/FkJedODbCuESpgSfQLPjMpTzi1JQjfOOQaJ5XzkR5IOtj85ajyXF8G4yMXq3pmItu+TCmflaqfFfmY7HE94Hju7e/FjgvRCSXIPDX9G6r7lWzpFLL0uaSc0jJBHkaHxGmDhe1RK1ckS1zqQBJi/XRucsw+uePYI7CHiPnswIUgWwJG4z1gtqAu0vKp8z8YHV8JwoomYvmUgsJrsii0TqfI0x792M8/WY0EzEA9GxxCT3vL+/t4s8/NlpeOAuoiAEhkCA5YywUK4tzacnn8LMjRyvHB8LR1bkcqYPxAPRIuVHQMXHXwSFbjzwAqW6QKl8ab+fTNp2iauMgOwrkMmQujqOWLFUS58YKGDbgqpPLwSE76W2tcphBtaDVNyTlD2OMqf1qtyAkpSnhuxbRk0xVsVL8TUSJMTTnHCUjAoATr5nqf37dYyJZ2A+VNds1/ZZYOLsvEDam8h2kY/BycymCMyROIFReUllAjD/hIKJIgB9noDKxO2DUAKPaJq1NjlLfz5s2dQabhLzjOg7I2ZTF+Zgi9plVgbfIUkgcOOOTmUA5ajGV5CKvNYcoM4N+3glLt0vXtHljZm/AH96/u5/qJw16x9jndPQka8Qps6+M+Y0B2OJTwJiqAcQiUMUS0z5cVYR+PkfQbX5Uxvoslray2JBmzqLNhd+P7s153yjLjcoORT2/Rg7sgdvPmoDz7pmNJ+ev8Iz9gWOA7+JcJflarD/DIbqT8vHdxJVRb7Hh/bsbMRdVZD0QMFExY23eceAQ8PLCsP1jcK+W1DaDuFLT/Flt3dtV/S0rM/WKiBZD7y5AHxhaKDEhi4FZX7qM4CKIFHnuyoQo8e4ypLfAXKKScCU4ZKdBuPLkPTBeKJttijg7pUlC3SygXQ2anGKZ5BYjoouIaAkRzfT+jvHam4noRiKaTUSvE9FE4ZxTiGgWEc0losuF9hYi+hcRLSCil4hoRCVzSwNTMVuUBAol5gd5jd8uHImsCqKUd3ZJyCLhoLqSnbpfMUKf3w1RfeMzDKIQUx2/XT889ZOJ2LpPK0Zv3QuAS6z5sUCglorbcY4Y0N3kkpTQq8U8qUuQXEjhLWZaqlZXMIqrpsSgP27Ql3H8nttmuhb4CPlcNgk5GQNamlJILn42gvA5f3jibb8wnKmnleo6ddnAs4AooYgMn4gyTfj5/cN3xMCeLbjmS+mT0Mcxl9cWpVeVlwOj1UBEuxHRyUT0//hfBmNfxRgb5/095LWdBQCMsbEAjgBwJRE5XhDnFQAmMcZ2BbAVEU3yzvkagE8YYzsAuArAZRnMzQiqBbvPiPjUFSXGsMkz8g/p2y3ym4iymEuqo9VQvdSqIFF+bGDQd/9zhgKIkkvYvVa8dz88Yifc9c39sftQtwibb3Pxrt2UiZ+w57ZGxwVziBJ+xyEhNYk7fkeBSy5h6CQSjgnbu44IOnrDr2ujqBbTMKymnJO67rr4HGTwZ9Wcc4wILy+QF4fWFOpLP4gyRtoxVZGZSC5ZSr+ihCK+n2lyi5mAjyE7J5ige5P+2d9SxbxxIkwSV/4CwJ+8v0MBXA7gs1WazxgATwAAY2wFgNVw0/1vD+Atxhi3BD8O4ETv82QAN3mf7wIwiapirYxC9VImDV0sMd+DrE+38AJQ21xSXkoWNpcUem8x5T6/9O4CUeP2EseR5iZcVlPOwfgRgUeYLLmYqjO+etBI43m781VILkKcCyd8PM5F5i5xj3rstn1w8vhhAPSPhF/W707aw2/LO6QkhPkcJTIzGTd+ZR/tb3zsppxjJLkcstMg7LpNfKmEVGoxbpOLWd8q25/yOMV6TSu5pJm7WHFVZPikccQoF5weiMyl1VA6TJIiaxF/YzLTkwBMArCMMfYVAHsAyKLY9dmemusGIuLb/dcBTCaiPBGNBLA3gGEAFgAYTUQjiCgP4HNeOwBsC+ADwM97tgbAAChARF8nohlENGPlSrXHUho059PzsFIJ2NThRVxLRjf5cZcYM45SDvqojlqM46hdtw59L5WiwYDdhIUtepKFJRf9+PycnIFaTEQaAsHnJEMkHHxX2unlFpOPj5uXQwFj1alKONPcb/sBPuHOO46yNK9DpExcqUPOIQztp1cZ8mtpypupxXJOvEMGQ1CbxQSBWilGcjFcyqpYEvncpP1JGmIr5scT50/I1qDPpySWt+7ZYuaO3pyg8aiF3cXkbdzEGCsBKBBRbwAr4EoSsSCix4lojuJvMoBrAIwCMA7AUgTOAzcAWAxgBoCrAUwFUGCMfQLgW3CzMT8HYCEArktQLRvlE2aMXccYG88YGz9oUPq0KDJ6K+IOVJMR24qMYbVXL7x/jzCPjlR/ZOkDv6rlLcaxjaTKK5RKEYO+qCrgMRyzFq8JvfBxO0lZnx1P1AKkIW7uHFT1boLx+a66s1jy0r+EkYsj9kRoTXAdFu8B/5jPkZJYLlvblspbLIlYBp59jlFam5wTZa4y0jB3Pr+4e2iqFlOlnJFrxidm00jBXEqlYBMnvp8M5b1//Xs0K9t99axwX3t3MwsoTpJcVHWAsobJaphBRH0B/B3AKwBeBTA96STG2OGMsd0Uf/cxxpYzxooe0/o7gH29cwqMsR96dpjJAPoCeNv77b+Msf0YY/sDeJO3w2VGwwDAk2r6AKhJPU9Z8nDnEH9OscQwZ4lbO2TcsD6h3+T1zRgrw4ib6nAl4giTrPcvMTF1fKDH54b2l977WDg26DfuPskZc+MIg3i98gv15y/uqR8EaoO+mLgykFxKkZ3vjV/ZJ1Fy4bvHft3VxEM8nxP4vENKb8K2ziKaMvQfFjcCJtKuiSt5GoO+z1wMNhlJuPGFhZE2OQo9S4N+iTFfupQzDJjcyx8cvqP/eXj/7tq5qdRipoG0SRutWgR3mlSi/DZjbDVj7G9wjeyne+qxskFEYlGT4wHM8dq7E1EP7/MRcKWWed73wd7/fgC+DeD/vPPvB3C69/kkAE+yWtRohdqXPGkXWGIMH61vR4/mHAb3ag39FlWLlWPQr65aTCVJcQLBX4J8jvDA9w7G9PMnhasrGkouQc14E8kl6FROj7Jd/4SKlw5F7heJkks+rBbjRw7o0YxDdx6cyPR4bixZ2hPHkj/nHFIGuXUUSqnjXOIgxiQVDdKB5B2KVS0xBrSmkBx95hLTaSV2gXVSXrJsmYterWdCeY7eLSB/bpJSvcs4EJYIexumQkpyAmhPKEyWBeKCKEczxuargimJaK8KgygvJ6JxcOnpQgDf8NoHA3iUiEoAlgA4TTjnD0TELZ8XM8Z48ZTrAdxCRAvgSiynVjCvVDD1jhGN/MUSQ3uhiNamXGSnJy+yQomhZ0s6VU8WhjquIvjOoaPwl6feCf2mYi5cmuEvQVPOQc+WPHq25EMqHlOeP6Cnu9OftIub2yuJqHGIuzWV67AMVf5BQhChzxl70Uu5zwkKf55xjgYdhRIO2XkQzjhgBL572A7KY0TCyj/lNUGz7YVSam+xOIjeTiaGc8dALSYTtHu+fQBeW7QaFz8wL3Is95iMYy7XPvuu9re0yNLFh7FgSxJSi7Hgvv70qJ1BIAzo2YybX1zoayuA8IZIlaafo+irZ4N13TPGA7B/j2a/HHmSirIWarE4NvgjAF+HOpiSoYIgSsbYaZr2hXCzL6t++4KmvQ3A58udS6W44qTdcc5ds/zvJmqx9s4SWvJORHRVJb1rbVKrVHSoNN0EAKxc51apVBlb5RT6QLArFJkLh6jLFolYHKHaqncrpp07CYN7uTapOA88rZuvwgCvmjc//VsTR+Gap99xXZEV0eOiZxmnJ9xzZ8SA7li4KgiUBdyXtyWfw0UxpXXDNpdALaZKeCo6CMThm4eMwt+eeSfxOM7AmnKO0YYkR8ku4fLvew7vh3Vt6szGN01dCCCeWGaJLAvQlYT8Z7I0yd/hfUf09z0gO4slnH9PUGlSXLMuc9GPA4RpisrOyzF1ymEY/XO3lo0uIzZHXdVijLGvex+PloMoARxT9Zl1EXx+/LDQ9wE94x3pSoyhvVBCS1MusrsQF91L767CghXrU+eAyoK5cKi8k1S7Z/5+tXjETzxPJFyXnrC7/znpXd+6T2tZBOGgHQa6/SvStcggAj7rZUX+/N5DAbjp3m9+cSGAsDeguMPkTGH3bd3Yj30ViTVNXl7x+gKDvvqV/NXk3TB66174+XFjcPguWyX2nQQuUTZLzEXnZZTsLaZmPkk76CT35qyQZbqbhas2Bu7q0vvA72RI5SmpysX3PJ/TSy47DOoJwK3Lw6ELjpywfX+0NuUwfrt+6N2aT7TNNIpBf6ph2xaPA0YNwKUnjMW1p+2N3Yf2UR7D1WIteSfyIos7GK5KWPxJeEechCxdDFWETqUWkyUXHSkzYvkAACAASURBVIE8dPRgfHqMSxizCkWSVW2f8wIpHUrWsztEOHT0YCy89FhsP8itpHjjC4EKQ5RcHAIG92rF0bttjb9+2dUUH7jDAFw8eVf89KhoUkoj5hIiQC503oGDe7eCiPC1g0Zim76tymMAM5vbrtv09plLU45C0qVO9ZZzkoMtVb+3aKSt9kIJ2/btVp0EmgrIc0vrti7ivY82YLrnqCJK9/+d9SHeXOauHfG65EcaVos5Sul7xgWHYxjPWOEQZl54BP555n5apsH7+Nc39seMC45IZKa1sLlo7zARbU1EewPoRkR7EtFe3t9EAOXn3NiM8eNP74QeLXkcuevWIQ8h8TmXmFvGtSXvJjMUGUyhWMLF/52HFevafMKW1hU5y/oRSuO9yubihJlL3IzlmJhKoYtncNVi8efKBEc+XlaL5RzCNV/eG3sN7+e3/b/9R2CgQlptV5RaeOLHh4SC4HIKtZhJNoK0Th4y/nnmfn6JBlktphs/5yTPTcVcdLE57YVSqridSiHPLU02ARXme0xEZMarN3bikofmAwi/A7oS4UA4CaYIuQ5T3+7NOHCHgVp1F+8h5xCa8/EbAYdqk7wyTuF5JIAzAAwF8HuhfR2A86o4py4LccHqXsRNHSXMX7YOQ/u5om5zPjDgPrfgI0x/72O8v2qD/+KpbBxxyJK5qIhYHMPhNiTxZTlu9yF4YNZS/3uQQTkjyUXh7cXHSRpDdoN1jxdURPmw5JIGbQrJZdSgnti6d6tvn3EUu1sTghvrCWTgM9G3e7Nf0yOfC++cdYxLFUQaGpap17xOzdZeKJaV1qRcyFOrdPlxDYHueVFo46DvR2fQ190bbX0mqQvxWWzduzWUmPTd3x6rn1CGiLO53OTZV86QbC6fZYzdXZPZdTGIBtfQYhY+T7l7FpauafPVJio1RFuh6Ivbad+BLEuYKu0rCmLBGQb3fhN3wr85fmy4TynJZaX4wr7Dw3MRJIAkhkDSpcjMpjmkFks3YxP1ZMjmAvPNRFL0tQm4u3NzjnDnN/f323Wbonwu+X6qpqVTkbZ3lsqSwK46ZY/kgxTI0hUZCHb+ugwD4r1KklxUajHdc+CSy7Fjh4TaZQYVbPgcPPezQ5V9VRsmT/cBIvoiEZ1HRBfyv6rPrAtC3FXods1zP3TFaZ4JV3zBOOHtLKZP+8JhIrn8+nO7GfWlqs+hCuTjRJmrxURVlazblqP5K8X3J+2obM8RJRKvNGqxpPnef/aBOHWfYbHHADG6eG7QNxCR4nb8Mp3iyTNl8KSWg3u3Yrdt+/ju0rpATZMgyjRqsbYyJRfRLVeFHQb3NJpbpeuvI2ZzCMTHu4nSdj7npMoOwLUj7YWSXynV7TMMMT9fpWrUcmEy6n1wk0MWAGwQ/iwkiJJLUj6sIMVIlIAVSyxEZH4Z484qQ8Vc5GR3unQTMlQR1yoCw9u6eXriDUIAW6RoE7e5ZLTeZSbO76HjUGLEuHwpMsFpkpISxmH3oX0xcefklEJiN8o4F4dwyvhhOHGvodo+0kgu/zxzAr6niLM5Ya+huOT4sfj6p9xMTn6urBhvsTjvvZ226qlcG9v06YZTxkeZbrmSS5IhXpfYUZZSdVeyrSbgVYavedAsZHG9yIRfFDIG9miOuILvoXEGAoLr7yyW8MgPPoU/fcHNQqGTXGoSTa6BydMdyhg7hTF2OWPsSv5X9Zl1QYg2F3HNqXYxzfmo2y5fkOJOhiidj75KLSYvXl4/JQmqF1ntiuy28Z0UL54ERJmsX/Y4M8VYGL7NhZLjQiK7WUdmLmnVYumuKSzFeDvNHOGyk3bHb08YqzsttOO/S1BpAYH7KkfOoVCWarH9i/sNjxB48fkeMzZIUpojdSXHbfq04vEffQqTdtlKeY8cx70eGe2eU0taJEk7uiwB8txUj/OrB47EtyaOUp7/0PcODn1PlFxE5sJjX0b2x1cPHImx2wbMY3DvqOfffWcfpOwTCGhMR6GEPt2atNkf+LPKujJmGphEME0lorGMsdlVn00XR0venBjxzMFb9wkMvEHN9kD6IDI3JuvyUonM6uLJuxrvGFUvssomwJnLKI+w7SeoYmSCLZdELhf9ujfhrE9F86f6NheiRI8geQ7yfW6pwKCvHzT4GPYW89q8geKIKH9+OYcirqmfHz8U7360IaQyMXneqkSMPzx8Jzw0exkA97mpJBMiwg6De/nHmKK9UKyK5KKTVk2m1px3tDt9WaJpTzLoKzYaw/p1x4WfGRNq27p3ugTzzYLkAgTrRlU6ot4wYS4HATiDiN4D0A739WCMseh2ZAuHaieqA/es2mNYX0x79+PQOXLSPdNdvuOoM+rGeQPtMqQ33li6Fiqokt8pbS5eW/8ezZh+3qRYtZtfo6XC1X/bWROwy5BoAB7v1THQNSfp4UOFoAyok8kliYeE4lw4cxE6ueDYXfDrB9+I9MEJTLHEIsSNiDDl6HDcTZqaQOLmQV7PaslEODfFM20vlMpyTEiSXHQOEXzuvz95D5QYcNX/3ooc05yL5vm69Wv7YdTgHhG1WkeCWkxsjpMdTBNRBnN0O+abSH7H5Xn7arE66sVMnu7RAHYE8GkAnwFwnPffIgYi4fhAEQj56V3dYELxheRpRwqlkr9YCOZlaHOUXMM7L+nOP7XjQO2xql2i0uYiTHBw79bYBItyYbFyoTN8p4kXkYnhKi8vE0dTCknUFCGCrfAWE3/X2XBEAmuy+08jIejsTA6pvcHEjU8aOxpj8VUodUjK9rt2U6eynT+/E/YaipP2HupHuo8RNijdmvOR96d7Sw5D+nSLPP8gc7aB5OL1qVpCaTJJu+N5zMVjbnxe8mvfJWwujLH34aa0P8z7vNHkvC0dIuGYLqSd54v6jANGuMeJzKXEmQvzd+VXnryH8Y7QhKA2553Q7jkuFkT1Iqt0zGnUIb7BvUJirbvWwOZSuaQhEg7TUstpoKrnIkL3bERJxIRxmHhlcaIqXqcjSS6qeyA2maxT8ZA4FZfOXiZeiyp49ZONHZE2IKoW4zm6Dt4p2FydccAIrY1C9/hNXJGDtEHR49LWIOLX3yGpxbSVNxtZcvHKHP8MwLleUxOAW6s5qc0BOuK2rq2APt2alJl1OXN5f9VGrG3rhEPATlv1Qk/DNNsm9C8vpfGIowfGNpcUjILTwrS8ZdSgHhjWP9B7619qzryS+0wTZGnCQE0uSa8Wi0p0uv7E52KSKTmV5OLdV4cgbULU90tcS3EMeDuvvs81X9rbb9MZo4EgM7YMkSHdf/aBuFRwfMg7FIl7Us0TCDZ5aze5no0Tdx6Ebs25SMYHX80qx6pwG5WJQd/vK3qsSTJSEVwt5ttcoFZ/BZJL/biLyao7HsBn4bkfM8Y+BGDmbrQFI44WhXeHQbtoa7l12iJ/QeuS1cX1q0NTjpQEjBf3EqH0FjNMCaMfnweHpuMuT/x4Iv5walD8S/dS86mUIxlt0yfsuSNeV9o0PDqEVE0KV2RVpmQZIrMwU4slz11MHwJENxaOxltMdz0yHvn+p/DEjw8Jrfdh/fTMRSzEd7yXLw4I4nMAlzmdKjCTBZcco83pJ099lBcP4292vHbZdqFi+gCExJW6+x+cwGLUYt2a0ymBfMmlEG/Q58+q0W0uHV7xLQYAvJiXRYAnfnwIHv5+2FUxLpZRp3qYvjBcQJO/rHFptnX9XntasEO88Yx9/M9NOUdZpGr/UQMi/ZnGuaRxd6wk5YdI3JJ2jOWosY7zMiQHfQWM2Ezlk25MlQQZkhYMzjNhHCZR/0Gci8dcck7ouerVYoLkEnP93ZpzGDWoZ+h4VbE9jj7dmnDZiWNxyfFj8cPDdwq1x0EfdxKe208+vTN+f/IemLjTYOXvHKMGueQuslnh90sbdBo7TR/dmtKVHBjcqwVf3G84bvDeaZ32i9+Ghra5APg3EV0LoC8RnQXgcQRVIC3guuDKnkv/eXWx9nhxQcYRQf6TadCj2JdoqDx09GCM385NtpjPUUhqUBmSOVTePCpilhRg3L9HM76wrxtIx5lLoZQ+B5p4ff01pYN9dWMZkovqFDEoMwvoiLGj2CHrrkFsNpFc0khdfhLSprCKyCH1fGS7TBJU3mW9FHE4/Xo045R9huOL+w3HcEGqTq4pkzgFAO46PEEIVOW9imq1K07a3ffmkhlnoZQU5yJILojWZeHQ5grTwHEIlxw/Frtu40poW3txMsftHk4HUw0bYVqYGPR/B+AuAP+BW8jrQsbYH6s9sc0ZOslFBv9tWP/uGDkwWWAUXwBZQhClkHDeI+93hUShSpeu2gUnFZt69edH4LdeLRduwCwnwaZ4r2SPNL7D9NVLZQhIqp03b0tjV4odQ6MS5R9NEh6GJDgDImLCGDkB5M+3OeeEVERuVujg+LMP3SEyR5N7rvIum/3LIyPHXfQZ86wUomE/baJXeeX2aMnjhL22jRynU4tpJSVxDF8tVrnNRcaAni2Yd/GR+LYU/Ol7kTVyECURXcYY+xmA/ynaLMpAXkNgZIhE5NR9huG3D8+P7VckIrKdpllwYQwtcu+znOIb0BQLExLi8UCyNAs40BmnX/S63dj8Xx0V8UIrhxmogwTd/1kZ9EPMxYkyEpM+wucZqMVS3Av+zFvyTlRyEcY9YNQA/PmpBaHx03qL6TZWx+4+BIN6hT3BtunTig+FzA/i+n7up4f6bvxpn7vv8h9i9Kp1oO43Lot0ZAzFcZUyFwDornh3dS7KtYSJwu8IuN5iIo5WtFkYwjFUi5m8iCJ8Y2zOiSxaLjG0F0rKHbOKeKoSV4oGX85cTGqwB/MIe7ukge5eiZH4aZiBDBUB4WNmlftPr0aKchfdNaS9NCMVCS/hLDzfkM3FCRv0d9q6Fw7aYSDOOTKoSm4yjol3mUod++RPJvrzufEr+2DnrQKfIlG1lFYdVIqRKkygT7kffFZVp+RIG+diikYIotQyFyL6FoBvA9ieiHiReALQE8ALNZjbZou84c7TlAn5x3h99RZcnTl23aY3Hpm7DP17NKtVL4pVqGQ4gsGXI0VS16owFxEVSS4VqsVMhsyFnmnQrvJy03WXlhCmIbhFz2OxOe+EpAPZ5tKcd3Drmfulnpd4iO6eqlR94gbi0J0Ha/tPb2vQSxUmyzop3ioJWXkhynDIzRDN1Zf1QJzkchuAhwH8FsAUoX0dY+xj9SkWJhD1wqZqMZN4Bs4M+nZ3jZAXHLuL7wX27UN3wD4j+2PC9gOwti2IYvb95FPOXWQ8adRiFTEXE+Ll/S/HoK96FpVIQirkNBsLpbeYZkj52s44YAT+MXVhzJjm8+PBec05B0P7BYZ0IlK6TqdFSHLT3NO4DA9JKNeQXa5JTacWCxn0uXSkuGvVMrwTER7/0SFV6dsU2qfIGFvDGFvIGPsCgL5wU758Bm60vkUCfnrUztrfTL3FxAWqMxyK4MSXxwKcefD2vldJziFM2N5lNKoKiKb8wVEQwa0UmV114DaXcsT1NLbacgz6KmLHn09WO0yRQYYSV/pee4i0xfUBAGcePDJ+TIObwdMRHeylA+rneSj29qQXOc6lXDVSSHIR7unw/gEjq6T8cVpirVqHaS5NX89FHEMfoS/fR9Htuho486D4tZIlTCL0vwfgnwAGe3+3EtF3qz2xro4v7bed9jeRiCXVuuZI4xWUdKxqd1xiDC+eexgmjdarHEJ9CJ1wpmWCNV7up636mDMkDhPC4Rt2y2AGca62Rgb9lGoxdZxLlPEnjZMkpZkkiNx7u/5YeOmxOH3/EfjxETvhshNd7z5Oex1Sq/HSInx9weeHvn+wX1vGZCOlQ9pNgCp6nkvXJn0lZYoAzO06l5+0O75/uLr4XVa44LgxyQdlBJOneCaA/RhjFzLGLgQwAcBZ1Z1W10fcC503ZC45Awln4aXHoqcnqfAh01QM5AueARjSp1vES+eL+4XTafg7bBAG9WrxKxiaYkgfNyr702O2SnUeYKYW4y9yOWox1SPjTKUqrsghYh2VXHS6J5PaJCJ2GdILP/m02Y7YcQjfnbRjEFsl3E9Vos20EJemeC96tuTRz4tdasqXf6/Tqi9V0fM/PWo0vnnIKHxGCqrlMJGyBgvvUZKQzjMQVMv+Ui+YMBcCUBS+F1G+ynWLQZxorzPqygjHdST3x1/4RMmFop/95HrSuZccHy5a5dsgCHj5/MPx40/r1X8qTNh+AJ7/2aGYPC4aS5AEE8JRKnEVRFaSi/vfKOW+wWuhjXHikotCVRaZk7RmksYlIpx9WGU74qharMx+YtzwuR1O5aVoinIJtHg9fbo1YcrRo7X2FDGGTFUz6OpTxinXakb7k7IgpzaqBUxckW8E8BIR3eN9/xyA66s3pc0DsWnnhVW227b6kqY6/XRkLO83HsyYRAjFfo/ZbQguf+RNnLCnG62ctP65s8DkPdMzBw7RUJwGJtKD6T1QIVYtlhFl0Lkik/QfMAuidPvJZGpK+GoxJxvjs9iDfE87/AJc5TOXtM+pnMSOoWBlZRaLcNu4YX0BAAeO0pe3AKrrNvzEjyf697dWMInQ/z2ArwD4GMAnAL7CGLu62hPbnCG+pLtu0wdzFRHKgCS5xOzmfuSpPHp3y3vHJtlcgt9HDOyBhZcei5290sdxL2f35hx6tTZhzi+PxDkpJZYskDMw9BZZ+ZJLHPHMyqsniW6axDZFdPc12BHLxcJ0t/eyE/Xlmd3z9KreRR+7dY+26Vv+Lju1zSXGk0sHVfBrHPberh/m/PJIHK5RBddCoOnWnEOf7ukKk1UK7VInolYi+gER/RnAPgD+yhj7A2PstdpNb/OEvCB7KHIrAeZZeb+033ZYeOmxfqBkkndQrPtzHIH1Jt6zJZ+Za24amEguzDfop++/3PvCMX5EP2yfkKZH7Ee1aw7nfVNDnkpWUlUc3DiX5DE/u0e8RKuzuQBuXAYA7OXlwSsHqW0u/EOK08SyESojveq59tS845sz4l7BmwCMBzAbbkT+72oyoy0ApsQgpBYz2LXz45NtLvrf46ZWD4YiwshbzJP8y7K5KPpnKRwEerU24cmfTIwfQxH/IEJlD4vrA6juzjdIjxLOiqwbM1klq7E5Afj2xFF4+icTMWpQz/ImCzPPOBFxqVl04Jsc3TnlqrfqmaqlGohjp2MYY2MBgIiuBzC9NlPafJFzCMUSM1axhNViJgblYJz4fvW/9e2mz8Bcb2cWfjvirq+ocUwwgYqB+EWhqhBEKRIhUhAs3SYg6i1WC8lF8hZLiK7XeVrFSS75nIMRBgla49BaZjqVNPeQ34fMJMY6v1fVQtyT8MO4GWOFGsxlswffVZnSqXD22+ij4qWSOUzrx8e9SN+c6MYa9FPoZ+udxrvJcTBiQHdc+fk9tMdwb7Gs0r/4kktG1z5x50HKdqZgijpmKi+Fmhj0SQ76VMNxCK9ccDiuOln9jEzrv5SLWjBaPm8dczGRQA7aYWDEo7SeGYyrgTjJZQ8iWut9JgDdvO8EgDHGeutPtVChOe9gU2fR+AUQD+OuyE05QmeRYVCvFlz02XBacn54JbvslnwOCy89VvlbLXT7cXAcwtPnHBp7TEkKovzrl/bC3a8uweNvLPeP2XFwT7y9Yn3k3LjLSxM0PqRPKz6/91Dlb8fvORS/e/QtLFm9KdFTiT/Hk8eH+4qqxeohueiPHaCoba/sszo5G1MhMOibw78PmpNMmISYly1tCqaugrj0LznGWG/vrxdjLC98zoSxENF3iehNIppLRJcL7ecS0QLvtyOF9qO8tgVENEVoH0lELxHR20T0LyIyq65VZQyU6oBz/3jThazKIMuzHZcU2SJN1EaVoN6Siwm42+dn9nCLJx0zdgj+7/TxoWP++92D8M4lx0TOjaMJaa79xXMn4Ucx3nT9ejRpxxOJdvfmPF46bxJ+NXm30DER5lIDIk2y5FLmRiNtcbFqI66Qlw6B5OJ+P31/fTYOEzTAbagK6rZ3IKJDAUwGsDtjbFd4DgNENAbAqQB2BXAUgL8SUY6IcgD+Ate5YAyAL3jHAsBlAK5ijO0I1136azW9GA1k4y5Xi5lmEQ4FW3orkNduUJcWdo8xSXJZDhqBGCRh+0E9sfDSY3HYaH0GAF3J3k5FdcwgziO7a//zF/bClycMD1UvDWwu4XG26t0aiZmSp1Jdgz4fU33P0iJUibIBNivlSC45yebyy8m74XuTgiDVzUy7VTbqKZh+C8CljLF2AGCMrfDaJwO4gzHWzhh7D8ACAPt6fwsYY+8yxjoA3AFgMrlv5WFwq2UCrpfb52p4HVr0bg3bLVr8pI1mq0+k5bysKq9doWJQaSWXtIWKGkGNkQV0PDKuomaW9oERA3vg158bG1ucLA6y1FArg34Ww+hyi2WJO74+AbefNcHo2GPGDsHJ44fi/GPNc26pDPqZXMlmxpTqSS52AnCwp856hoj28dq3BfCBcNxir03XPgDAasHpgLdHQERfJ6IZRDRj5cqVGV6KGbhaTCV1/Opzu0XaRNsJL8zF002o+ghsLsmP9Q+njsMjPzg48TgR1TDANhI6i9F7mrVBX4e4zLkyIgb+GjwWN3FlBpJLjLdYVpiw/QC/1EQSWptyuPykPSI59eLAFQPi7EMu5psblygTVWUuRPQ4Ec1R/E2G60zQD24izHMA/NuTQlQrjpXRHm1k7DrG2HjG2PhBg9ReO9VEwFyiv31x3+GRNlEdwnd8PCFenL7eZEc4edy22G6AmdsnV990BbVYEr7xqe21Dg8FocaMHGleK8ZqYpyPqsWqN7fAJkGZ3IO0ZZEbEX4uP2H64jNJqxbrmnchGVVlLoyxwxljuyn+7oMrYdzNXEwHUAIw0GsXa8YMBfBhTPtHAPoSUV5qbzgENhdF1UfFChNdFSds3x/nHj0avz1hrLYPqpLN5Y+njnPn2AA68kpx7jG7aNVIBYHrn7IPZ/ZenEuV7Fgcvs2lLLVYNWYUhkPZP/+upma95Wv74oHvHqQswxAunlfzqTUk6pmT4F64tpKniWgnAM1wGcX9AG4jot8D2AbAjnADOAnAjkQ0EsASuEb/LzLGGBE9BeAkuHaY0wHcV+uLMUFcoSwVwRPVW0SEbxwyCm2dboLq7s1Re0m1vMX8GilddKdpikKcWqyBrj0iudSEuWQjuYhKhWob9G/66r4Y1q8bhvUvL1GqjIN3dLUdQRbyAGG1WDps7WUs7t2ttrm/qo16MpcbANxARHMAdAA4nbmK57lE9G8A8wAUAHyHMVYEACI6G8CjAHIAbmCMzfX6+hmAO4jo1wBeQ4NmbY6zuaigStvf2pTDBcfugkMVRb1M07+kBTd0bw6SSxyKntPEF/aNFlutlWeTiXFenktV1WIZe4uJqDbDPmSn6qi+lQb9kFosHXv57mE7YofBPXHkrulrHDUy6sZcPI+vL2t++w2A3yjaHwLwkKL9XbjeZA2NOLWYCjrD/JkHbx97XtZEICAwmXbbcOj0mKhYEiCIUK+VzSUZ8lxq8VzIyWZzEU55U3F3dYGfWywUsxP8nlZyac47ZdU3anRseak664g4g76InbfqhTeXr0ut5y9WkPokDrwq4UE7xtejaGS8fuGnfXduHcZ72Xj3VmTlrZXkYjKM/Hg5s9lvZP/Y807dZxj28IJM08+LMmFi4tpPm2SyUaA26AdfVOtnS4RlLjWEH1eSwFz44k1bNIkzl6Z8ti/tNn274bmfHopt+nbLtN9awqSWxaRdtsLMC49A3+5BZoVKUviXAxO1WERycQgPfe9gDB8Qb1u49MTdU89HlVus3OSQbn9ujztt1bMm8TnVQFBALmj77B7b4KYXF+KWr+5XcfLNzQWWudQQg3u7bsTFBLUYl1jS2k54nEYllfx0yMoo2ugQGQuQfSXKJJgMo5rLmG2qm+pPzC2mKu1rinKKczUa8gqby+DerXjup4fVa0oNia4pl3YhHLXr1v7nrXu7XiFJNhcuucSVSlaBq32aq+w2uyWh2vnaIuMZHFMv2xe/B5Wos3zm0oWXaOYp9zdTWOZSZfzttL0x0hOTB/VymUuSPT/vq8XKk1zSMiULPUzLGFQ8TjBg8rG1JGrCWuXEtBLJmNsdh/Qpv5RxvWFfLzNYtVgNcNNX9sW0d1ehd2tc0skA/CU2SeMigkeYV0MttqWC85Rqx/iIto1Gwp++uCeuefodtOQdQaIuf5I7DO6Jq07ZA4ft3HXdbv2syPY1i4VlLjXA8AHdMXxAd0x/72MA5mqxtJILjzBPe56FHk6NJBeORrNFHLnr1jjSU+3m/E1PZXM8fk91rZuuAqsWM4PlvTUEJ1BJarFyd4iFKhr0t1SoUn1UE2no1eG7RANpqwk+ty19fflxLnWeR6PDSi41BN/xJcW5+MwlrVqsZNViWYNqpBYj6X8SZl54hF/bp1bgknFzxq7uXQ1yPRcLNSxzqSG4JGKaHiKt+oF3W+0ki1sSaq4WMxxGdpmuBbhNL+v0Ql0NjiKI0iKKLXsLUmNwAhVXlAoIdq9pVTG8kFhXjXxuRPBHUO1dalCWvXEplvVGdJGrcexTV8WWvUpqDL4okwz65bq/8gwAds1nh5q5IqdIuV8v8FvQvw5SUyPBSi5msGqxGoITkCSlGF+zafX8PHKap+W3qBz+s6jyNqwrEKp9R/bHlKNH49R9olmjtyTw99LWbYmHZS41BN/5lZIs+vz4MiWXTR3xCRotzMGJfrUJidMFCBYR4ZuHjKr3NOoOvtFo4EfVELBqsRrC8dVi8ccFHkrp+j/9gBEAgD2Hl5f51iIK02dWKUxVphb1h+OHFNhnFQcrudQQ/u7UUDGWVs+//6gBWHjpseVMzUIDp0ZEn3udV5uJWVQOqxYzg5VcaghTvTo/bnOv/NgVUCu1GN9QWMml8eEHQ9d5Ho0Oy1xqiKH9uuHMg0bixjP2MTreujrWH8eMHQIAGNiruh5Sjs/ELMlqk070UQAACeJJREFUdAT2Mfus4mDVYjUEEeGC48Zof3fIVYukjda2qB7OPnQHnL7/CKNiY5UgZ5i9waL+sM/KDFZyaSDw4DQusFjBpf5wHKo6YwFqZ9uxqBy+Qd8qxmJhmUsDoUmysTRytLZFtuAbCbsbbnz4OQKtx38sLHNpIDR5CQGT0sNYbH7gGwmrx2989GhxrQkbOwp1nkljwzKXBsKhO7sp1P1syFZw2WLgB9ha5tLw4EX/1rdb5hIHy1waCJeeOBZP/2SivzOy2HJgjcRdB71aXRvc/2/v7mPkqso4jn9/3aULLdAXK7j0hbZQwEIUcKUtqAGE0iKx/oFKY8KqxEZCFIxG29RINCEBJYAaJRStJca0BWygaSJNqUUiLy1FC5SX2iWAFCpt5UUl2tDy+Mc9U8fNlu7s3tk7c+f3SSYz99yzs+eZM5lnzrl37qlcyNP65uTSQDra25g8bmTRzbACyAf0m8aoI/zlrz+cXMwawLAh+7GmDdbRh9f/7MEycHJpYD7k0joqpyL7ZI7GN2Zkay850F9OLg3I58+3Hh/Qbx5eRrx//Co1MPlXlC1DvhiilYyPTDUif8C0HF+vqrk8vPB82mtdE6PFOLk0oMrHi9+6rWOYf6HfVI4bfUTRTWh4hU2LSfqapG2Snpb0w6ryRZJ60r6LqsqXStolaWuv5xkraZ2k7el+zFDGUU+eFWsdletV7ffIxUqikOQi6TxgHvChiDgVuDGVTwcuA04F5gA/l9SW/mxZKuttIbA+IqYB69N2U/PUSOuRD+hbyRQ1crkSuD4i9gJExK5UPg9YERF7I+IFoAc4K9V5EHi9j+eaB9yRHt8BfKaeDR9KHrm0jqPSVRlGDvdMtZVDUe/kk4CPS7oO+A/wrYh4DBgPPFpVb0cqey/HRsROgIjYKemYejR4KPm7a+uZf9Yk/v3OfrrPnlx0U8xyUbfkIul+4AN97Fqc/u8YYCbwUeBOSVPp+xh2bp+1khYACwAmTZqU19PmriNdHfnABSyt9NrbhrHgEycU3Qyz3NQtuUTEBQfbJ+lKYFVkBxc2SXoXGEc2UplYVXUC8Ooh/tVrkjrTqKUT2HWwihGxBFgC0NXV1bADhMUXT2fsyA7mntZXbjYza3xFfTW+BzgfQNJJwHBgD7AauExSh6QpwDRg0yGeazXQnR53A/fWpcVDaNSIw1g495QDK1OamTWboj69lgJT02nFK4DuyDwN3Ak8A9wHXBUR+wEkLQceAU6WtEPSFem5rgculLQduDBtm5lZgdSqp712dXXF5s2bi26GmVlTkfR4RHQdqp7nXczMLHdOLmZmljsnFzMzy52Ti5mZ5c7JxczMcufkYmZmuWvZU5El7QZeGuCfjyP70WcrccytwTG3hsHEfHxEvP9QlVo2uQyGpM39Oc+7TBxza3DMrWEoYva0mJmZ5c7JxczMcufkMjBLim5AARxza3DMraHuMfuYi5mZ5c4jFzMzy52TS40kzZG0TVKPpIVFtycPkiZK2iDpWUlPS7o6lY+VtE7S9nQ/JpVL0k/Sa/CkpDOLjWDgJLVJ+rOkNWl7iqSNKeaVkoan8o603ZP2Ty6y3QMlabSkuyU9l/p7Vtn7WdI30vt6q6Tlkg4vWz9LWippV1rGpFJWc79K6k71t0vq7ut/9ZeTSw0ktQE/A+YC04H5kqYX26pc7AO+GREfJFt6+qoU10JgfURMA9anbcjin5ZuC4Bbh77JubkaeLZq+wbg5hTzG0Bl3aArgDci4kTg5lSvGf0YuC8iTgE+TBZ7aftZ0njg60BXRJwGtAGXUb5+XgbM6VVWU79KGgtcC8wAzgKurSSkAYkI3/p5A2YBa6u2FwGLim5XHeK8l2zhtW1AZyrrBLalx7cB86vqH6jXTDeyZbTXk62KugYQ2Q/L2nv3N7AWmJUet6d6KjqGGuM9Gnihd7vL3M/AeOBlYGzqtzXARWXsZ2AysHWg/QrMB26rKv+/erXePHKpTeWNWrEjlZVGmgY4A9gIHBsROwHS/TGpWlleh1uAbwPvpu33AW9GxL60XR3XgZjT/rdS/WYyFdgN/CpNBf5C0khK3M8R8QpwI/BXYCdZvz1Oufu5otZ+zbW/nVxqoz7KSnO6naQjgd8C10TEP96rah9lTfU6SLoE2BURj1cX91E1+rGvWbQDZwK3RsQZwNv8b6qkL00fc5rWmQdMAY4DRpJNC/VWpn4+lIPFmGvsTi612QFMrNqeALxaUFtyJekwssTym4hYlYpfk9SZ9ncCu1J5GV6Hc4BPS3oRWEE2NXYLMFpSe6pTHdeBmNP+UcDrQ9ngHOwAdkTExrR9N1myKXM/XwC8EBG7I+IdYBVwNuXu54pa+zXX/nZyqc1jwLR0pslwsgODqwtu06BJEvBL4NmIuKlq12qgcsZIN9mxmEr55emsk5nAW5Xhd7OIiEURMSEiJpP14+8j4gvABuDSVK13zJXX4tJUv6m+0UbE34CXJZ2cij4JPEOJ+5lsOmympBHpfV6JubT9XKXWfl0LzJY0Jo34ZqeygSn6IFSz3YCLgb8AzwOLi25PTjF9jGz4+ySwJd0uJptrXg9sT/djU32RnTX3PPAU2Zk4hccxiPjPBdakx1OBTUAPcBfQkcoPT9s9af/Uots9wFhPBzanvr4HGFP2fga+DzwHbAV+DXSUrZ+B5WTHlN4hG4FcMZB+Bb6cYu8BvjSYNvkX+mZmljtPi5mZWe6cXMzMLHdOLmZmljsnFzMzy52Ti5mZ5c7JxazOJC1OV+V9UtIWSTMkXSNpRNFtM6sXn4psVkeSZgE3AedGxF5J44DhwMNkvy/YU2gDzerEIxez+uoE9kTEXoCUTC4lu87VBkkbACTNlvSIpD9Juitd5w1JL0q6QdKmdDsxlX82rU/yhKQHiwnN7OA8cjGro5Qk/giMAO4HVkbEH9I1zboiYk8azawC5kbE25K+Q/aL8R+kerdHxHWSLgc+FxGXSHoKmBMRr0gaHRFvFhKg2UF45GJWRxHxL+AjZIsy7QZWSvpir2ozyRafe0jSFrLrQB1ftX951f2s9PghYJmkr5AtgGXWUNoPXcXMBiMi9gMPAA+kEUfv5WMFrIuI+Qd7it6PI+KrkmYAnwK2SDo9Iv6eb8vNBs4jF7M6knSypGlVRacDLwH/BI5KZY8C51QdTxkh6aSqv/l81f0jqc4JEbExIr5Htlpi9aXSzQrnkYtZfR0J/FTSaGAf2dVmF5AtKfs7STsj4rw0VbZcUkf6u++SXX0boEPSRrIvg5XRzY9S0hLZFW+fGJJozPrJB/TNGlj1gf+i22JWC0+LmZlZ7jxyMTOz3HnkYmZmuXNyMTOz3Dm5mJlZ7pxczMwsd04uZmaWOycXMzPL3X8Bdr+nedk8UWEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x2af85a3841d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(ham_eq[\"output/generic/energy_pot\"])\n",
    "plt.xlabel(\"Steps\")\n",
    "plt.ylabel(\"Potential energy [eV]\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJztnXmcFNXV93+ne4Zh3wdFVhEUd0UEFFdI3BWfNxpj4k4ekrjHuJtFTczz5ImJiUtMTNRoNGo0GhF3UVxQUVB2RBFQkB3Zl5np7vP+UXWrb1XfWrqnq3um+3w/H5ju6lpu1b11zz3LPZeYGYIgCILgJVHuAgiCIAgtExEQgiAIghEREIIgCIIRERCCIAiCEREQgiAIghEREIIgCIIRERCCIAiCEREQgiAIghEREIIgCIKRmnIXoDn07NmTBw4cWO5iCIIgtCpmzJixjpnrw/Zr1QJi4MCBmD59ermLIQiC0Kogoi+i7CcmJkEQBMGICAhBEATBiAgIQRAEwYgICEEQBMGICAhBEATBSKwCgoi6EtFTRPQJES0gosOIqDsRvUpEn9l/u9n7EhHdSUSLiGg2EQ2Ls2yCIAhCMHFrEH8E8BIzDwVwIIAFAK4HMJmZhwCYbH8HgBMBDLH/TQBwb8xlEwRBEAKITUAQUWcARwG4HwCYuZGZNwIYB+Ahe7eHAJxufx4H4GG2eB9AVyLqHUfZFq7agt+9shDrtjbEcXpBEISKIE4NYhCAtQAeJKKPiehvRNQBwC7MvBIA7L+97P37AFimHb/c3uaCiCYQ0XQimr527dqCCvb52q246/VFWL+1saDjBUEQqoE4BUQNgGEA7mXmgwFsQ9acZIIM2zhnA/N9zDycmYfX14fOFDeSTFiXSmUyBR0vCIJQDcQpIJYDWM7M0+zvT8ESGKuV6cj+u0bbv592fF8AK+IoWI0tINKZHPkjCIIg2MQmIJh5FYBlRLSXvWksgPkAJgI43952PoBn7c8TAZxnRzONArBJmaKKTVaDEAEhCILgR9zJ+i4D8CgRtQGwGMCFsITSv4hoPIAvAZxp7/sCgJMALAKw3d43FmoSllwUDUIQBMGfWAUEM88EMNzw01jDvgzgkjjLo3A0iLQICEEQBD+qciZ1TVJ8EIIgCGFUpYCQKCZBEIRwqlJASBSTIAhCOFUpICSKSRAEIZyqFBASxSQIghBOVQoI0SAEQRDCqUoBkfVBiJNaEATBj6oUEEqDaJJ5EIIgCL5UpYCQeRCCIAjhVKeAsJ3U4oMQBEHwp0oFhK1BpMUHIQiC4EdVCohkUqKYBEEQwqhKASEzqQVBEMKpSgGRINEgBEEQwqhKAUGmxU0FQRAEF7EKCCJaSkRziGgmEU23t91MRF/Z22YS0Una/jcQ0SIiWkhEx8dWLnv5a2sJCkEQBMFE3CvKAcCxzLzOs+0OZr5d30BE+wD4DoB9AewG4DUi2pOZ08UukNIgRD4IgiD405JMTOMAPM7MDcy8BNbSoyPiuJDyQYh8EARB8CduAcEAXiGiGUQ0Qdt+KRHNJqIHiKibva0PgGXaPsvtbUVHuSAyokIIgiD4EreAGM3MwwCcCOASIjoKwL0A9gBwEICVAH5n72tyHef04EQ0gYimE9H0tWvXFlQoMTEJgiCEE6uAYOYV9t81AJ4BMIKZVzNzmpkzAP6KrBlpOYB+2uF9AawwnPM+Zh7OzMPr6+sLKheJiUkQBCGU2AQEEXUgok7qM4DjAMwlot7abv8FYK79eSKA7xBRHRHtDmAIgA/iKh8AUSEEQRACiDOKaRcAz9ij9RoA/2Tml4joH0R0EKwB/FIAPwAAZp5HRP8CMB9ACsAlcUQwKYhEgxAEQQgiNgHBzIsBHGjYfm7AMbcBuC2uMukQRIEQBEEIoiWFuZYUIgKLDiEIguBL9QoIiAYhCIIQRPUKCPFBCIIgBFLFAoJEgxAEQQigegUEJFmfIAhCENUrIMTEJAiCEEj1CgiQaBCCIAgBVK+AIIliEgRBCKJ6BQTExCQIghBE9QoIiWISBEEIpHoFBCAzqQVBEAKoWgEB8UEIgiAEUrUCQi07KgiCIJipWgFBJEuOCoIgBFG9AgJiYhIEQQiiegWEpPsWBEEIJFYBQURLiWgOEc0koun2tu5E9CoRfWb/7WZvJyK6k4gWEdFsIhoWa9kgGoQgCEIQpdAgjmXmg5h5uP39egCTmXkIgMn2dwA4EdY61EMATABwb5yFklxMgiAIwZTDxDQOwEP254cAnK5tf5gt3gfQlYh6x1cMmSgnCIIQRNwCggG8QkQziGiCvW0XZl4JAPbfXvb2PgCWaccut7fFghXlKhJCEATBj7gFxGhmHgbLfHQJER0VsK9pYkJOD05EE4hoOhFNX7t2bcEFEx+EUEw+W70FA69/Hm8sXFPuoghC0YhVQDDzCvvvGgDPABgBYLUyHdl/1Ru1HEA/7fC+AFYYznkfMw9n5uH19fUFl02yuQrFZPoXGwAAL81ZVeaSCELxiE1AEFEHIuqkPgM4DsBcABMBnG/vdj6AZ+3PEwGcZ0czjQKwSZmi4iAhYa5CEZF5+UIlUhPjuXcB8AxZxv4aAP9k5peI6EMA/yKi8QC+BHCmvf8LAE4CsAjAdgAXxlg2EICMyIeKZeGqLRhU3wG1yaqd6iMIzSY2AcHMiwEcaNi+HsBYw3YGcElc5fEi6b4rly/Wb8Pxf3gL44/YHT87ZZ+SXlu0UqGSqOrhlbzMlcm6rY0AgI++3FCya6rcjzLoECqJqhUQJEvKVTyl7KxJvBBCBVLVAkLkQ2UimdwFoThUr4AAgcUeIBQZaVFCJVG9AkI0iIqnpPUrWotQgfhGMRHR7yMcv5mZby5ecUqHzKSuXMrZV0ubEiqJoDDXbwG4NeT4qwHcXLTSlBBrPQhBKA6iQAiVSJCAuIuZ7w86WK3l0BqxUm2IiBAEQfAjyAfxmN8PRHQiADDz7UUvUYkQE1MVUIYKlrk1QiURJCBeI6L+3o1EdB6Ae+IrUmmQJUcrFypDnKtzTWlSQgURJCCuhSUkBqkNRHQNgOsAHBNzuWJHNAihmIgPQqhEfH0QzPwcETUAeJmIxgG4CMARAI6y8ym1aiTdtyAIQjCB8yCY+RUA/w3gLQB7w1pfutULB8CeKCf2gIqmHLUrLUqoJILmQWyA1d4JQHtY2sNXZBlbmZm7l6aI8SAaROVSDnNPNlmfNCqhcggKc+1ZslKUCXmVhWIh+Z+ESiTIB5EuZUFKjawHIQiCEIyvD4KIPgg7OOI+SSL6mIgm2d//TkRLiGim/e8gezsR0Z1EtIiIZhPRsHxuJF+sAZ9IiEqmHAMAaVFCJRFkYtqfiD4K+J0A9IhwjSsALADQWdt2DTM/5dnvRABD7H8jAdxr/42FREKWHK1UymHukfUghEokSEDsF+H4VNCPRNQXwMkAbgNwVci5xgF42F569H0i6kpEvZl5ZYRy5I2k+xbiQJqUUEkE+SA+L8L5/wBrwl0nz/bbiOjnACYDuJ6ZGwD0AbBM22e5vc0lIIhoAoAJANC/f85E78hIum+hmIiTWqhEYlsPgohOAbCGmWd4froBwFAAhwLoDmtmNmCOTszpw5n5PmYezszD6+vrCy8fZLRX6ZRjnos0KaGSiHPBoNEATiOipQAeBzCGiB5h5pVs0QDgQQAj7P2XA+inHd8XwIrYSifpvisW8QcIQnGIJCCIqC8RHWt/riOiDmHHMPMNzNyXmQcC+A6A15n5HCLqbZ+HAJwOYK59yEQA59nRTKMAbIrL/wAoDUJERCVSzhny0qaESiLISQ0AIKKLAFwKoAuAPQAMAPAnAN8o8JqPElE9rD56JoAf2ttfAHASgEUAtgO4sMDzR0JsxpVLOfrocmSQFYS4CRUQAC6HZQaaBgDM/CkR9crnIsw8BcAU+/MYn30YwCX5nLc5iA+iclHVWsr6Fc1BqESimJh2MnOj+kJESVRAdmNZD6JyKWdnLS1KqCSiCIipRHQtgLa2H+IJAJPiLVb8JCRZX8VSliyujtpShosLQkxEERDXAtgC4BNYs6InA7gpzkKVAgIhIxKiIilPig1pS0LlEeiDsM1JDzDz+bBSX1QOokFUMFbFlrJ+M5nSXUsQSkXYgkFpAL2JqLZE5SkZBLEGVCrlTNInmoRQSUSJYloM4G0iehbANrWRme+MrVQlgAhgGfVVJOXxQZReaxGEuIkiINYCeBXWqnLt4y1O6bCWHBUJUYmoTrqUfbXIBaESCRUQzPyzUhSk1MiSo5VLOcJcZR6EUIlEmUn9KsxJ846LpUQlQrK5Vi7lDHMVOSFUElFMTD/VPrcF8C0ADfEUp3TIehCVi6wkJwjFIYqJaZpn05tE9GZM5SkZokFULuyEuZauhmVOjVCJRDEx6UuFJgAcAqB3bCUqEUQk5oBKpRwahOMYl0YlVA5RTEzzYL1yBGuJ0SUA/jvOQpUCSfdduZTFB1GGawpC3EQREIOYuUnfQERRjmvRiImpcimLD0LmQQgVSJRcTF4fBAB8UOyClBpJ9125lGWp0TLMvRCEuPHVBOw1H3oDaEdE+yOb4rsz8pgwZ+dzmg7gK2Y+hYh2h7UEaXcAHwE4l5kbiagOwMOwfBzrAZzFzEvzv6XI5RJ7cYVSTg1CECqJIA3iZAB3w1ob+k8A7rH/3Qggn8lzVwBYoH3/DYA7mHkIgA0AxtvbxwPYwMyDAdxh7xcbokFULuWo1oy0JQHA4rVb8fRHy8tdjKLhKyCY+UFmPhLAeGY+Uvt3EjM/GeXkRNQXlqD5m/2dAIwB8JS9y0Ow1qUGgHH2d9i/j6UY13GUmdSVS1lmUjvXLvmlhRbEcXe8hav+NavcxSgaUeZB/IuIjgewL6yJcmr7ryOc/w+w1pPoZH/vAWAjM6fs78sB9LE/9wGwzD53iog22fuvi3CdAiAxMFUo5V1yVFpVNZOqMFUy1ElNRH8CcD6AqwC0A3AOgMERjjsFwBpmnqFvNuzKEX7TzzuBiKYT0fS1a9eGFSOgfGI3rlikWgWhKESJYjqCmb8LYL2duG8kLL9EGKMBnEZES2E5pcfA0ii6amGyfQGssD8vB9APcMJouwD42ntSZr6PmYcz8/D6+voIxTDT6hfVFnwpR/CBzKQWKpEoAmKn+ktEu9rfB4YdxMw3MHNfZh4I4DsAXmfm7wF4A8AZ9m7nA3jW/jzR/g7799c5xiF+QmZSVyzliWIq37UFIS6iCIgXiKgrgNsBzASwFFkncyFcB+AqIloEy8dwv739fgA97O1XAbi+GdcIhUhGfZVKOdJeSEsSdCrFfB22JnUCwIvMvBHAk0Q0CUA7Zs4x/QTBzFMATLE/LwYwwrDPTgBn5nPe5iAzqSuX8oS5SmsSsmQYSFaAHTtsTeoMgD9q33fkKxxaKpLuu3Lxq9dZyzbGVucyk1rQqZQBQxQT06tENC72kpQa0SAqFlO9vjBnJcbdMxXPfPxVycsjVB+VIiCiJN27FEAXImoAsANOIlTuHmvJYoYAkRAVislhvHjtVgDA5/bfYpPJlH4NCqHlUinNIIqA6Bl7KcqAlYtJqEzKN5NaEIDK0SBCTUzMnIblPL7O/twbwEFxFyxuZD2IwmFmfLp6CwBgZ1Maa7e0rBVoyxrmWvpLCy2QSplQHWUm9d0AjgVwrr1pO4A/x1moUlCOKKYVG3dg4PXPY+5Xm0p85eLy5PTlOO6Ot/DWp2sx/qEPcehtr5W7SC7Ks2CQrAfRknnv8/UYc/sU7GxKl+R6VaNBADicmX8Ae8KcHcXUJtZSlYByTJR7Y+EaAMCj074s7YWLzLwVloBbvHYrpi5aX+bS5FKK0XxTOoNUOuN8r5QRY6Vyy3PzsHjdNixeu60k1+NM+D6tgSgCosmeD8EAQEQ9ALT62yeUXsoX83Kvf7Iaf3t7cfFOWAAttU8sxQS5g299FaP+Z7LzXTmpK2Xk2Bp4Zd4qLN+wvdzFMFIp7SCKgLgHwL8B1BPRLQDeQcxrNZSEMqb7LkYS84v+Ph2/en5B+I4xEGMW9qIQVK/FqvOtDSms29rofPfrEJas24abJ85zBIhQPCb8YwZOu3tqpH1V9ZSq6VaKgIiS7vthIpoB4Bv2pjOZeW68xYofKkO6vspoMlla6juQTfddugKq/t/bMfzokRn4ZNUWnHVoP+zdu3PJylMtfL2tMXwnjdIJiNJcJ26iaBAAkATQBKAxj2NaNOVM992yx9+F0ZJGyGVZMMi+ZsZjfFUCo4UrXRVPqTP8VkqEZJQoppsAPAZgN1jpuf9JRDfEXbC4IVTeiL7U6M8vXSEvRKGkxQdRUgodkJTKctCCxkvNIoo2cA6AQ5n5p8x8E6xEe+fFW6z4aa1LjqbSGdz2/PzYr/Ph0q99X0LTaDjdgt6IoCimuEbymYBrCs2DmfHQu0uxZWeTsy1fQVzqd71SBgpRBMQXcPsqagCUN3ymCBAoNrVz0/YmLF1nCKcrQqN5bcFq/PXtJc0+TxBTFq7BmX9+Dw9MDb6Orka3pBfCqVdDkeIqprr/cpsWPvpyAzZrHWmhMDOe/mg5djSWZt5AEO9+vh6/mDgPN0/MDozy1VidZSsrwEn97xnLcdljH8d2fp0oAmI7gHlE9Dci+iuAOQA2EtHviej38RYvPuLUIE66820cc/uUwGsXyvYSvLArNlprRC1aY85bZFLTW4sGEReqQ/B7DKUwbTSk0vh/f3oX339oerPPNW3J17jqX7Nw66R5RShZ81Cag65BtKDxiJE4y/eTJ2fhuVkrwncsAlFyMT1v/1O8H1NZSkqcM6m/2rjDuL0Y12tMxT8FpRAB5nXO+sHMaExnUFeTzP8iEckm64unhvUJcoqsgCjeNZvSGWzZmUL3DtHmpSohPXv5xoKux8x4ed5qfGPvXtjWkAIArN5c/jQqKfu+arQFFvI3MRVWLw+8swTHDu2F3Xt2yOu4lqRRN4couZjuD/rndxwRtSWiD4hoFhHNs+dQgIj+TkRLiGim/e8gezsR0Z1EtIiIZhPRsOLdprF85ZsH0YzRZEMMAuIf73+B1Zt3Ot9V6fyej9EHEfFh/vODL7HXT1/Cyk1mIVoM2PO3ufz0P3NwzG/fcL6b6iAb5lqkiwK4+slZGPbLVyN3bo4fpMAyvP7JGvzwkRm46/VFTh0HXXv91oaSpK5IpW0Bkch2V4VqrPl03I2pDG6dNB9n/vm9Aq4Tfd9Fa7bgf15cUHbzpIkoUUwnENGHRLSGiL4mog1EFGXRoAYAY5j5QFjJ/U4golH2b9cw80H2v5n2thMBDLH/TQBwb/63E51Ck/Wt2bwz1s4tDL8X8uQ738ZTM5bnfb4l67bhZ/+Ziysez7Vp5uOjifrCPvORtR7Dsq/dz/C9z9fjr28txiPvf4EnPmxeKhJ2/AHNOo3DI+9/iaXrszN2jQLCp3duzvKnz860zAhROxtVB4XetppTsHxDtm6CznXIr17D9/42raBrNaYyuODBD5y0LUE4GkRC1yDyu57aPaqmC2SfZyE+nXwE0fkPfIi/vLkYq7RBWkshionpbgDfhuV7iPx42XpLlRG71v4X9NTGAXjYPu59IupKRL2ZeWXUa+ZDoSamEb+20iss/d+TQ/dlZjw7cwVO2G9XtK1NFqXD8uuc5q3YjKufnIUzDumb1/k2bLc6Bd0ZmR095u6/aUcT5izflPN7VGHbZDAXAMDZf3VbLs86tH+k85lwOgOtTKbZ3w2pdEGmroZUrpAO80Hk0zEplJ8sw4xkBK0z08xQqoT9jJjZ0XLDqnXGFxsKutb8lZsxZeFabNjWiGcvPcLZns4w9rjxBdx86j64YPTuALImPdVmdjal8ZN/zSrouvl03E2FVJpNPoNPR7AbDkmlM6hJlm/qWZQrLwcwk5mbmDmt/kU5OREliWgmgDUAXmVmNdy4zTYj3UFEdfa2PgCWea7bJ+J95E3UJUfH3f0O/lPgKmTTlnyNK5+YiV9OcoelNsdJbeqcGg028ahs3WnZmju2zY4VnM7BsP95D3yAD5bmKpA3PD0n0vVmLbPs47WJGBt9BFPLGwvXYK+fvuSUJx8amoJMTOaLNscmHfVYZeaLqq08O/MrrN+a9TGoKkkzO3bGuIweqlNMJNwvg2rfv3lpobMt64OwCjhp9kq8tmB1fhcswPyWThcucPPRcPz6g6c/Wo7BN72IZV+XL99UlLf0WgDPEdE1RHS5+hfl5LYwOQjWBLsRRLQfgBsADAVwKIDuAK6zdzc9ppzHTEQTiGg6EU1fu3ZtlGIYURpEJsN4d9E6/3WMl2/ClU/MNP4WhjIHfVnECjbJAlOHFZVNOyz1uVNdbXZjgAahd6h6RzT5kzWh19LNY/HKh/C3882FVtspZAS8sxANIo+eKZNhnPnnd40r44UdF2X/9VsbMGXhGlzx+Exc+s+saVFpEBk2v4zFRD2PpKd3VMJAlxuOBmFvrEnkXzr1SPIJj22OBhFW3yqMWB/weY94ce4qAIhkhouLKCamW2Cl2eiKArO4MvNGIpoC4ARmvt3e3EBEDwK42v6+HEA/7bC+AHJiuZj5PgD3AcDw4cMLHuBYPgjg4feW4ubn5uPP5xyCbu1r0aNjGwzu1Uldq9DTAwDa1lrmi+Z04GEws1GriMr2RkuDaFOT7bEdJ3WRx4/KnAXEGwZoimJy/BKwOlL1vRBtzqhB+Cw5mjV3RT//zlQaHy7NCq58NYgwTrt7qhNpt07XIBwBwY5JLi7HqaNBeCqgyTahqu3bGlK4+TlLA1dO6mQBAkKRj6BWzvFC3oMw2bJw9RZc9a9Z6NS2VgsKcV+nQxur//ALbWetnuIiyjiuFzOfxsw3MfPP1L+wg4ionoi62p/bwUr29wkR9ba3EYDTAajEfxMBnGdHM40CsCku/4N9fTCz45D78uttOOu+9/GN37/l7BPV8bpm807ji6Q6XTXijONlY871S7yxcE1ktbQpHVCmkOLmezvrteynqWKG+3gIimK6d8rn2P/ml7MTpwo4f1AUk98zidKWvt7WiIHXP4/3F7vX2NiyM4V/RwhAyHZowehh2HW12S7A7YOIF9VRezVJ1R5Vv7dxR9ZBrHwQtR7/lXFSqg/F8g2EESaI1CBje2PK6eS9QqVdG2v8vs1HQJRi6lEUATGZiMYUcO7eAN4gotkAPoTlg5gE4FEimgPL6d0TwK/s/V+ANUN7EYC/Ari4gGvmBSP7guwMsCsH8cmqzRjx68l45P0vcn5TDcw74szn5du0vQkbAjJWpg0axIUPfojj7ngL339oOm6eGDzRSanvDOCKxz/GpNkrIo9KwmZae9m4Pfuyp5uhvpuY+9Um3PTMHDBzjmnm1ufm4/ZXPnX23daYdgRUISOwYCe1udFE6ZiU+e7u1xe5tt/0zFz85MlZ+OjLYHNYIX4O3UmvBuZ61TRH4AWhrpGjQdjtUfkmdGUha2Jyd1veSambtjdh8VrzJM98it2kvRs6v335Ewy8/vncAzTCqkK1P31ek1cDbK80CHtOipdSTE6NYmL6bwBXE9F2WNlc7QhR7h50EDPPBnCwYbtR2NjRS5dEKE9RIDtbX1v7BQl66YNQs43fW7we5x420PWbqnxHg/Acu2VnEzq1rUUQB976CgD/qKkMs2tEq0wdO5rSjiPv5tP29T2/rkE8O3MFnp25Ar8780Bjeb2YJlG99ela9OnWDnvUd3Rtf/PTtfjhP2Y435vhV3exaUcTNm1vwnf/+j4270zh6uP2ckwCqv5Mguyf9qp++cgHpdKbndTBAmLh6i0Y2LMDenasM/4O6FE67vOv2myN+Lf5dBSKtI+ZSzF7+Ub8ThOUAFCnmxbth5Fmzkay+bSCy7Ww6DVbdqJXp7bG/d76dC0+tIMafnLcXs52xwfhMRepgAvlm9AFiHJSJ5PBlXbyXW9j+YYdrncmm203DxOTz773vPE5AEuA1CYT2LS9Cbe9MN8lfML6DvW7HmDiHTS1tQevfpNjSzEZL4qA6Bl7KcqAlYspq0EEvfQK0xwE9VImDV5XVfneCiYiLFm3DcfePgX/960DcNpBu2Hhqi04sF9XqyypNFZvakD/Hu2dY773t/dxzsgBuWXMuE0eO/KcuGRyxEWZJOXHeQ98ACBXoJ1vb1ekiqRBjLv7HSxdvx3tbH+PnkIlSunz0R8yDCTJx0mdcV/zq407cPo9Ux3t76Zn5uK25xdg/q0nYM7yTejSrtZVv0B2ZOw9f2QnteZnMXH1k7Pw6Wr3yFoXEKqzNoW5fr2tEd3a1zpC5PnZWevviNsm44ObxrqExOadTejcttZpD4BbQKjRsq7BXfbYx04HTpo/RFGTIOxoTOPCBz/0ewQA3PM4vHj7/NWbd6JT2xq0b5PbFTomO58KULPcj/3dlJx1KcI677RJg8hJFW/9JcqG/171zT1zzhEnUWZSpwGcCeA6+3NvWBPfWjUJez2IOkeDyO2wvBVgMiNlZ3nmdjVewaC3mS/WW3bTpz9ejhufnoNx90zFqk3WRJkb/j0HR/32DWzVRoxTF63Hjx79KKexZphdwm1bY/Ao06/8+nlnq3kOeZ3JH5ONOGrj/vULCzDw+ud9X1I1gU2ZA1KZ7Jg3UsdKhEyG8cKclS6HrQkl1EwJ7LLJ+qzvE2euwNotDa5RqHI2nnr3OzhKm5mtULZ172Al6gxp1cF495u6aB1ueW6e8XiTiSmdYddcmPVbGzDsl6/ij5M/8722bj58Y+EaHHDzK47moOOYNJ0opuxvz81agefnWIKHtLIoapKET1dv8S1DENk24X4II3892XemdNggRuWGMi1aFNa8dQGh7tV7vYwmRNXg9E9TFuX8HidRZlLfDeBYAOfam7YD+HOchSoFRFYlqpfCpMZ5K9k4ysi4Q/B01DmZgdPvmYrXtVBQVblbG1KYadueVef+9qJ11neDScHbsZ561zsuQbK9IU8Nwn5h9fP+/d2lTrmLgSlxYVQBcd9bVuLgsN1VR5xKs1bw8GsQgI+XbcTFj36En/0neKHEBSu3YOD1z+Oh95YCAHbpnDUXeU1MXkdqFJQJxWvu5BDNQKE/03Pvn4a/vGmZQr73t2l4cOpS4/G6iccU5spgJ/pM1YUJ/TyPvm/OoQNxAAAgAElEQVSZ71Z4cpIt+3o7Bt/0Ip75eLkjzLw+iGxZrL96n0kg3/3/+tZizPjCP8GDahLqEaUz7AzI5q3YbDwmLJBi8w7/wViY9m0SEDnvhFZWJSD0sOAiu/GMRDExHc7Mw4joYwBg5q+JKFr2sBYMkZXuW9ngdQfRiX98G5ePGYxRg3q4jtFDQRWm47O/WTWoCwGFmn+gJqoB2ZeyTdLf9ui9zuJ121wOuXw1CFX+UmSJ1clXPf5k1WbsbEpj1rJN6NutHfbu3Rn9urfP2S+VyWjnDu+kieAI2C/WB0d+PfORFUk09yurQ6nVZrh6J8qZ2gpg2esVX6zfhgE9skngVKfo9UF8ssoaNYd1OvqI8u3P1uHtz9ZhzlfRY+hdvhvt0W21Bx2qjZiSFeod19fbLE2ss8e/Nn+l9dxenLMK37Jn/HsnynnPp7d3Bvv6jG57wVqfPSzDgXpG1zw5C08bJsBua0jhjYVrcPL+vUOjwoJScOjNe2dT2gl5B6y2ryaW6j6IVZt24gAtEYJqx6kMO6ZjXRCXYpGuKFFMTUSUgP2ciKgHCpwP0ZJQ8yCUBqC/XAtWbsZlj32cUwG6E0m9rKoSn/4ot7GpDt7UGSqVfGtDymmAanSkRp8mf4LJyaZv2panBqFedtO1GNYo8P53luD/XvokdG2Ad23NR2fCw+bU0/kKiFPvegffuvc93DppPib8YwaO/L83nJQfOql01mnfJsIonkBosO/dm/7Dy0PvuU2MevPwmpj8RrojbpvsfD76t1Nw3B1vIp1hbGtI4bS7pwIwB0wA4RqEacQ7SfMVGB20WjFVX/Xp6i2Oz4g5O4hRTtNlATZ+wD1K11Hfa5LklMVvSoPjMHe9c+FBBfpEzsZUJqvFewIXTMIBsEw4l/7zY0xZuNYoCHWCMiur67z56VoM/dlLrgmZc77a5IQaN6Yyjr9nwj9m4CV7chyQrc9Mhp1Bg552o6w+CCJS2sU9AP4NoN7OyPoOgN/EXrK4sWdSOxqA52GnMuwSGvNWbMJD72Y7CHVck9aIvC+g6qiaDA1NjcY2bm/KsSWqkampQ359Ye6MZf0l2rQjOLHYlIVr8MA7SwAAG7Y1Yov98psc8M/NWoFz/jYNv5w0H3+a8jkenZbrg9H5rpa4jdlSi1+Zb06JkM4w1mzZiev/PTvSRD/Tu7B4XW4oYyqTFRBRc9io/fOdgKXXm9dPEDXL6aert2LRmq1YuSmrWfjOTYlo1/bDGJCgfVb3s3pzg1MGRlbDUoEApvpKu55FduTrur6KUEoknP2TCcIJf3gLj3/gTtDopP3QipzJZJ3nfoy7Z6rzefivXsWht73m+p0Z+DJEUwSAj7/c4JTfb6DekErj1ufMqzuqZzDVHjR969538aAdTacLlgbNxAQA05Zk58AoU29KNzFpbbQU2V+DTEwfABjGzA8T0QxYE90IwJnMHGysbQWQLSFSho5eob9PJ9/5juu3VCaDNki4XsqLH/0If/hO1n//xIfL7H1zK1JXH1MeM5USECazjzcLKuAe/QfNmQCAC+wIkFs9+aH8TEwrNumTqqIntrvssY/xqo9wAKz7vmXifDw/ZyWO2rM+8nnDSGey80Ki9PdE2RdW9yNt2tGE52evxNkj+vkdCmbg5Xmr8IN/zECfru0AZDuGfFZiS1C0dT5Um9HDoxtTGXzj92/i5tP2CQ2ZTvkIHmbGrZPmOwLA/WOugDCdRx8cqd/9ogBrE+TcCxHhk1VbcL0nl5cSBPr7leH8ZjVv1sy3WR8EGwMEFF3bWdZza65Mtk4aU5kcs+H8FZt95wKZNMn731mC/zesr8u36M2jZuyHWDMxUWlNTEECwikJM88DUP6lpYoIkVU5d7xmxYWbJ8r5V4AaYemd/0vzVrk6RWVzNV1bV1/VKF69ZLV2Q9zRFM2foHfuejqLfPDr0PTn0q19cAeko5s2TOjzNwpNnWCqnqZ0xokCiqKB6zPR9XL89uVP8Mj7X6Jvt3a+x2aYMdFe2UuZDFSb2Z5HuLFqi2GkMhk88/Fy/PiJWXjpyiMxdNfOWLlpB778ejt+MXEefnvGgYHHmzofwLr/B6cuNf7GYCcdixogmM5j0iC8Go1qS8kEhUbgqPDbjOe8hTpmvU5qE/NXbHbWRUlrAzcAOOa3b+DdG8a69m8MyEKQDX/PbmtTk8CBt7ziPkcq49KJTMI3lTZrEOWeKFdPRFf5/cjMrXa5USDXfbnd4NwNqoCGpjTQrjbnZYkaeqYLFjVCU9vqAjQIE7o5Y8P2XBNTKp3B4nXbsOcunSKdww89sZuiNknB6Tp8SKXZGaEVEvEDmEeT6Qw78wii1EU6k9HWG0jgPx9/hTsnf4ZRe1gBCksC0jiw85+2zf6ez8p/TWmOtP/Fj36EUw/cDQDw96lL8b2RA5wsvAmi0ElgfvXkJzgA637U4EHVk+m9cI/0zQJCaXY1yUTWtOtTJjWLWD8HMxc0f0af9RzUJk66823nc1M647r2ik25azUEmXOdGeHaiN9Ux1YUU+4+3vXeVR24VtUrgSc4yEibBNARQCeff60ar7PL1BkH9S9r7Zh57/yJKJ0SgQLNTu3slyPMn6DQy26KrLjz9UU47o638FlADHmhUUyFCAfAGnEGTTKMgukF2dmUdjSIKLJa91kkE4RrnpqFxeu2ORE4a7f4z41g5pz69nPQBtGUzjjzYoKwJupZDffxD5fh1LvfyeY0Igo1OZicrkThAl6ZN2oSCWQyjM8Ma5XrdeGEHHsFRFPWlJf22Ueh3gG3ZtL8UfMUgw/PRCrNoc9z0w5/bd0kIExlb0q7NQilSeq7pjIZ5/1s44qcK68GsZKZb429BGXC6+wyzTlYsMpsIgKsmPjde3bIGXlHleomVVI1oE72qDDMn6DQzUOm+5hjr1EcFMaZ7wzs5pLOsPMS1RpMTJkM+4ZAOvsYXpCz7nsfp9mj7GgaRHai4ZufZtPH77A1yiCbd4ZzhZCT0iGPl/fpj75yHJJheJ+Jut6Sddtw2/MLAo81CYIMc7AGgaxpKMOMO1//DH94LXfCnKsjt9ux15m9M5UdBSth5ZeTS5lSMh7NpLkC4rEPloXvBOt98F5qwsPT8c19dnG+bzRo64pVm3Zi0ZqtLpOQaTJuYyrjMmekDEEz6Ux2Up6eXLHF+CAqEW/fYxpB/0DLHeTl6idn4c9vfo4Ru7tTUkWpNKJgO65qVOujCoimYAGhRmPFXJeiuaQyugaR29S+2rgDc0Ni+F+eZ3aCq44pSmeSzjC2G3w9KvY/KGqGmXMEiJ95JYi3P4u+ron3UemXUfMl/DD5OdIh5i0VjQZY9+TNNOucR+voVTv2+vX0GeJKc1hpMN2oa+n7AcXRIKKyoymdY7J7Zf5qV1RekIBQKcp/cNQgZ5vJjO3npM54BO4W+71uq818zyevVKEECYixAb+1fihcg/Cjc9sabN6ZwqI1W3FA3y6u36KEnt3/jjnywftSrAkwb+joGsRWw32oSTqfrSksTUEcZDLmSYqKB6cuDc0W67eqmBqpRXl/lq7fZhxVPv2xNSmOyJ3fSSfDuddg57foL28+5j2vMPWLTIpKKsOhM4ZV+0pn/J3E7nBU669Xu3bMJ5oD2E+oqXPonaDlg4inU/R2tjsa06HCKIoJWH8fTYEwDakMFq/NmhcbHc0qe+0npmfbZ6KlTJRjZv956xVArpM6+kvat1t2Bq9XbWxOllJHQNgniSq0trtMTLn3ocIT9fUYys0fXvvUmV1u6uSCoofCyPogOFRg+6VZUIcxAz06mDOw6qnFFaqjyWekaxLqfvilxy6USCYmWyNLZfzt8rOWbcTA65/HojVbsmnuU36j4/A0FsoZ7Y2OikuDeMQzx6cxnQkV8iaNwEuYAPemI2lM5WoQ7vPp865CL99syrcadpnxOqnzGZm47IqeUdKNz0Rbm9lE2tO5RI2l10dqplQbSoOI6vQuBXq0lSkyxWSvjYpKErd+W6Mr/5WJsE6gKZNBuzbm14Q5V2NUX/Ppt6MOTvbdrXOuBtHMXiKVCTMxuTUIPVWIznOzrXDfl+etdjp17zoGjamsxhg2S9n7LgDxmph+/qw7ij+dyQ1A8OK3kI9OmG/P+447Ob18Ho/uR2oRyfoqlbAZmUEkQhxPUWlb6378qvGrRhAWejqop5XHZ3uIk1qFxukTh1oSJudpc5zmurAf/5A51YcirH9Npf3NKlsaUjkZS7M+iOjtImqnl85wjgbhXVxoWP+uka8LWCPScA3C1mgbU8aJmkB2cmdDU9Y087VnTk6TbmIK1SBynf0ZZiwsMJtrvsz4YgNmGVK56EQJTQ7TMrZ43kn17JZtMPsL9bpqEem+C4WI2hLRB0Q0i4jm2Wk6QES7E9E0IvqMiJ5Qif+IqM7+vsj+fWBcZbOuV/ix+iAuakoFEx082WGzPgirEYSNLFXHr0eLGHMx2e3Iq+20FEwNvTnPNR/CRuCpdCbQTOUVumrPZroGjDDnmpjeWOh2cKuAhKg0pTMhy84ydmppYfxwUpWnMs660t4oPD03Wdhzz2oQrqLgty8vDDyumKhFpZpDmPnQ66ROZxhvf7YWp9z1Ts6+e/fu7BYQrVyDaAAwhpkPhLV+xAn2WtO/AXAHMw8BsAHAeHv/8QA2MPNgAHcg5nxPzQnR0uul0HkAANC+zv0yX/j3D7F8w3bHbhkmINT8AX0kY2qQTvqHFiogTCPYfFJVNAeT41CnKcDubiJrIijuy9upbQ3SzKHpQ9rVRknQnGVnUwZn//V9398Z0drNR19stM+Xdjo976RNPbV8mAZhNjG5j1HhzC2ZfOcXpZnx7MwVxt+6tKtx9TetWoNgCzWjptb+xwDGAHjK3v4QgNPtz+Ps77B/H0uFLBgckeacWW+o6Qw78xbyxatBANb6ycpsFTaKVnNmwoSUehlLndI7jI511v2bHHmlEmZBKZsB4M2Fa41Lq/pRyES5KBy9Zz0yGQ5NS9I+Tw1iR1NwtA5zNG1OCYWdTVmTlTfti9rn0zVbQp236Qxj0uwVmL08m53VKyCMuaNaGPkEIACWxuQ3/6mLJ3NDPrP1C6Wwni0iRJQEMAPAYFhZYT8HsJGZ1VNbDqCP/bkPgGUAwMwpItoEoAeAdZ5zTgAwAQD69+/fnLIVfKz+QjWlM9h71874wLB6Vhh1tUkktVmlgBVr3buLtXRjqAZh30NYHh81Sm5pGkSHuiS2NqSM5oZSlTXMcf/VxuDU1l6cNPBFVv+TCWumdNjkwXw7zTWbzU5nRUMqnVcI+PptDY6Q9B6nOrS5X2121tTwozGdyUnt4pVjXh9eSyTfBbwyGfZt+14BUQozbKxPmJnTzHwQgL4ARgDY27Sb/dfU8nPeMma+j5mHM/Pw+vrCs4A2T4PIft7emEbbPEdtipoEuabOK9TkIdMiITpRk9wpH0UpRhz5oDSIaUtyhevOEmk7xezHO7RJZhcOyrDvokGKfNpgkqyBRDLkoHx9EGGBC5+u3mrMQ+SHPltfhRv/0c5wnE9Irmldaa8vqG0ewvAv5x6C0YOzC4CN9ExwjQvTcqRBpNlfQLRvU+MKimlOgExUSiKCmXkjgCkARgHoqq010ReAMrgtB9APcNai6AIgtrkYzYli0hvql19vx8YCM6jWJCjSSNkvmZ0uILoGZFqNoyH5jVTrQjpFIJtPRnVmpsWWgpLktVTa19W4NIi2Ic/CtEytH4kEgTl8USNvp9nZY/68+dR9Il+zEHSNS82fUAMBP1OoaZBkwrtWdz7p59skE66cX0cO6Rn52OYQJUuvTibDRv/bQf26oiZBHgHRijUIIqonoq7253aw1pNYAOANAGfYu50P4Fn780T7O+zfX+cYV8Qolg8CgJMiOF+CXvZ9end2Pusv0OF79MC+u1m/6REtPTr4rwLrfbGKQZ2Pel8b4WVXL01djf8LvrgAAXH2iMJNjsWgbW3C5YMIG+Ey+y9N6iVBaiZz8CvhHUwM8WTwPWO4//oWQQzokbu8qwndLKoEgnoOfhpsh7rwjr4mQTlzfKIMRpzjk+Tk1wIKTxAZNymPiUm91/v16YzamoTrGYYFWBSDOJ9SbwBvENFsAB8CeJWZJwG4DsBVRLQIlo/hfnv/+wH0sLdfBeD6GMvWrCgmr1MvSqdoIqiRdtRGfnoncnD/rjionxXrrguYbu0DBETISCOfkazC757zOZXJx7JLZ/Os5SjU1SQcTSrfW1JzSpqDEnhsz/gNFRCIPnpWPogg38YPjhqUY3as7+h+noXUNQCceUhf4/axQ3uFHqt8BX4mpigTztq3STrt5bsj++PSYwc7WuwZPmXTqU0m8OHS7LKf/boXPlM/TtIeDWJob0vAp9Kc886VImw9ziim2cx8MDMfwMz7qcywzLyYmUcw82BmPpOZG+ztO+3vg+3fF8dVNqB4PgigcAERVATdhKN3Iu1qk871dAHTNUBAhDmz8rVbe8uk07NT9A7eO3o+d9QAvHnNsXmXRZEgchyjPTrmJ2gKeQZe1DNRM37DnKjM4X4KRUNTBmu3NGDhqtxU24obTto7x0fhjWryWys7DD9hFyXYQ7VXP61wz106hp4jmSC8/ZkVrzJ8QDdcffxezrPTZ2Xv5bPmifcdPXn/3qHX9KI/yw45zzXv0znoC3F9tXGHKwebalOpDGOlJ2CiYnwQLZHm+CC8JqZCV0QLQm+MtVon0rY26fymX7Z7B38fRFgOJlO4rR8/PdmKM/Dr2C4bM9j5fPmYwYG+Ee9otm+3dnk5Hr0kE1mzhnfkHEZY9M//O7hP4O9A1uym1onwM6HdcZa18hsRRV4s6d3PrSyqfgkKFUnP+byCryZB+OXp++Uc94OjB+Vs0/Grlyh28LD3Y8xe4VqIPqdClUV1nk3aiG2vXf0EBLnquJAoxi7tsm25needCbrHobt2wikH5C+QgKxgy2Q4J2S2tZuYWjRR28fjE0blbPMKiELV9iD0xqz7F9q1STqdrq6K9u/e3veewuz53gl7QahV6fw6tna1SccpeVD/rphwlH/H431szc3UqY+Oo2gyavc2yYTLXGca+e/Xp0vONi/ZZTKBD5ducNXHGM0UU2NrfoToPoioazF726K3Y08kCKd4Rs9jhvbCDSeaAgyz+AnQvTVfmaJLu1pX3dZomm7Pjm2wXx/3MX6m1m/sbRYcqiy1NfYKd5rzu0OdebBTm0zk1Gu+MkJvX16/SZCAGLl7d2clQBNBzd5Z5tUgIFq1k7pSMDnCvGH7+YxG/nzOIdHstvbIr0u7WtRrnV3bmqSzsLreYPr36IDaRGEvQD6Tq5T/xa9jSyYSTueQIAoMI/U+t+ZOLtPnCPTs6G9yU3S3zXKN6YzLBLHvbl0wdNdOuGj07lpZw6+vNIbpX1jBd7OXb8KvTt8Pb197LB644FBnP124RjVPRg3X8Ha2pnrympmimJ38NIhj9qrHOaPcwQFj9+7l6jD1z+u2NmL0YHcEkV+wht81sxqE9VcfWLTxOVdtknJEbL7mNr3e2ns0iJoAfyIjeFndoKR7SktKZzIuk9LuPTtgrI8ALSZVKyD0zuma4/fy3c/0gnk7MmbGKz8+yrXtrrMPNp5vUH0HHL1X+PwN1TB26+p2pukahJ7o69i96nNC6qJ2Pt7GHoS6d/3c3x2Z7SBqEuR0CGHrJHsHXSYBoXfSlxy7R2Ckkn4+tWSo4qT9d8V5hw1wbdPNX/r9JBOEl648Cj8/dR+MP8K6vnfmr5rMqKPayibNHHLOqAHo190dAaR3JqqeD+wbrKHoj6auJuF6Lu5zux+qyVdEPs3iD2cd5Ht9U0bbRbediMP36Ok8uw5tkjhkQDeccUhfl1Dwlsm7DrVfR62fQ9fAlCagOl19oqXfYK02mchpi6q9dQswg+p00zR576AqSINgdtf5tSe4+5ugd0S1qXSG8dszDnC2n3fYABwyIP65HNUrILTPQWsPmDpZr8RnzppeFKYOROE32nj4ohHada0Sdu9Q6/KXtKtNopPd+ekzVTu1zW3kUSNk8tEgUgYB8Z1Ds6GTCY+ACIq68XYMJgEx/shsR3jN8UPRv7t/uKXuoD16T7cQ/p//OgCXjx3i2qbb510jXsrt3PRR6j+/P9I4olfaZthsZ92npKLP8vFjNaQyONgna6vXFGQaueZ2yNbNnH5wHxzjM3gx+VNq7Dag2lmvzm3x7x8djsP36Olq495787YJPxMtATj/sAE4ckhP/PyU7PwNVW965+l/bxa1yWwI8qTLjnD9ds3xQ3MGDyZOPSBrJvK+M/o9mN47XUu6+JjBrt8CTUyOI54xoEcHjBpkCYU4/J4mqlZA6M83qCM11UPuKmK5NdzRJz8TIdtY9KO6tq/FUVqnpkwF3vDVtrVJ56UPS+fgZ2aZcvUxru/5OKmVcNRNb3pnUJMgZxSXSGRHR1d+w905A4bV0QxvirfzCFLH9V+8tvFEIvdceqend6Q1hs96pMyg+o7GZ6/O92hIFlC1BjcRMLiXFcET9sJ7pwR5Q0avPm5PALkTJk0j6qBL+aXV2NGYxsRLR2PcQbm2dDVY0J+hywfhEVLegUDQvd8ybj/8Y/xIt4/Ifs7quk3pjDMg83uVa5LktJ3unjlDCUJoPrXffGt/l/bq1VD1ezDNEQrS5oPa9B52+zjCntinrhtjmjoXVSsg9Acc5Cg0VYS3Qk0jgI4+zrI2NYmsM1M7MHduhXXdbu3buOzfdbUJ44LuJnbpbNZivKPwfEI8VSfu17kmE+SMwHUNwpQiItcHkRuV4R0RnrDfrr5l02203pFzTSKRM7LXhdwlx2ZHdfo1laBOuUapwD3fHZZzfdWO3vrUvMb0lKuPwaTLjnA0COasCbExJHmd99dendx1e+kYSwB7w51N/Yj3merN2bs+gYpu6tGxDQ7o29Vo4uthD0RWaSk5apIBGkREAaHvpXewqr0qgZ/OMF7+8VF4/4axruvqtEkmnPdWnevbw/tqvwe/A/Wd6lxObq9A0e/B6zthBCdZDBIQg+s74sObvoELDh8IAOhsR1KVRjxUsYBwaRABAkJ19Pra07kCwvr+9MWHO9s6ta3N6Yh/d+aBGNCjg9O56i/K7Wce6NpXjcq9ZbNSBtgaSIjjsiZJ6KP5MJ6YMAqvXXVUTkdpUvH9nGqqE9eFiv5Sen0QaqAbZnYBcjWIm0/dJ6dse9R3xNL/Pdl4vD5xqLYm4TIdJhK5L5X+0u67WxfHVKZvVwECeogjEWHE7t1x3QlDQ+9JZ2DPDtivTxfHDJTKMDq3s9pXWEI8rwZxxJCexgg7rz3dFM4dNPjMERBH7YEXLj8SB/fvZp8vl+G2LVzP66TXm9ek6icgurSrxdF71jtCSb9lXUCoDli1u1SG0bltLXbt0tbXXKWbmNQ7pQRlJsKM9gSRa0DjFQJuAeE+FzNyAkh0gsZ5yQShvlOdc23VdkqxmhxQxQIiaXAUmtilc1s8MWEU7vxO1umc66S2/g7r3w1D7TjsDm2SeOvaY/GvHxzm7Pcte8anUkH18xy/r3tkrFRqgrtzrdUERJiJiUCYev0YRxiOHNQDg3vlxomb7em5I6qhu3ZynLVtPEJBkUgQ1KNNJrIN2TSCqqtJuGy56nnccOJQnHxAb5x32MBIgkWhz8xuk0zg+cuO1MqYyOmovGUiTfNRfHdEf/zftw5w2ajVYRccPhAXjd5dE8LuB/naVUcby6kLVyV4toYkzTN1IqMG9cjZ1rtrO9foNooGccqB2bBXbzqLJBH22S1rrjNp1CZzqt7hRtUgurWvxUMXjcAutnak72UKQ3b8Q5r25effq0lmAyaU5kiOgAifsOi9B290o/4OtK3xahDBaXWCOnvvcep7c9ahyYcqFhDZz6bG8eIVRzq2+pGDerhsu9761EeXj3x/JB4ZP9IZ3RxgiE5xHE8Bq2qpBkkEdNLMVTXJ7Ag9aljoxz87Dh/e9A3f35u0cihbtumZvHTlUc41azxCyymfZmIiyqYyN5mYmIE5Nx+PS47dA0BWg/jB0Xvgnu8OQyJBkeaYKCGjr+1QkyB00eosQVbHfPuZBzoLzXjLpL7q10wmCN8+tJ9LS1IdbLs2Sfz81H3wv9/aH/26t8PAHu50HX5pQ3Sh2NkQcGBC1yAmXjrad7+OdTWY/YvjQuafZO9v0W0n4r8OzppavILKG/Fkqg5TkIM+dycnismu54tG747JPznaac+qHZk601rDgM7xD2nvgV9HXJtI4LA9eriuo4plmtF+00nueSFeoepNFBhkYrKu6d+Og8Z53mSLut+lFFStgNArXB8tX3fCUDz1w8Owd+/OGKjl59E7QdXAlbbwl3MPcX7r2bHOcSgB5saim5geumgEbjltX99yEpFrRNgmmXBeuDAtU91il/buuRRedG3g6D172WU0N41TDtwNY4b2wpXfzDqdazw+CDXqT2oCIpEgvHjFkbjr7IO1sGLLNnvkEMs5v79hMlqUaA0loPW1Hfy0gzMO6evk4fFqJ+pr2DW9ncWRQ+rx9rVjcjpKv8ll7bVV35RNeUuYgIhwXgUROQJFL6myubsdyO569pr5OGfOT+71TG28pz2TnSj3eapJhycfsCv2qO/odJ5qP6UF6CNrvYNVdan20wMI/AYUiQThL+cegld/fJTLBApYz3YPTy6us0e6fS1q33+MH4FHxo/MMSPpmotpomVQrjQA+PePDnMl6FR0buc2GWa1ptIIiFgXDGrJ6J2aPnr40TF7GPfX91EN9+GLRqCXjyNY5/ffPhBLtTz5en6Vo/esd4Vk/vTkvV0pkwlARy1ioiZJSGVsE1WR7JB6VIfSavwERMe6GtekL8A9ustxUmsax969O2Pv3p3x9bZG3PPGIudZjxrUA5N/crQxYV5UAbFy007XmslBUR6qfEmno1HmPLvcYdf0+9lzTT+HabuAnD46HwMxQVcAABR+SURBVNw4FiN+PRmAlU1VLbKjP5NJlx1hXPRINQ0iq5326NgG++7Wxd4W/kyH9e+Kj77cmHOvpmNNGsQvTrUGPXd99+CcTvvC0QMxenBPJy2GmoeTdSLbph9NWJnagUmT9j7zCw4fiL+/u9S5jp7dVp0yk2EcPrgn/jF+BM69/wMAuWZnta8azCxcvcX9e5CTmnMjp7wcMqA7BvfqiPkrN+MPZx2EK5+YCSA3WiqrQYiJKVb0UWCU+QJ6I1cvn18H4OX/DeuLq765p/M9yET0/SMH4Ren7uvSDrq51PWE07GFOaqihsLpI/chu3RCx7oaXJuHA9YVFqo5qYmyQkx/gbp3aIP5t57gmuizR31HY3nDFsgB4Jh2Lhw9MFJ5HQ3HY9pwNIiQS/rJj6juEn2woQuL1646GpN/crRjntEHH387T5+JnT1+vz5dcmYmA/oqXISj9qx3hENUHrxgBP5+4aEu86l1vlxMIZz9e7TH/RccivZtaozanJ4zSQWCeCdh6u+HqW3UGkxMXlPOzaft6xvUkPVBWN8PGdDN9zxhgwavyfWZiw/H9Seqd4hBRPjluH1x//nDfc+hblEPm/cKG9PkwDipWg3CNdPTfuhB2Tf1Bqo6vbDFW/xQxwWuBawaCQFnDe+Hn/1nLgBLmDUmrcbhlQ/6aCkqT/3wMAwf2B1/PucQ7NK5Dh3rajD3luPzWs5QN3kkEwnXs8oE+CCiEMVJ3aGuxukErvv3nND9HQ1HmTRUQIBnu2+ZfO4l39QNw/p3dc1iV3MiXv/JMTnOYn1OSxStStcg8uGf3x+Jj5dtRJf2tTjGkESvkGywYQMVlT9JvQ5OgrqQQbIyY+mht1nzEXB1QIYEax/3QEsXdN4yex/519vc65QnXQNIxsH9u2HBSreWce5hAwPLo84QNO5zIrdKpEFUrYDwOiKBaBPGjt93F7w8z8qoGXWmspcoUUjOCw5yjThrkuTb2Y4a1N0lIMLmSQDA8IHWKN47vyDf1c7049RjydjrIgD+E5iKQd5J1xwfifXdW49h9+4vIKKXYdqNY9GpbY3xXF3a17oc7IBnQBPhQt/YuxcemLrEGOkEWGtHjN17l5zthw/uicMNGoki6FlHzUzrpWOdO+xbnce0Xpge9KEPDBQq3fevTt/flQLGhGNiUgO+gOfqFRjetUz0Y70JA6NagpWlICizgWlmf5zEJiCIqB+AhwHsCiAD4D5m/iMR3QzgvwGo2UQ3MvML9jE3ABgPIA3gcmZ+Oa7y6S+m+rh/SD6cBbeegNokYfBNLwIoPIvr7j064PzDBuDcCNP7vehRTF6O2asXTtp/Vxw1pB7XPz2nWYua5zuVv4292lVC80GkM9kFbgpdh8CPvt3aOesW53vmrIC2vnvj4sPu3e9W8kkhryYxRl00Ue+gopg2Dx/c09e0AlhrRxSCXz0+d+kR6NnJ386+5H9Owu43vGD8zWtScjRsz7P54KaxOTZ5L8MHdseUq4+JtAKeGihkHG3L+q5mZT944aH48RMzsXF7U859XzZmCJjhDMhUmxkxsDt+ZadTz7fJX3v8UAzo3h7H7eM/GbSSnNQpAD9h5o+IqBOAGUT0qv3bHcx8u74zEe0D4DsA9gWwG4DXiGhPZo4lp63eCfTq1BYPXTQCw3zy2yiUvfib++yCV+evLjgfSiJBuGVcbk5+E95GVptI+F63bW0Sf/reIfjoS2vlrJ3NSAesXpbTD9oN/5m5ImRvy5m2bmuDa95GhrNLZBY7d8yb1xyLJz5chhufmeN6Ro9PGIUvtYAAE0qAqRDCy+xZyOo8YcLM7/dCZGCY+eVXp++HIb3cC+qUKg9PEN5ihw2ugu5TOXDVynAJbYCh45097sfAiKsDqiLpvrxHxo/EEHsBo2P36oW+3dph4/amHK29e4c2uPm0ffHI+18glWGnzJeMGZwzm900BqhJEPp1b49TDujt5L9q1yaJC+wkjA9cMBx9uuYKOTU4CJt5XyxiExDMvBLASvvzFiJaACBo1ZVxAB63V5hbYi89OgLAe3GUz/uSeZO7BXHX2Qdjw/bGWPOhmMIUAXcyPD+Uqcy0+Hk+zLvleNTVJPCfmSuwR33wS3fE4B74z8wVrnkQ6Ux2lF7sTi2ZIMdEoI/cRw3q4WtWUSgBVptMuEbZUcNc/apdFxxhUStROWdUrpZZqCmnGEQVovnQqW0tFtx6guMDdMJPY+4DTdfRQ9QBOJkA/G73mYtHY9KcFZi/YnPOb45PwZCrbeGvTgTg39bGDM01/wHZsNeg6LdiUhIfBBENBHAwgGkARgO4lIjOAzAdlpaxAZbweF87bDmCBUqzaE6H1bY2id5dSrOmralhhpm2lA2zuStOKVvq3FuOD73m/37rAJx/+ED06tzW6YDTGcaFowfi+dkrnElKxaTQ/iNpGDkC5pnUJnx9EJrlR58bE8Ytp+3rjFqjUE4NQt17sUtgyqzb3PVBwtDDXP249vi9cMXjH2OQzwBp/75dsH/fLjj3/mk5vwU1o0Lr8JT9e2PN5p2h/pViEbuAIKKOAP4N4Epm3kxE9wL4Jaz3+5cAfgfgIpjbXE7NEdEEABMAoH//wh9SS1DTgwgaPYWVvZ0jIIpjnfNLPPjBTWOd0Xvb2qSTr+e6E/bC5Y/NxP59u6BjXQ0W/4+/Lbw5KEGo8hlFxa8DMs2kNuH3s67J5NO+zrcTsUUlKK9P3KhnFGcyUcfEFLMKoedi8uPYob0w++bjm3WdYt5GIkH4/pH+s+SLTawCgohqYQmHR5n5aQBg5tXa738FMMn+uhxAP+3wvgByjN/MfB+A+wBg+PDhBT/6QsMuS4Uex+4lrPNRM69Pj7COcnPwswkfMqA7pl4/pqjXMq2vccoBu2HtlgZ8b2R+zn6vU1qhnrXf8732hL3wfy8t9DUt6pvj7MTzyU9V9Gs7GkT+ZZh249hI+0XNVtxc9FxMsZzffkal8RbEQ5xRTATgfgALmPn32vbetn8CAP4LwFz780QA/ySi38NyUg8B8EFc5WvpGsRZh/bDi3NX4Xv2co5/+t4wTFtsLVwfJtzqapKYf+vxOUnDWivvXHescUGkZIGjqdMP7oOFq7biCs8CQqpJ+HXAFx8zOGexF/fxhWkQrQnK+RAdv/TzXrzhp3Gh52JqLtccvxe+2rDDFeiikhhGXbGuJRKnBjEawLkA5hDRTHvbjQDOJqKDYAnWpQB+AADMPI+I/gVgPqwIqEviimACWv4LvEvntnjximw20pP2742T7MXmo4wg81lGtKXTt1t4yKIfT198ODZub3Rtq6uxkux5Uf17wZP6XKGoLbt9FQrF5IMwXSPuUH89F1NzOaBvV7zuWYjrhH13xS9P3w9nHtLXfFArIM4opndgbkfmYGjrmNsA3BZXmXRauoCIm2k3jo28ZnVrZlj/buE72YSZmMLQD6vU9hVHFJMXNWvclLyxmMStqSQShHMNUWiticoZZuZJnA28NRBV3a9GCu3c9SZV6CTKlo7jg4jx9gb06IBJlx2Rs857sSmVptKaqVoBUakmAKFwghY3ikK+s51bI+T5Gxf7xaw9AMChdpqZkbt3D9mzeqlaAdHSo5jCqO9Uh3NHDcCZw/tWhamoFDiRYwU2DZcPokI1iFKYmErFiN27Y94tx+fkThKyVO2TKWeoYDEIWiFOKIyMM3tdfBB+OIKhQm5PhEMwVTv0rNQRnlA4yldZaNOoBh+EorLvTlBUrYCoBBVZKC5qYlahbYOqYB6E0rzjzEMmtByqVr+q71SH4QO64aIjdi93UYQWQjbtc2HH67On4/ALJROE4QOih+3GgXo0FSr/BA9VKyDa1ibx1I8OL3cxhBaEyrpZqAahR8bFoUF8/uuTin7OfMmGuYqEqAaq1sQkCF6yJqbCjte1htYeJedHhfmohRBEQAiCTUPKSo9eV1tYDit9nYbWHiXnhzMPojJvT/AgAkIQbJSAKHSt8eqajyISohqophYtCIE0OhpEoQKi8jtNNZmwQhUkwYMICEGwabDX8K4rME16NWgQzmRCERBVQeW3aEGIiOODqCnstajU/Es6ajJhobPNhdZF5bdoQYiI44MoUEBUk4lJNIjqQASEINg0NlODKOda0aVCrU8+du9eZS6JUAriXHK0H4CHAewKIAPgPmb+IxF1B/AEgIGwVpT7NjNvsJco/SOAkwBsB3ABM38UV/kEwYua3Na2wDDXakgh36VdLd67YQx6dqwrd1GEEhDnkCcF4CfMvDeAUQAuIaJ9AFwPYDIzDwEw2f4OACfCWod6CIAJAO6NsWyCkMMfv3MQLh87BEN3LWyhmmpwUgNA7y7tquZeq53YapmZVyoNgJm3AFgAoA+AcQAesnd7CMDp9udxAB5mi/cBdCWi3nGVTxC89O3WHld9c8+C00hIpylUGiVp0UQ0EMDBAKYB2IWZVwKWEAGgjJl9ACzTDltubxOEVkGlZnAVqpfYk/URUUcA/wZwJTNvDhidmX7IWS2WiCbAMkGhf//+xSqmIBSFa0/YC6P36FnuYghCUYhVgyCiWljC4VFmftrevFqZjuy/a+ztywH00w7vC2CF95zMfB8zD2fm4fX19fEVXhAK4OJjBuPAfl3LXQxBKAqxCQg7Kul+AAuY+ffaTxMBnG9/Ph/As9r288hiFIBNyhQlCIIglJ44TUyjAZwLYA4RzbS33QjgfwH8i4jGA/gSwJn2by/ACnFdBCvM9cIYyyYIgiCEEJuAYOZ34J/ycaxhfwZwSVzlEQRBEPJD4vIEQRAEIyIgBEEQBCMiIARBEAQjIiAEQRAEIyIgBEEQBCPEnDNZudVARGsBfFHg4T0BrCticVoDcs/VgdxzddCcex7AzKEzjVu1gGgORDSdmYeXuxylRO65OpB7rg5Kcc9iYhIEQRCMiIAQBEEQjFSzgLiv3AUoA3LP1YHcc3UQ+z1XrQ9CEARBCKaaNQhBEAQhgKoUEER0AhEtJKJFRHR9+BGtAyLqR0RvENECIppHRFfY27sT0atE9Jn9t5u9nYjoTvs5zCaiYeW9g8IgoiQRfUxEk+zvuxPRNPt+nyCiNvb2Ovv7Ivv3geUsd6EQUVcieoqIPrHr+rAqqOMf2216LhE9RkRtK62eiegBIlpDRHO1bXnXKxGdb+//GRGdb7pWVKpOQBBREsA9AE4EsA+As4lon/KWqmikAPyEmfcGMArAJfa9XQ9gMjMPATDZ/g5Yz2CI/W8CgHtLX+SicAWsNc8VvwFwh32/GwCMt7ePB7CBmQcDuMPerzXyRwAvMfNQAAfCuveKrWMi6gPgcgDDmXk/AEkA30Hl1fPfAZzg2ZZXvRJRdwC/ADASwAgAv1BCpSCYuar+ATgMwMva9xsA3FDucsV0r88C+CaAhQB629t6A1hof/4LgLO1/Z39Wss/WCsPTgYwBsAkWCnm1wGo8dY3gJcBHGZ/rrH3o3LfQ5732xnAEm+5K7yO1Xr13e16mwTg+EqsZwADAcwttF4BnA3gL9p21375/qs6DQLZxqZYbm+rKGy1+mAA0wDswvbqfPbfXvZulfAs/gDgWgAZ+3sPABuZOWV/1+/JuV/79032/q2JQQDWAnjQNqv9jYg6oILrmJm/AnA7rAXGVsKqtxmo7HpW5FuvRa3vahQQpkWMKiqUi4g6wloL/Epm3hy0q2Fbq3kWRHQKgDXMPEPfbNiVI/zWWqgBMAzAvcx8MIBtyJodTLT6e7ZNJOMA7A5gNwAdYJlYvFRSPYfhd49FvfdqFBDLAfTTvvcFsKJMZSk6RFQLSzg8ysxP25tXE1Fv+/feANbY21v7sxgN4DQiWgrgcVhmpj8A6EpEarVE/Z6c+7V/7wLg61IWuAgsB7CcmafZ35+CJTAqtY4B4BsAljDzWmZuAvA0gMNR2fWsyLdei1rf1SggPgQwxI6AaAPL2TWxzGUqCkREAO4HsICZf6/9NBGAimY4H5ZvQm0/z46IGAVgk1JnWwPMfAMz92XmgbDq8XVm/h6ANwCcYe/mvV/1HM6w929VI0tmXgVgGRHtZW8aC2A+KrSObb4EMIqI2tttXN1zxdazRr71+jKA44iom615HWdvK4xyO2XK5Ag6CcCnAD4HcFO5y1PE+zoCljo5G8BM+99JsOyvkwF8Zv/tbu9PsCK6PgcwB1aUSNnvo8B7PwbAJPvzIAAfAFgE4EkAdfb2tvb3Rfbvg8pd7gLv9SAA0+16/g+AbpVexwBuAfAJgLkA/gGgrtLqGcBjsHwsTbA0gfGF1CuAi+x7XwTgwuaUSWZSC4IgCEaq0cQkCIIgREAEhCAIgmBEBIQgCIJgRASEIAiCYEQEhCAIgmBEBIQgRICIbrKzic4moplENJKIriSi9uUumyDEhYS5CkIIRHQYgN8DOIaZG4ioJ4A2AN6FFX++rqwFFISYEA1CEMLpDWAdMzcAgC0QzoCVF+gNInoDAIjoOCJ6j4g+IqIn7ZxYIKKlRPQbIvrA/jfY3n6mvb7BLCJ6qzy3Jgj+iAYhCCHYHf07ANoDeA3AE8z8pp0Dajgzr7O1iqcBnMjM24joOlgze2+19/srM99GROcB+DYzn0JEcwCcwMxfEVFXZt5YlhsUBB9EgxCEEJh5K4BDYC3MshbAE0R0gWe3UbAWoJpKRDNh5c0ZoP3+mPb3MPvzVAB/J6L/hrUIjiC0KGrCdxEEgZnTAKYAmGKP/L1LORKAV5n5bL9TeD8z8w+JaCSAkwHMJKKDmHl9cUsuCIUjGoQghEBEexHREG3TQQC+ALAFQCd72/sARmv+hfZEtKd2zFna3/fsffZg5mnM/HNYq57paZoFoeyIBiEI4XQEcBcRdYW17vciWOamswG8SEQrmflY2+z0GBHV2cf9FFbWYACoI6JpsAZlSsv4rS14CFamzlkluRtBiIg4qQUhZnRndrnLIgj5ICYmQRAEwYhoEIIgCIIR0SAEQRAEIyIgBEEQBCMiIARBEAQjIiAEQRAEIyIgBEEQBCMiIARBEAQj/x//MdvyTOg/vwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x2af85a368358>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(ham_eq[\"output/generic/temperature\"])\n",
    "plt.xlabel(\"Steps\")\n",
    "plt.ylabel(\"Temperature [K]\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Structure analysis\n",
    "\n",
    "We will now use the `get_neighbors()` function to determine structural properties from the computed trajectories. We take advantage of the fact that the TIP3P water model is a rigid water model which means the neighbors of each molecule never change. Therefore they need to be indexed only once"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_struct = ham_eq.get_structure(iteration_step=-1)\n",
    "\n",
    "# Get the indices based on species\n",
    "O_indices = final_struct.select_index(\"O\")\n",
    "H_indices = final_struct.select_index(\"H\")\n",
    "\n",
    "# Getting only the first two neighbors\n",
    "neighbors = final_struct.get_neighbors(num_neighbors=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Distribution of the O-H bond length\n",
    "\n",
    "Every O atom has two H atoms as immediate neighbors. The distribution of this bond length is obtained by:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEQCAYAAAC9VHPBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAEJtJREFUeJzt3XmQZWV9xvHvI0hYREGnMRYwaUyQkuBGmpQR4wLRGkElWkSlXFApJxUrxBUdiyQaK6nCipWgaKRGJLgQNCIoihsqSBKBclhEFncHHCXOIIkasYJDfvnjXpyeZmb60NP3nJ55v5+qrr5n6fv++q3ufvo9y3tSVUiS2nW/oQuQJA3LIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMbtOnQBXSxbtqymp6eHLkOSdihXX3317VU1Nd9+O0QQTE9Ps2bNmqHLkKQdSpJbuuznoSFJapxBIEmNMwgkqXEGgSQ1bmJBkOTsJOuT3LCFba9PUkmWTap9SVI3kxwRnAOsmLsyyYHA04BbJ9i2JKmjiQVBVV0O3LGFTf8IvAHwGZmStAT0eo4gybOBH1bV1/psV5K0db3dUJZkT+BU4Okd918JrARYvnz5BCvTzmp61cW/fr32tGMHrERa2vocEfw2cBDwtSRrgQOAa5L85pZ2rqrVVTVTVTNTU/PeIS1JWqDeRgRV9XVgv3uWx2EwU1W391WDJOneJnn56HnAFcAhSdYlOWlSbUmSFm5iI4KqOmGe7dOTaluS1J13FktS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1bmJBkOTsJOuT3DBr3d8n+UaS65NcmGSfSbUvSepmkiOCc4AVc9ZdAhxWVY8GvgW8aYLtS5I6mFgQVNXlwB1z1n2+qjaOF68EDphU+5KkboY8R/By4DNb25hkZZI1SdZs2LChx7IkqS2DBEGSU4GNwLlb26eqVlfVTFXNTE1N9VecJDVm174bTHIi8Ezg6KqqvtuXJG2u1yBIsgJ4I/Dkqrqzz7YlSVs2yctHzwOuAA5Jsi7JScC7gL2BS5Jcl+TMSbUvSepmYiOCqjphC6vfN6n2JEkL453FktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklq3MSCIMnZSdYnuWHWugcnuSTJt8ef951U+5KkbiY5IjgHWDFn3Srgi1V1MPDF8bIkaUATC4Kquhy4Y87q44D3j1+/H/jjSbUvSeqm73MED62q2wDGn/fruX1J0hxL9mRxkpVJ1iRZs2HDhqHLkaSdVt9B8OMkDwMYf16/tR2ranVVzVTVzNTUVG8FSlJr+g6Ci4ATx69PBD7Rc/uSpDkmefnoecAVwCFJ1iU5CTgNeFqSbwNPGy9Lkga066TeuKpO2MqmoyfVpiTpvluyJ4slSf0wCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXETm2JCWsqmV13869drTzt2wEqk4TkikKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktS4TkGQ5Mgu67pK8pokNya5Icl5SXZf6HtJkrZP1xHBGR3XzSvJ/sBfADNVdRiwC/CChbyXJGn7bXPSuSR/ADwBmEry2lmbHsjoD/j2tLtHkl8BewI/2o73kiRth/lGBLsBD2D0h3vvWR8/A45fSINV9UPg7cCtwG3AT6vq8wt5L0nS9tvmiKCqvgx8Ock5VXXLYjSYZF/gOOAg4L+BjyZ5UVV9aM5+K4GVAMuXL1+MpiVJW9D1eQS/kWQ1MD37a6rqqAW0+UfA96tqA0CSCxgdftosCKpqNbAaYGZmphbQjiSpg65B8FHgTOAs4O7tbPNW4PFJ9gR+CRwNrNnO95QkLVDXINhYVe9ZjAar6qok5wPXABuBaxn/5y9J6l/XIPhkklcCFwL/e8/KqrpjIY1W1ZuBNy/kayVJi6trEJw4/nzKrHUFPHxxy5Ek9a1TEFTVQZMuRJI0jE5BkOQlW1pfVR9Y3HIkSX3remjoiFmvd2d0pc81gEEgSTu4roeGTp69nORBwAcnUpEkqVcLnYb6TuDgxSxEkjSMrucIPsnoKiEYTTb3SOBfJ1WUJKk/Xc8RvH3W643ALVW1bgL1SJJ61unQ0HjyuW8wmnl0X+CuSRYlSepP10NDzwP+HrgMCHBGklOq6vwJ1iYtmulVFw9dgrRkdT00dCpwRFWtB0gyBXwBMAgkaQfX9aqh+90TAmM/uQ9fK0lawrqOCD6b5HPAeePl5wOfnkxJkqQ+zffM4t8BHlpVpyR5LvBERucIrgDO7aE+SdKEzXd453Tg5wBVdUFVvbaqXsNoNHD6pIuTJE3efEEwXVXXz11ZVWsYPbZSkrSDmy8Idt/Gtj0WsxBJ0jDmC4KvJnnF3JVJTgKunkxJkqQ+zXfV0KuBC5O8kE1/+GeA3YDnTLIwSVI/thkEVfVj4AlJngocNl59cVV9aeKVSZJ60fV5BJcCly5Wo0n2Ac5iFC4FvLyqrlis95ckddf1hrLF9g7gs1V1fJLdgD0HqkOSmtd7ECR5IPAk4KUAVXUXzmYqSYMZYr6ghwMbgH9Ocm2Ss5LsNUAdkiSGCYJdgcOB91TV44BfAKvm7pRkZZI1SdZs2LCh7xolqRlDBME6YF1VXTVePp9RMGymqlZX1UxVzUxNTfVaoCS1pPcgqKr/BH6Q5JDxqqOBm/quQ5I0MtRVQycD546vGPoe8LKB6pCk5g0SBFV1HaM7lCVJA/MpY5LUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaN9RcQ9KSMb3q4s2W15527Ba3zV4v7UwcEUhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklq3GBBkGSXJNcm+dRQNUiShh0RvAq4ecD2JUkMFARJDgCOBc4aon1J0iZDjQhOB94A/N9A7UuSxnqfhjrJM4H1VXV1kqdsY7+VwEqA5cuX91SddnRzp5Qe6j2kHckQI4IjgWcnWQt8GDgqyYfm7lRVq6tqpqpmpqam+q5RkprRexBU1Zuq6oCqmgZeAHypql7Udx2SpBHvI5Ckxg36qMqqugy4bMgaJKl1jggkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGtd7ECQ5MMmlSW5OcmOSV/VdgyRpk10HaHMj8LqquibJ3sDVSS6pqpsGqEWSmtf7iKCqbquqa8avfw7cDOzfdx2SpJEhRgS/lmQaeBxw1Ra2rQRWAixfvrzXurTjmF518WBtrT3t2N7aliZpsJPFSR4AfAx4dVX9bO72qlpdVTNVNTM1NdV/gZLUiEGCIMn9GYXAuVV1wRA1SJJGhrhqKMD7gJur6h/6bl+StLkhRgRHAi8Gjkpy3fjjmAHqkCQxwMniqvp3IH23K0naMu8slqTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxg06DbW0EH1OPb0ts+twSmrtyBwRSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWrcIEGQZEWSbyb5TpJVQ9QgSRrpPQiS7AK8G3gGcChwQpJD+65DkjQyxIjg94HvVNX3quou4MPAcQPUIUlimCDYH/jBrOV143WSpAEMMQ11trCu7rVTshJYOV78nyTfnGhV81sG3D5wDUuFfbHJMuD2vG3oMpYEfy42WSp98VtddhoiCNYBB85aPgD40dydqmo1sLqvouaTZE1VzQxdx1JgX2xiX2xiX2yyo/XFEIeGvgocnOSgJLsBLwAuGqAOSRIDjAiqamOSPwc+B+wCnF1VN/ZdhyRpZJBHVVbVp4FPD9H2dlgyh6mWAPtiE/tiE/tikx2qL1J1r/O0kqSGOMWEJDXOIJijy/QXSZ6X5KYkNyb5l75r7Mt8fZFkeZJLk1yb5PokxwxRZx+SnJ1kfZIbtrI9Sd457qvrkxzed4196NAPLxx//9cn+UqSx/RdY1/m64tZ+x2R5O4kx/dV231WVX6MPxidvP4u8HBgN+BrwKFz9jkYuBbYd7y839B1D9gXq4E/G78+FFg7dN0T7I8nAYcDN2xl+zHAZxjdJ/N44Kqhax6oH54w63fjGTtrP3Tpi/E+uwBfYnRO9Piha97ahyOCzXWZ/uIVwLur6r8Aqmp9zzX2pUtfFPDA8esHsYX7QXYWVXU5cMc2djkO+ECNXAnsk+Rh/VTXn/n6oaq+cs/vBnAlo/uEdkodfiYATgY+BizpvxMGwea6TH/xCOARSf4jyZVJVvRWXb+69MVbgBclWcfoP56T+yltSXLqlHs7idEoqUlJ9geeA5w5dC3zMQg212X6i10ZHR56CnACcFaSfSZc1xC69MUJwDlVdQCjQyMfTNLqz1SnqVNakeSpjILgjUPXMqDTgTdW1d1DFzKfQe4jWMK6TH+xDriyqn4FfH88B9LBjO6Y3pl06YuTgBUAVXVFkt0ZzbGypIfBE9Jp6pQWJHk0cBbwjKr6ydD1DGgG+HASGP1eHJNkY1V9fNiy7q3V/962psv0Fx8HngqQZBmjQ0Xf67XKfnTpi1uBowGSPBLYHdjQa5VLx0XAS8ZXDz0e+GlV3TZ0UX1Lshy4AHhxVX1r6HqGVFUHVdV0VU0D5wOvXIohAI4INlNbmf4iyVuBNVV10Xjb05PcBNwNnLIz/tfTsS9eB7w3yWsYHQZ5aY0vldjZJDmP0eHAZeNzIm8G7g9QVWcyOkdyDPAd4E7gZcNUOlkd+uGvgYcA/zT+T3hj7UCTr90XHfpih+GdxZLUOA8NSVLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINA6kGSM5Jck+SIoWuR5jIIpAlLshewH/CnwDMHLke6F4NAGkvyliSv346vn07yyyTXzV5fVb8AHgZcBrxzvO8eSa5Lctd4zippMAaBtLi+W1WPnb0iyUOAPYGfM5qfiqr65Xi/Jmco1dJiEKhpSU4dP5f5C8AhE2rmL4G3AzcyeqSntKQYBGpWkt9jNL3244DnAot+IjfJNKPn+H4EuBn43cVuQ9peBoFa9ofAhVV1Z1X9jFnPW0jyJ0nekeRdSf5uvO4Ts7Z/NMkuHdr4W+Ct4+m5DQItST6PQK271zzsSY4EZqrqVePlM5M8GZj9oJn7zfcIwiSPZTTSeGKSdzN6cM/XF61yaZE4IlDLLgeeM76CZ2/gWeP1JwFnzNn3cODQcSi8n24ned8GPGvWU6oegyMCLUGOCNSsqromyUeA64BbgH8bb7o/45FCkoMY3QOwAXhdVX01ybHA1LbeO8lRwF5V9cVZ7f04yV5JHlxVdyz+dyQtjE8ok+ZI8ijgVGA9o1D4K+C9wPOr6q4kfwOcX1Vfn/N108Cnquqw+9DWWkaHoW5fnOql+84gkBZJkgOBrwA/mXsvwRb23QO4gtHI4lGOEDQkg0CSGufJYklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLj/h/z0Biizlha3AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x2af85aa1fe48>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "bins = np.linspace(0.5, 1.5, 100)\n",
    "plt.hist(neighbors.distances[O_indices].flatten(), bins=bins)\n",
    "plt.xlim(0.5, 1.5)\n",
    "plt.xlabel(\"d$_{OH}$ [$\\AA$]\")\n",
    "plt.ylabel(\"Count\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Distribution of the O-O bond lengths\n",
    "\n",
    "We need to extend the analysis to go beyond nearest neighbors. We do this by controlling the cutoff distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "neighbors = final_struct.get_neighbors(cutoff=8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "neigh_indices = np.hstack(neighbors.indices[O_indices])\n",
    "neigh_distances = np.hstack(neighbors.distances[O_indices])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One is often intended in an element specific pair correlation function. To obtain for example, the O-O coordination function, we do the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Getting the neighboring Oxyhen indices\n",
    "O_neigh_indices  = np.in1d(neigh_indices, O_indices) \n",
    "O_neigh_distances = neigh_distances[O_neigh_indices]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEcCAYAAAAlVNiEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAFhNJREFUeJzt3Xu4XXV95/H3h4Q74WbSDgLxgFKmYEVsFAxqFdoRCUq1rcVqC1Um06l1UBxsLI6X1k5xxloU++hDba0KIoI4TyVeEWjHKlESiYiICIRLQa7KRSgI850/1gpstueys3PW3ifJ+/U8+znr/Nbl9z3rrLM/Z132WqkqJElbtq3GXYAkafwMA0mSYSBJMgwkSRgGkiQMA0kShoG2AEm+kOS4cdcxiCTHJ/naRsy/yfysmlsMA22U9s3riiQPJPlRkg8l2XWGebZN8ldJbkzyYJJrkpycJF3UWFUvqaqPdbHscUryziRn9rZtrj+rumcYaGhJ3gy8BzgZ2AU4FHgK8JUk20wz67nAEcBRwALg94HlwPs7LXgSSeZ3tNwk2WqmNmnOqCpfvjb4BewM3A+8sq99J+B24LVTzHcE8O/A3n3thwCPAk+bYr51wFuB7wE/Bj4KbNeO2w24ALijHXcBsFfPvJcAJ7TDxwP/CvwNcDfw7kn6mgf8GXAtcB+wen29wFLgW8A97delff38Zbv8B4GnTdG2C/D3wK3AvwHvBub11Pe1nmW+H7gJuLet4/lt+5HAw8DP2t/D2kl+1q2AtwE3tL+TjwO7tOMmgAKOA24E7gROGfd25Wt8L/9L0bCWAtsB5/c2VtX9wBeA35hivt8AVlXVTX3zrQJupgmLqbwaeDHwVOCXaN7ooHnT+yjNXslimjfdD06znEOA64BfoHmj7ncS8CqaPZedgdcCDyTZHVgJfAB4EvA+YGWSJ/XMu34vZwHNm/BkbR8DHqEJhoOB/wScMEWt3wKeCewOfBI4N8l2VfVF4H8C51TVTlV10CTzHt++XgTsSxPU/evlecD+NOv97Ul+eYo6tJkzDDSshcCdVfXIJONubcdPNd+tU4ybbj6AD1bVTVV1N82b+KsAququqvpMVT1QVfe1435tmuXcUlWnV9UjVfXgJONPAN5WVVdXY21V3QUsA66pqk+0854NfB94ac+8/1hVV7bjf9bfRvOm/hLgjVX106q6nWYv5djJCq2qM9uf75Gq+mtgW5o370G8GnhfVV3XhvRbgWP7Do29q6oerKq1wFpgslDRFqCT46XaItwJLEwyf5JA2KMdT5L7e9oPaNv3m2KZj803hd69iRuAJ7d97EDzhnokzSEjgAVJ5lXVozMsZzJ70xwi6vdkHv9vv7eOPWdYdm/bU4CtgVt7zpdvNVVN7XmZE9q+i2ZPZbrAnK7eG2j+5n+xp+1HPcMP0Ow9aAvknoGG9Q3gIeAVvY1JdqT5z/erAO0hjPWvG4ELgUOS7N0333No3oQvmqbP3nkWA7e0w2+m+W/5kKraGXjB+sVOsZyZbtV7E82hqH630LyZ91pMc9x/umX3tt1Es94WVtWu7Wvnqjqwf6Ykzwf+FHglsFtV7UpzrmL9zzXTz9Ff72Kaw1O3zTCftkCGgYZSVfcA7wJOT3Jkkq2TTNBcKXQz8Ikp5ruQJig+k+TAJPOSHAqcBXyoqq6ZptvXJ9mrPXb/Z8A5bfsCmvMEP2nHvWMjf7yPAH+RZL/2CqBntOcFPg/8UpLfSzI/ye/S7O1cMOiCq+pW4MvAXyfZOclWSZ6aZLLDWgto3rzvAOYneTvNnsF6twET01yhdDbwpiT7JNmJx88xTHZoT1s4w0BDq6r/RfOm/F6aq11W0fzne0RVPTTNrL8FXAx8keZKmDNprq55wwxdfpLmjfS69vXutv00YHuaQ0yXtsvdGO8DPt32dW9b2/bteYOjafZE7gLeAhxdVdMd2prMHwDb8PiVUefRHCLr9yWak/E/oDnE8+888XDSue3Xu5KsmWT+f6AJ5X8Brm/nn2kdawuVKh9uo7kvyTqaSyYvHHct0ubIPQNJkmEgSfIwkSQJ9wwkSRgGkiTm2CeQFy5cWBMTE+MuQ5I2GatXr76zqhZt7HLmVBhMTExw2WWXjbsMSdpkJOm/RcpQPEwkSTIMJEmGgSQJw0CShGEgScIwkCRhGEiSMAwkScyxD51JeqKJFSunHLfu1GUjrESbO/cMJEmGgSTJMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJAwDSRKGgSQJw0CShGEgScIwkCTRcRgkeVOSK5N8N8nZSbbrsj9J0nA6C4MkewL/DVhSVU8H5gHHdtWfJGl4XR8mmg9sn2Q+sANwS8f9SZKGML+rBVfVvyV5L3Aj8CDw5ar6cv90SZYDywEWL17cVTnqyMSKlVOOW3fqshFWImljdHmYaDfgGGAf4MnAjkle0z9dVZ1RVUuqasmiRYu6KkeSNI0uDxP9OnB9Vd1RVT8DzgeWdtifJGlIXYbBjcChSXZIEuAI4KoO+5MkDamzMKiqVcB5wBrgiravM7rqT5I0vM5OIANU1TuAd3TZhyRp4/kJZEmSYSBJMgwkSRgGkiQMA0kShoEkCcNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJAwDSRKGgSQJw0CShGEgScIwkCRhGEiSMAwkSRgGkiQMA0kShoEkCcNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEh2HQZJdk5yX5PtJrkry3C77kyQNZ37Hy38/8MWq+u0k2wA7dNyfJGkInYVBkp2BFwDHA1TVw8DDXfUnSRpel3sG+wJ3AB9NchCwGjixqn7aO1GS5cBygMWLF3dYjjYHEytWTtq+7tRlGzzPTKZb5jCmq2O2+5I2VJfnDOYDzwI+VFUHAz8FVvRPVFVnVNWSqlqyaNGiDsuRJE2lyzC4Gbi5qla1359HEw6SpDmmszCoqh8BNyXZv206AvheV/1JkobX9dVEbwDOaq8kug74w477kyQNodMwqKrLgSVd9iFJ2nh+AlmSZBhIkgwDSRKGgSQJw0CShGEgScIwkCRhGEiSMAwkSRgGkiQMA0kSA4ZBksMGaZMkbZoG3TM4fcA2SdImaNq7liZ5LrAUWJTkpJ5ROwPzuixMkjQ6M93Cehtgp3a6BT3t9wK/3VVRkqTRmjYMquqfgX9O8o9VdcOIapIkjdigD7fZNskZwETvPFV1eBdFSZJGa9AwOBf4MPAR4NHuypE0GyZWrJxy3LpTl42wEm0qBg2DR6rqQ51WIkkam0EvLf1ckj9OskeS3de/Oq1MkjQyg+4ZHNd+PbmnrYB9Z7ccSdI4DBQGVbVP14VIksZnoDBI8geTtVfVx2e3HEnSOAx6mOjZPcPbAUcAawDDQJI2A4MeJnpD7/dJdgE+0UlFkqSRG/YW1g8A+81mIZKk8Rn0nMHnaK4eguYGdb8MfLqroiRJozXoOYP39gw/AtxQVTd3UI8kaQwGOkzU3rDu+zR3Lt0NeLjLoiRJozXok85eCXwT+B3glcCqJN7CWpI2E4MeJjoFeHZV3Q6QZBFwIXBeV4VJkkZn0KuJtlofBK27NmBeSdIcN+iewReTfAk4u/3+d4HPd1OSJGnUZnoG8tOAX6yqk5O8AngeEOAbwFkjqE+SNAIzHeo5DbgPoKrOr6qTqupNNHsFp3VdnCRpNGYKg4mq+k5/Y1VdRvMITEnSZmCmMNhumnHbz2YhkqTxmSkMvpXkP/c3JnkdsHqQDpLMS/LtJBcMU6AkqXszXU30RuCzSV7N42/+S4BtgJcP2MeJwFXAzkNVKEnq3LRhUFW3AUuTvAh4etu8sqouGmThSfYClgF/CZy0MYVKkroz6PMMLgYuHmL5pwFvobmn0aSSLAeWAyxevHiILiSYWLFyZMtcd+qyWa+ji/qlDdHZp4iTHA3cXlXTnluoqjOqaklVLVm0aFFX5UiSptHlLSUOA16WZB3wKeDwJGd22J8kaUidhUFVvbWq9qqqCeBY4KKqek1X/UmShufN5iRJA9+obqNU1SXAJaPoS5K04dwzkCQZBpIkw0CShGEgScIwkCRhGEiSMAwkSRgGkiQMA0kShoEkCcNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEiN60pmkuWNixcoNnmfdqcs6qERziXsGkiTDQJJkGEiSMAwkSRgGkiQMA0kShoEkCcNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJLoMAyS7J3k4iRXJbkyyYld9SVJ2jjzO1z2I8Cbq2pNkgXA6iRfqarvddinJGkIne0ZVNWtVbWmHb4PuArYs6v+JEnD63LP4DFJJoCDgVWTjFsOLAdYvHjxKMrZok2sWDnluHWnLhtqvtnua64bZl1szjbX3/OWpvMTyEl2Aj4DvLGq7u0fX1VnVNWSqlqyaNGirsuRJE2i0zBIsjVNEJxVVed32ZckaXhdXk0U4O+Bq6rqfV31I0naeF3uGRwG/D5weJLL29dRHfYnSRpSZyeQq+prQLpaviRp9vgJZEmSYSBJMgwkSRgGkiQMA0kShoEkCcNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJNHhk8606ZlYsXLcJWgDjPL3tSlsG9PVuO7UZSOsZNPknoEkyTCQJBkGkiQMA0kShoEkCcNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJDoOgyRHJrk6yQ+TrOiyL0nS8DoLgyTzgL8FXgIcALwqyQFd9SdJGl6XewbPAX5YVddV1cPAp4BjOuxPkjSk+R0ue0/gpp7vbwYO6Z8oyXJgefvtQ0m+22FNs2EhcOe4ixjAnK4z73lscE7X2cM6h9Dze+430jqnqWMmc2p9TmH/2VhIl2GQSdrq5xqqzgDOAEhyWVUt6bCmjbYp1AjWOdusc3ZZ5+xJctlsLKfLw0Q3A3v3fL8XcEuH/UmShtRlGHwL2C/JPkm2AY4F/qnD/iRJQ+rsMFFVPZLkT4AvAfOAf6iqK2eY7Yyu6plFm0KNYJ2zzTpnl3XOnlmpMVU/dxhfkrSF8RPIkiTDQJI0gjBIsneSi5NcleTKJCdOMk2SfKC9bcV3kjyrZ9xxSa5pX8eNuc5Xt/V9J8nXkxzUM25dkiuSXD5bl3ptRJ0vTHJPW8vlSd7eM24ktwgZsM6Te2r8bpJHk+zejhvV+twuyTeTrG3rfNck02yb5Jx2na1KMtEz7q1t+9VJXjzGGk9K8r122/xqkqf0jHu0Zz13dhHHgHUen+SOnnpO6Bk3qr/1Qer8m54af5DkJz3jRrI+e/qbl+TbSS6YZNzsbZtV1ekL2AN4Vju8APgBcEDfNEcBX6D5bMKhwKq2fXfguvbrbu3wbmOsc+n6/mlus7GqZ9w6YOEcWZ8vBC6YZN55wLXAvsA2wNr+eUdZZ9/0LwUuGsP6DLBTO7w1sAo4tG+aPwY+3A4fC5zTDh/QrsNtgX3adTtvTDW+CNihHf6v62tsv7+/6/W4AXUeD3xwknlH+bc+Y51907+B5gKYka7Pnv5OAj45xd/0rG2bne8ZVNWtVbWmHb4PuIrm08m9jgE+Xo1LgV2T7AG8GPhKVd1dVT8GvgIcOa46q+rrbR0Al9J8dmKkBlyfUxnZLUKGqPNVwNld1DKddpu7v/126/bVf1XFMcDH2uHzgCOSpG3/VFU9VFXXAz+kWccjr7GqLq6qB9pvx7VtDrIupzLKv/UNrXMs2yZAkr2AZcBHpphk1rbNkZ4zaHdhDqZJ4l6T3bpiz2naOzVNnb1eR7M3s14BX06yOs0tNjo3Q53PbXeDv5DkwLZtTq7PJDvQ/OF/pqd5ZOuz3Q2/HLid5g1pyu2zqh4B7gGexAjX5wA19urfNrdLclmSS5P8Zhf1bWCdv9UezjovyfoPpo502xx0fbaH2/YBLuppHtn6BE4D3gL8vynGz9q2ObIwSLITzR/7G6vq3v7Rk8xS07R3ZoY610/zIpo/uD/taT6sqp5Fc/jo9UleMMY61wBPqaqDgNOB/7N+tkkWNfb1SXOI6F+r6u6etpGtz6p6tKqeSfPf9HOSPL1vkrFvnwPUCECS1wBLgP/d07y4mlsq/B5wWpKndlHjgHV+DpioqmcAF/L4f7Uj3TYHXZ80h17Oq6pHe9pGsj6THA3cXlWrp5tskrahts2RhEGSrWneEM6qqvMnmWSqW1eM9JYWA9RJkmfQ7LIdU1V3rW+vqlvar7cDn6WDwwWD1llV967fDa6qzwNbJ1nIHFyfrWPp2w0f5frs6fMnwCX8/OGJx9ZbkvnALsDdjOGWK9PUSJJfB04BXlZVD/XMs35dXtfOe3CXNU5XZ1Xd1VPb3wG/2g6P5fY1063P1nTbZtfr8zDgZUnW0RzSPTzJmX3TzN62uaEnMzb0RZNQHwdOm2aaZTzxBPI36/GTStfTnFDarR3efYx1LqY59ra0r31HYEHP8NeBI8dY53/g8Q8UPge4sZ1vPs2JuX14/ATygeOqs51u/ca745jW5yJg13Z4e+D/Akf3TfN6nniS7tPt8IE88STddXRzAnmQGg+mOUm4X1/7bsC27fBC4Bq6u2hgkDr36Bl+OXBpOzzKv/UZ62zH7U9zIUPGsT77ankhk59AnrVts9MfoC3qeTS7J98BLm9fRwF/BPxRO01oHoRzLXAFsKRn/tfSvAH/EPjDMdf5EeDHPeMva9v3bVf8WuBK4JQx1/knbR1raU4mLu2Z/yiaK3uuHXed7XTH05zo6p13lOvzGcC32zq/C7y9bf9zmv+wAbYDzm23wW8C+/bMf0q7Lq8GXjLGGi8EbutZ1//Uti9t/6bWtl9fN+Z1+Vc92+bFwH/smX9Uf+sz1tl+/07g1L55R7Y++/p9IW0YdLVtejsKSZKfQJYkGQaSJAwDSRKGgSQJw0CShGEgScIwkGZVktOTrEny7HHXIm0Iw0CaJUl2BH4B+C/A0WMuR9oghoG2OEnemeS/b8T8E0kebO96+Ziq+inNcxwuAT7QTrt9+xCUh9v7Q0lzkmEgDefaau56+ZgkTwJ2AO4DHgWoqgfb6Tq/6Zq0MQwDbRGSnNI+/u9CmhuQdeFtwHtp7r1zQEd9SJ0wDLTZS/KrNHd0PBh4BTDrJ3fbB/gsBc6hearbgdNNL80188ddgDQCzwc+W+1jIXsfYp7kd2jusDoPuKeqTpmsbYA+3g38eVVVEsNAmxzDQFuKn7s9b5LDaG6XfmL7/YeT/NokbftX1dVTLTjJM2n2OJ6X5G9pbit8RRc/hNQVDxNpS/AvwMvbK3sW0DxmE5pHl57eN+07J2l7eIblvwd4aVVNVNUEcBDuGWgT456BNntVtSbJOTQPfbmB5slWAFvT7jEk2YfmMwI397dV1fVTLTvJ4TRPaftqT3+3Jdkxye71xOc6S3OWD7fRFivJr9A8Dep2mmD4HzSfE3hCW1Xd2TffBM1Tp6Z6iPpkfa2jOfx050zTSuNgGEgbKMneNM9lvqv/swaTTLs98A2a5+7+insKmqsMA0mSJ5AlSYaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJOD/A/2g10fS+d02AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x2af85aad3748>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "bins = np.linspace(1, 5, 100)\n",
    "count = plt.hist(O_neigh_distances, bins=bins)\n",
    "plt.xlim(2, 4)\n",
    "plt.title(\"O-O pair correlation\")\n",
    "plt.xlabel(\"d$_{OO}$ [$\\AA$]\")\n",
    "plt.ylabel(\"Count\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "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.5.4"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}