{
 "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 import Project\n",
    "import ase.units as units\n",
    "import pandas"
   ]
  },
  {
   "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": "89aa5565a8384b5693970a7a5664fba9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "707bbf57e7174bf68263179508d40234",
       "version_major": 2,
       "version_minor": 0
      },
      "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, pbc=True)\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": [],
   "source": [
    "water_potential = pandas.DataFrame({    \n",
    "    'Name': ['H2O_tip3p'],\n",
    "    'Filename': [[]],\n",
    "    'Model': [\"TIP3P\"],\n",
    "    'Species': [['H', 'O']],\n",
    "    'Config': [['# @potential_species H_O ### species in potential\\n', '# W.L. Jorgensen et.al., The Journal of Chemical Physics 79, 926 (1983); https://doi.org/10.1063/1.445869\\n', '#\\n', '\\n', 'units real\\n', 'dimension 3\\n', 'atom_style full\\n', '\\n', '# create groups ###\\n', 'group O type 2\\n', 'group H type 1\\n', '\\n', '## set charges - beside manually ###\\n', 'set group O charge -0.830\\n', 'set group H charge 0.415\\n', '\\n', '### TIP3P Potential Parameters ###\\n', 'pair_style lj/cut/coul/long 10.0\\n', 'pair_coeff * * 0.0 0.0 \\n', 'pair_coeff 2 2 0.102 3.188 \\n', 'bond_style harmonic\\n', 'bond_coeff 1 450 0.9572\\n', 'angle_style harmonic\\n', 'angle_coeff 1 55 104.52\\n', 'kspace_style pppm 1.0e-5\\n', '\\n']]\n",
    "})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/srv/conda/envs/notebook/lib/python3.7/site-packages/pyiron/lammps/base.py:170: UserWarning: WARNING: Non-'metal' units are not fully supported. Your calculation should run OK, but results may not be saved in pyiron units.\n",
      "  \"WARNING: Non-'metal' units are not fully supported. Your calculation should run OK, but \"\n"
     ]
    }
   ],
   "source": [
    "job_name = \"water_slow\"\n",
    "ham = pr.create_job(\"Lammps\", job_name)\n",
    "ham.structure = water\n",
    "ham.potential = water_potential"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The job water_slow was saved and received the ID: 1\n"
     ]
    }
   ],
   "source": [
    "ham.calc_md(temperature=300, \n",
    "            n_ionic_steps=1e4, \n",
    "            time_step=0.01)\n",
    "ham.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4eb876d7df744ab4a8e47da346cb02b2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "NGLWidget()"
      ]
     },
     "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": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The job water_fast was saved and received the ID: 2\n"
     ]
    }
   ],
   "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 = water_potential\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": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cd96be9a5a654dfcb69546d8a77d0e0a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "NGLWidget(max_frame=1000)"
      ]
     },
     "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": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEHCAYAAABSjBpvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA5y0lEQVR4nO3dd5xU5fX48c/Zzi69SVk6iAWlbbB3E8CGxoaJGluILdGf0aghidGEfG0xtkTFJLbYuxEwggI2lCK9L70sbSkLC2w9vz/undk7de/WWXbP+/Wa187cuXfmubO7c+7TziOqijHGGONHUqILYIwx5tBhQcMYY4xvFjSMMcb4ZkHDGGOMbxY0jDHG+GZBwxhjjG8piS5AVYnICOAJIBn4p6o+GG//9u3ba8+ePeujaMYY02jMmTNnh6p2CN9+SAUNEUkG/g78ENgIzBKRj1R1SaxjevbsyezZs+uriMYY0yiIyLpo2w+15qlhQK6qrlbVYuANYFSCy2SMMU3GoRY0ugIbPI83uttCiMgYEZktIrO3b99eb4UzxpjG7lALGhJlW0QeFFUdr6o5qprToUNEk5wxxphqOtSCxkagm+dxNrA5QWUxxpgm51ALGrOAfiLSS0TSgNHARwkukzHGNBmH1OgpVS0VkVuB/+EMuf23qi5OcLGMMabJOKSCBoCqTgQmJrocxhjTFB1qzVP15sWv1/DRfOsuMcYYLwsaMbz63XomLcxLdDGMMaZBsaARQ0pyEqXltqqhMcZ4WdCIISVJKLOgYYwxISxoxJCcJFbTMMaYMBY0YnBqGuWJLoYxxjQoFjRiSE4SSsuspmGMMV4WNGJISbY+DWOMCWdBI4bkJBs9ZYwx4SxoxGCjp4wxJpIFjRhs9JQxxkSyoBGDjZ4yxphIFjRisJqGMcZEsqARg/VpGGNMJAsaMSQnJdk8DWOMCWNBIwaraRhjTCQLGjEkJ1ufhjHGhLOgEYONnjLGmEgWNGKw0VPGGBPJgkYM1qdhjDGRLGjEYKOnjDEmkgWNGJIEFAsaxhjjZUEjBhGw1iljjAllQSOGJBFULWoYY4yXBY0YBKtpGGNMuAYXNETkjyKySUTmubdzPM/dKyK5IrJcRIbXcTkArLZhjDEeKYkuQAx/U9VHvRtE5ChgNHA00AWYIiKHq2pZXRTAjRmoVtw3xpimrsHVNOIYBbyhqkWqugbIBYbV1ZslBWoadfUGxhhzCGqoQeNWEVkgIv8WkTbutq7ABs8+G91tdSJQuSi35iljjAlKSNAQkSkisijKbRTwDNAHGATkAX8NHBblpaJ+o4vIGBGZLSKzt2/fXq0yJiUF+jSqdbgxxjRKCenTUNWz/ewnIs8DH7sPNwLdPE9nA5tjvP54YDxATk5Ojb72raZhjDEVGlzzlIh09jy8CFjk3v8IGC0i6SLSC+gHzKyrciRZ77cxxkRoiKOnHhaRQThNT2uBXwCo6mIReQtYApQCt9TVyCmoGDFlNQ1jjKnQ4IKGql4V57lxwLj6KEeSZ8itMcYYR4NrnmooxO13t5qGMcZUsKARQ3ByX2KLYYwxDYoFjRiCaURsxVdjjAmyoBFDsE/D6hrGGBNkQSOGwJBby3RrjDEVLGjEYENujTEmkgWNGCpSoye4IMYY04BY0IghMB/c1tMwxpgKFjRisNToxhgTyYJGDNanYYwxkSxoxGBpRIwxJpIFjRgsjYgxxkSyoBGDWE3DGGMiWNCIwYbcGmNMJAsaMVgaEWOMiWRBI4aK0VOJLYcxxjQkMRdhEpGCSo4VIE9VD6/dIjUMwXka1j5ljDFB8VbuW6Wqg+MdLCJza7k8DY7VNIwxpkK85qmLfRzvZ59DUqCmYXPCjTGmQrygcYeInBTvYFVdXcvlaTCsT8MYYyLFCxorgUdFZK2IPCQig+qpTA1Ckg25NcaYCDGDhqo+oaonAKcBO4EXRGSpiPxBRBpl57dXoHHKZoQbY0yFSofcquo6VX3I7RT/CXARsLTOS5ZgNrnPGGMiVRo0RCRVRM4XkVeBScAKGnEHeIBluTXGmEjx5mn8ELgCOBeYCbwBjFHVwnoqW0JVjJ4yxhgTEK+m8VtgBnCkqp6vqq/WVsAQkUtFZLGIlItITthz94pIrogsF5Hhnu1DRWSh+9yTInX7rZ5kNQ1jjIkQryP8DFV9XlV3isjJInItgIh0EJFeNXzfRcCPgS+8G0XkKGA0cDQwAviHiCS7Tz8DjAH6ubcRNSxDXDbk1hhjIvnp07gPuBu4192UCvynJm+qqktVdXmUp0YBb6hqkaquAXKBYSLSGWipqjPUyevxMnBhTcpQGbE0IsYYE8FPwsKLgAuAQgBV3Qy0qKPydAU2eB5vdLd1de+Hb68zFUNu6/JdjDHm0BIv91RAsaqqiCiAiGT5eWERmQJ0ivLUWFX9MNZhUbZpnO2x3nsMTlMW3bt3r6Sk0VkaEWOMieQnaLwlIs8BrUXk58B1wPOVHaSqZ1ejPBuBbp7H2cBmd3t2lO2x3ns8MB4gJyenWt/61qdhjDGR/EzuexR4B3gX6A/8QVWfqqPyfASMFpF0t7O9HzBTVfOAvSJyvDtq6mogVm2lVlgaEWOMieSnpoGqTgYm19abishFwFNAB2CCiMxT1eGqulhE3gKWAKXALapa5h52E/Ai0AxnkuGk2ipP1DK6P23IrTHGVIg3ue9jVT0v3sF+9olGVd8H3o/x3DhgXJTts4EBVX2v6rI0IsYYEyleTeNkEfkozvMCHFXL5WkwAn0aNuTWGGMqxAsao3wcX1xbBWlogn0aCS6HMcY0JDGDhqpOr8+CNDSWsNAYYyL5mdzXJAVqGnsOlCS4JMYY03BY0IjhqM4tAbj1tbmUlpUnuDTGGNMw+Mk9dZ6INLng0iwtOXh/wsK8BJbEGGMaDj/BYDSwUkQeFpEj67pADdHeg6WJLoIxxjQIfmaEXwkMBlbhrBM+Q0TGiEhdJS1scA6WlFW+kzHGNAG+mp1UtQAnjcgbQGeczLffi8gv67BsCdcmMxWAnYWNdmSxMcZUiZ8+jfNF5H3gc5y1NIap6khgIHBnHZcvoSbfcRpgzVPGGBPgJ/fUpcDfVDVklT1V3S8i19VNsRqG9s3T6dEuk4KDNuzWGGPAR9BQ1avjPPdZ7Ran4WmZkUqBzdUwxhjAR9AQkb1EZtPYA8wGfq2qq+uiYA1Fy2Yp7LagYYwxgL+O8MeAu3CWV83G6cd4HqdT/N91V7SGoUe7LOau383363cluijGGJNwfoLGCFV9TlX3qmqBuyreOar6JtCmjsuXcBcOcpYi/yZ3R4JLYowxiecnaJSLyGUikuTeLvM81+iz+eX0cOLio5+uoMzWfjXGNHF+gsZPgauAbcBW9/6VItIMuLUOy9YgJCVJ8P6aHfsSWBJjjEm8uEFDRJKBm1T1fFVtr6od3Pu5qnpAVb+qp3I2CGc/9kXlOxljTCMWN2i463MPraeyGGOMaeD8TO6b6y77+jZQGNioqu/VWakamEuGZvPOnI2JLoYxxiScnz6NtkA+cCZwvns7ry4L1dCcdniHRBfBGGMaBD8zwq+tj4IYY4xp+PwkLDxcRD4TkUXu42NF5Hd1X7SGw7tOeLkNuzXGNGF+mqeeB+4FSgBUdQHOwkxNhidmUFJuS78aY5ouP0EjU1Vnhm1rUrnCvTWN4lILGsaYpstP0NghIn1wZ3+LyCVAk1o0+0dHdwret6BhjGnK/ASNW4DngCNEZBNwO3BTTd5URC4VkcUiUi4iOZ7tPUXkgIjMc2/Pep4bKiILRSRXRJ4UEYn+6rWveXoKf7noGACKyyxoGGOaLj+jp1YDZ4tIFpCkqntr4X0XAT/GCUbhVqnqoCjbnwHGAN8CE4ERwKRaKIsvaSlOfC0ptY5wY0zT5Wc9jXTgYqAnkBK4wFfVB6r7pqq61H1tX/uLSGegparOcB+/DFxIAoJGcVlZfb2lMcY0OH6apz4ERuF0fhd6bnWll4jMFZHpInKKu60r4J2SvdHdFpWIjBGR2SIye/v27bVSqLRk56Mqsj4NY0wT5ieNSLaqjqjqC4vIFKBTlKfGquqHMQ7LA7qrar6IDAU+EJGjgWhVkpjtRO6aH+MBcnJyaqU9KT1Q07CgYYxpwvwEjW9E5BhVXViVF1bVs6taGFUtAorc+3NEZBVwOE7NItuzazawuaqvXxOpyRY0jDHGT/PUycAcEVkuIgvcEUwL6qIwItLBTceOiPQG+gGrVTUP2Csix7ujpq7GaTarN8GO8DLrCDfGNF1+ahoja/tNReQi4CmgAzBBROap6nDgVOABESkFyoAbVXWne9hNwItAM5wO8HrrBAfrCDfGGPA35HadiJwM9FPVF0SkA9C8Jm+qqu8D70fZ/i7wboxjZgMDavK+NZFmzVPGGOMrYeF9wN04+acAUoH/1GWhGqJATcNGTxljmjI/fRoXARfgDrNV1c1Ai7osVEMUqGlYn4YxpinzEzSKVVWpyD2VVbdFapjSbMitMcb4ChpvichzQGsR+TkwBSddepNSETSsI9wY03RVGjRU9VHgHZwO6v7AH1T1qbouWEMTCBp//O+SBJfEGGMSx8+QW1R1MjC5jsvSoKUm11tSXWOMabD8NE8ZKjrCAUosPboxpomyoOGTNyPv/01clsCSGGNM4ljQqIIBXVsC8P36XQkuiTHGJEbMPg0RWUj0TLICqKoeW2elaqDevelE+v/uE4Z0b5PoohhjTELE6wg/r95KcYhIT0mmU8sMCotKE10UY4xJiJhBQ1XX1WdBDhVZ6cm8OXsD955zBK0z0xJdHGOMqVd+ck8dLyKzRGSfiBSLSJmIFNRH4RqiVdudRQvHvr8owSUxxpj656cj/GngCmAlTlryG3DSmjdpu/YXJ7oIxhhT73yNnlLVXCBZVctU9QXgjLotVsPVLDUZgHK1xIXGmKbHz4zw/SKSBswTkYdx1vFukkkLAZqlJXOgpIxyixnGmCbIT03jKiAZuBUnPXo34OK6LFRDluRO8lOraRhjmiBfK/e5dw8A99dtcRq+wMRwq2kYY5qieJP73lLVy2JN8muKk/vAmdkI1qdhjGma4tU0bnN/2iQ/j0BNw5IWGmOaoph9Gqqa5969WVXXeW/AzfVTvIZH3LrG7v0lCS6JMcbUPz8d4T+Msm1kbRfkUPHMlUMASEmy9TWMMU1PzKAhIje5/Rn9RWSB57YGWFB/RWxYBndvw1XH92D3AatpGGOannh9Gq8Bk4D/A+7xbN+rqjvrtFQNXMcW6ezeX8KB4jKapSUnujjGGFNv4vVp7FHVtap6BbARKMEZRdVcRLrXVwEbou7tMgHYsGt/gktijDH1y0/CwluBrThrhE9wbx/X5E1F5BERWeY2d70vIq09z90rIrkislxEhnu2DxWRhe5zT4p3Kb161qOdMyF+Xb4FDWNM0+KnI/x2oL+qHq2qx7i3ms7RmAwMcF9nBXAvgIgcBYwGjgZGAP8QkUD7zzPAGKCfextRwzJUW4+2Tk1jXX5hoopgjDEJ4SdobAD21OabquqnqhpYyehbINu9Pwp4Q1WLVHUNkAsME5HOQEtVnaFO/o6XgQtrs0xV0TozlRbpKWzcdSBRRTDGmITwk7BwNTBNRCYARYGNqvpYLZXhOuBN935XnCASsNHdVuLeD98elYiMwamV0L177Xe/iAjtmqexs9DSoxtjmhY/QWO9e0tzb76IyBSgU5Snxqrqh+4+Y4FS4NXAYVH21zjbo1LV8cB4gJycnDrJ99EqM83W1DDGNDl+EhbeDyAiWarquxFfVc+O97yI/AwnRclZWpEydiNOFt2AbGCzuz07yvaEaZOZajUNY0yT42f01AkisgRY6j4eKCL/qMmbisgI4G7gAlX1DkH6CBgtIuki0gunw3umm9Jkr7v0rABXAx/WpAw11SYzjfx9FjSMMU2Ln47wx4HhQD6Aqs4HTq3h+z4NtAAmi8g8EXnWfe3FwFvAEuAT4BZVLXOPuQn4J07n+CqciYcJ06dDFpt2H2CX1TaMMU2Inz4NVHVD2LSIslj7+ny9vnGeGweMi7J9NjCgJu9bm07s2x4+XcG732/khlN6J7o4xhhTL/wEjQ0iciKg7rKvv8JtqmrKhnRvY8NujTFNjp/mqRuBW3CGuG4EBtGEU6N72bBbY0xT46em0V9Vf+rdICInAV/XTZEOHW2y0thScDDRxTDGmHrjp6bxlM9tTU63NpnMWruT3TZfwxjTRMRbI/wE4ESgg4jc4XmqJWD5wIGRAzrx0fzNDHpgMqv/cg5JtjCTMaaRi1fTSAOa4wSWFp5bAXBJ3Ret4evmJi4EKDhoizIZYxq/mDUNVZ0OTBeRF1V1nYi0cDbrvvorXsPmDRq79pfQOtN3lhVjjDkk+enTaCEic4FFwGIRmSMiDWa+RCK1apYavG95qIwxTYGfoDEeuENVe6hqD+DX7jYDvDnmeADrDDfGNAl+gkaWqk4NPFDVaUBWnZXoENO1TTMAthUUVbKnMcYc+nytpyEivwdecR9fCaypuyIdWg5rmUGSwObdNjPcGNP4+alpXAd0AN5zb+2Ba+uyUIeS1OQkOrbIYPMem+RnjGn84s3TyMBJIdIXWAj8WlVtXGkUXVpnWE3DGNMkxKtpvATk4ASMkcAj9VKiQ1Dn1s3Is5qGMaYJiNencZSqHgMgIv8CZtZPkQ49XVplMGXJVlSVsBTyxhjTqMSraQSbolS1tB7Kcsjq2CKDotJyCg7ax2SMadzi1TQGikiBe1+AZu5jwZkZ3rLOS3eI6NgyHYDtew+GTPirjsKiUspVaZFRs9cxxpi6ELOmoarJqtrSvbVQ1RTPfQsYHn06NAfglRnrUFUA1ufv52BJ1RY4fGv2Bo6+738c88dPg9u+X7+LwiKrwRhjGgY/Q25NJQZ0bcWwnm15acY6bn1tLgdLyjj1kanc9sbcKr3Ob95ZELz/yaI8bnxlDj/+xzfc8da8Wi6xMcZUj681wk3l2mY5yQonLMwjPdWJxdOWb/d9fFFpaK3kxv98H7y/eHNB+O7GGJMQVtOoJfeec0Tw/qeLtwKQluzv412fv58FG/fEfL4q65DP27CbV2asBWD5lr30vGcCs9fu9H28McbEYzWNWtKjXUU6rn1uH0RKsr/ht6c+MrXSfXbsK6J98/RK97vw784qvFed0JMvVzo1nU8WbSGnZ1tfZTHGmHisplGHUnzUNDbs3B/y+OIh2VH32+tjOO/+4tB9yt1OeZs6YoypLRY06lBZuVa6zykPh9YyLhjUJep+ew5UnsHlixU7gvdLysqD923CoTF14+UZa8n585SQbdv3Nu6M1xY06tCx2a2qtP9Fg7tyQu92XDwkm2evHBrynJ/1Olo2q2ht7Dd2Eg9OWgY4E2sADpaU8fiUFRGd7saY6vnDh4vZsa8o+D81cWEePxg3hVm11I+450AJny3dyoOTlkW0JCRKQoKGiDwiIstEZIGIvC8ird3tPUXkgIjMc2/Peo4ZKiILRSRXRJ6UQ+Dy+cuVO5i3YXfM5wNzOgLuHN6ftJQk/nrZQIb1Cu2D2La3iAPFZcxYlR/z9cJrNoGH5aps3n2Af365msenrOTVb9dX7USMMXEFmo9fn+n8b9W0trFy617y9xXx85dmc/1Ls3l2+ipmr91V43LWhkTVNCYDA1T1WGAFcK/nuVWqOsi93ejZ/gwwBujn3kbUW2l9uvNHh9OhRUVndVm5Bjumo/nblJXB+49eOpCurZsFH2ekVvxqOrfK4O3ZG/jDh4u44vlvmbp8W0TAef6L1dz/3yVR3+f5L9dw4oOf8+inKwA4UMVJh8aYSN7af86fp5C7bR/5+5xtKUk1u6b94d++4IKnv2amp8ayafeBiP/7REhI0FDVTz35rL4Fovf+ukSkM9BSVWeo86m9DFxYt6WsulvP7MdfLjrG9/5Tl20L3g9vyspISQ7ev3BwV2at3cXbczYCcO0Ls/ho/ubg86rKuIlLyd22r7pFN8ZU0UvfrAt5/PGCzSzJc+ZULdi4hwPF1bs4Cxy3KWy5hXvfW8gvXpmDqlJertz34SJyt+2t1nvUREPo07gOmOR53EtE5orIdBE5xd3WFdjo2Wejuy0qERkjIrNFZPb27f4n2NWGaFcC5TE6xL21icMPaxHyXJJ7pdL/sBYc1iJyqO2SvAI27tpPz3smcPZj06tURm95Nu0+wOrtFmyMqarv1oQ2FT/5WUXLwdNTc7nnvQXhh8R1sKSMyUu2snlPaLAYd9GA4P1Pl2zlw3mbWZNfyEsz1nHzq9+Hv0ydq7N5GiIyBegU5amxqvqhu89YoBR41X0uD+iuqvkiMhT4QESOpqIv1ytmPU1VxwPjAXJycuq1PtervTNfo3f7LFbvKASc9s5Wmalc9+IserTL5L7zjwYgIzU55usATPzVKXRt3YwpS7dGPDd5yVaem74agFXbC6tURu8HctKDnwOw9sFzq/QaxjRV36zawber8tlZWEx2m2bk7yvmQEkZbbPS2LGvoskqXv9jNC9+s5YHJy3jtMM7hGz/ybDujH1/UfDxnz5eQn6h8z5FpeXUtzqraajq2ao6IMotEDB+BpwH/NRtckJVi1Q1370/B1gFHI5Ts/A2YWUDm2mA+h3WgtxxI7lwcEVFaKfb9vn5sm288PVaAP4+NZcvVzpDZGONsjqqS0taZaayP0ofxOoqBgqvxyavqPax8Uxbvo3vVlftH8WYQ81Pnv+OJz/PZdmWvbTJTGPCr04GYH9Yc1TzdOea/J05G+l5z4RKO8cDjRTTV1S0jhzVuSUiEpzYe1LfdsGAARVzsepTokZPjQDuBi5Q1f2e7R1EJNm93xunw3u1quYBe0XkeHfU1NXAhwkoui8pyUmkp1R8tDsLQ4fLXvvCTB753/Lg41euPy7u653arz0A/To25/qTe/kqQ9+Ozf0Wt1bsLCzmmhdmcfn4bxk2bgpbC2wlw0PRwZIydhVWPry7qdmxr4jlW/YGsz0E7CwsJjPNCQ77i8s4rGVFU/LqHYU8NnkFd749H4D1O/ezNK8g5vD5rPTIlod2zZ2cdj8/xfm/H3F0aOPN7v2R87cKi0p5ZcZatu2tm//BRKUReRpIBya7I2e/dUdKnQo8ICKlQBlwo6oGhg/cBLwINMPpA5kU/qINyXG92wXv3/ra98GEhgBTwxIZVrYGR492WcHmo7Jy5V9frQl5vk1mKrvcP55fnNabkQM6071tJkP+NDm4zx/PP4o/ekZXXfWv74IjPcBJmJieEr+5LB7vVdS2vUV8OG8TY07tU63X2rBzP9ltmtmkxErMXruTldv2ccWw7gDs2V9CYXEpXTyj8Krq4me+YfHmAmb/7mxfaWtiufudBUxalMeUO06jY8uMar9OTSzNK6B3h6wa/V0H/OhvX7CzsJi0lNDrbBFo5mlm3loQWpvw9nMkJwkjn/iSIzq14JPbT2XF1r3069icAyVlrN5eSFFJZFNTS3ddnTGn9ua8gV3ISEni9x8uDj6/92ApOwuLQ75fHpu8gn99tYavc/N59qqhEa9ZU4kaPdVXVbuFD61V1XdV9WhVHaiqQ1T1v55jZrvNW31U9dZAk1ZDNahba/7x0yEA5O05WGuZapPDhvL1aJdJv44VnejXnNiTQd1a0zYrjdd/fnzF9pNCayhfrtwRHOkBsK2gZuPKw4fxBoLY4s172BPlaiiWeRt2c8rDU3ltps0leWPmeqYsiezPCrjk2Rnc+97C4ONfvTGXEx/8nL0H/X/e4QJ/p+GznKvqzdkbKDhYyrC/fJaQYaLb9h5k5BNf8vsPKvoC8vcVVass479YFWwtKA7rQygpKycjzd/XaKAGt2zLXo6573/86G9f8OasDdz59nzOe+qriA5wINhiISJ0bd2Mds3TQ1oxwJnTAU4aoUWb9vDmrA0A5NbRAJeGMHqq0To1rEMrmsHdW1f79b++50xeveE4endwOt+vGNadzq0qrjJP6NMu1qERAv0r1RU+W1VwmjrOffIrrnnR//Lyy9xANj/OpMi6UlauCR0H/+3qfL5fXzGB6573FnLDy7N9Hx9oC18fls/Mr/Cml7dmb+DbKvRRFZeWc80LMyMmtD7/5epqlacmAsNWv3E7o9flFzL0z1P4t9unWBV/mbgs5nOlZRqSzfrHQ7pyTNfofZRr8yv6Ife6n/WEhXlMXLgFgK+i/A9mpEXWkpY+MILMtGR6d8giNVn4n5tV+9npqznvqa/YV1TKYS3TyaujeR0WNOpQVpRfuNfMsWfxxpjj4+4TT9fWzchukxl83KdDVtT9zjqiIwCXDI09HWbzbv/p18EJCIEOvplrdkasUlhcWs4M9wtn7vrd9LxnAt/kVh6YAq8TPrJs1tqdnPfUl1VeDdGvgyVl9PntRJ78LLdOXt+P0eO/5cf/+KbKx6kqizZVpNav6qidwGu8M3tDyLbfvLOA0eO/9XX8uvxC1uwoZNry7cE2/IDPlm6LcVTdCeReC4wuWrHVuer+2sffYEB5ubJ8S+Q8iFRP9uq7hvdHRPjZCT0AuPG0Pvz3lycHL+S8ok2+9V6srYwyz6q9p9kpIClJmH/fj/jvrSdzbHZrluTtobxcQ5rCzju2C8Vl5cHgVJssaNQhEQn5A/O6YGAXOrbIqHZ7a3abihrFqEHOSK0z3ODgteLPIxl/dQ4Aj1xyLBdGSYiYmixs8XRcl5SV8+z0VXy7Oj/q1Q/AfR8uDn45vD93EweKnX/Ot288AXD+WcObVv4Z1hcTzUH3nzwQNCYsyOPqf8/kzrfns2hTAevyq3cVXZlAh+LL7lokifT3qbnc+pr/8fc3/mcO5z31VfDxnycsrVJ+sZ/9eya97p1IgY9MytF8s2oHpz0yLZhCI/zqtqZdU/n7ivjf4i1MX7Hd95Vz4O+xyL3I2Ffk/H4DI5ri2VVYzJ79Jbzy7TqGP/5FxPNJnhMa7fYn/fbcI3ljzPHB+VYd3blVnVtVvT9n+l2nc9fw/gAM7tEm6j6pyUlkpafQo20m6/P3M/aDimbKji3S+c2I/iz/08hgn0htsvU06tiL1w6jpKyca16YBcARnVrw23OO5Pje/puOws0ce1ZwxAY4zVCx5ll4O+5EhF+c1ocP5oWOVu7bsQXb3I7sqcu2ce2Ls0Kej/baizZXXNm2zkwNNk91aplB51YZFJWWRaRz/3zZNj5ZlMeIAZ1jnlugJhFot70l7Mtz+ONf8NmvTwuuy15TZeXKc1+s4sQ+zgi1+u57Ly0r52BpeciXmXdknR+B5gmvTxZtCV5MVCbQrBU+ys9bxnhp/pe4/SAvfrMWIGIm9Pr8/Xy1cgdtslLp27G5rwulj+Zvpl1WGif1bc9lz80IzkX62+UDuWhwNsf+8X+MHtadu0ccQZJAablSVq5kpCaTt+cA+YXO33PgImRfkVOmrEqCRlFpGYPdASSX53SLuk9SlD+S9JTkkP/pp64YwpSlW7locFeufWFWsNbtR492Wdx8eh/O6N+Ro7q0jLtvt7aZvDd3E6/PdGqJ5x7bmUcuObZWOv9jsaBRx07q63wZPXrpQMpVuSzGH2JVdGxR/dEoR3ZuyczfnsWwv3wW3FZUUsYXK7YzdXlkwIjFmxwxNUmCV6kZqcmkpyTx9pyNtMuKHH3z8ox1lQSNcvdnGQ99Er0tefHmgloLGp8t3crDnyxnQNc8oH7SyN/9zgKy2zTjl2f14653FvD+3E3c8cPDfR8/bfk2xk1YGnef/H3FlJdrMLNALN6O3cCXfrjxX67m5tP7xnyNkrLQq//New5GPL7yX99VvN5VQ/nR0dHm/Vb41etzAeeCxTt5dcueIkrKyik4WMr4L1Yz/ovV/O7cI/lo/mYWbNzD7N+dzQn/9zmd3BFbgZrJQTeQNYszoXb73iJ+MK5iAMDivNDVNH937pFMXJjHmFN7hyzHHE2HFunBUW2PXHosJz9UsQRCTo82zF4XP/mgiFQaMAC6t80MefzE5YN8reNTExY06km8/oT61iKsyhqYuX7tC7EDhqqyv7gseKXmDRpPfp5Ly4wU2mWl0TozlT0HSlB1xraf3r8DowZ14f+96TRlFVQysidQ03j+y9hNWSu3OuPl/TQ1VCYw6iswWTLwFTtt+Tb6dmwe0mdUFeXlylmPTeem0/pw2Q9CLxTedPsOFmzaw2S3CS/ehMvwcx37/qKIvEThHvh4CQ98vKTSmf7hw7ejefiT5Tw+eSWTbj8larCu6nLCr3y7Lm7QCF+YzCs1WdgSFpRe+Hpt8PNYluf0QXibW/fsL6HY7eOIFkOLS8tJS0nijbARe4s2FQSHs/dun8X1J/fihlN6V7lfrX3zdNKSkyguK+fY7Fb846dDQi7awJnE958bjqOwin0Q3duF/n3WdcAA69Nokrw5rz779WlkVtJhf7CkjHETlnLMH//Hlj0HWbl1b0SnXcHBUk49vAOpyUkhM2O7tm7GRYMrAmZpWfw2aT9t8U99nstVnivXqtq8+wC7CotZvHkPt70xD6iYzSviXHFe88IsTn5oqq+FtKJZuGkPa3YU8oePFsXcZ3Kc4bRe731fkXZt7PsLowaMd9y+pPsvOLpK5Yw1AOL6k3ux6P7h9HbT4hSXlXPWX6ezyh3G+fepuYxw2/vnb9wdPO7HQyqaxG46Pfo8nS9X7uCbVbE7pL0Lk01cmBfyXGpyEos3h9YAvJ/HlWF/FyVlysAHPuW/86MnkPhw3iYO/90k1uwopH2UHG8n93NGQA7u3iZYCw0f8lqZjNRkVowbyZIHhvP2jSfQsWUGD19ybMg+ZeVK26w0urWt2kVKD8/+U+44tUrHVpfVNJogbxNMnw7NyUhNjkiB4LXb7RQsV1i6pYBVMbLpdnM75735cDqE/SNWlvbAz7K24IzIWp+/n4y0pCo1163evo8z/xo7wePWgiLe9XxJ/23yCu50OyXj2V9cyubdB+nbsTm52/Yyyk2J3ylsYtvlz83wXdaAiQvzuPqEnqgqr34XOX8lp0cbcnq2ZfH9w8lKT+HxKSuC82QC9hWV8v7cTVx5XHfKypVSt/2/MGyo9GmHd+C2s/sxpLvTARuohQYs3LiHPh2aB/tdFm3aw459xVx/ci8uy+nGnHW7eO/7TQDcdlY/vludz/frd0eU+RevzGHhH4eHbNuxryhi1Fx4Qr77PlpMdSxzR0EVl4XOs3jXLesZj07jqM6RzUGXDM3myM4tQpqVRYR7Rh7BcWFr3lTG2w95WU43jurckq9yd/DgpGX06Rh95GNlAv9fA7u1pq9nvlZdsqBhePm6YSGjb8IdKCkjMy2ZotJy1mwvjGjDDugTJXVJ+KziQED5y8SlHN2lZUhn7Vcrd/DxgtAry3hOfcS5IvWbbPFgSVncgBEQWPEQ4O05Gxjasw3XvjCLb+89i04xRsOc88SXrM3fT+64kQx//Mvg9g27nLHyIsKKrXv5bo3/ppwurTLYW1QabE6MlpwuKy2Z+0c5tYtA06F3v0DTy7gJS3l95no6t8zghW+c2cLRPreycg0GjGhe/W4dt785L/g48HdzXK+29O/UIti0NObU3mSkJsdMqBd+cfDQJ8t4ZtqqYMLPulJUUs7+4lLKypV1+fv5wpPnKTDZ9c4fHR5ce6Zr62YRCQTBGVpbUwO6tmJA11Yc2bklQ2OMkqqMiDD1ztNp3zxyaG5dsaDRRPVolxm8shoQYzLS6f07MG35dg4UlwUDxQMfO2PNM9OSmXTbKdzx1nzmuJ165xwT2cEdXtPYvb+Ef0zLZfwXzoQvb9AIb1rwapaaXOPFo3bsiz3rvVvbZvTt0DwixcvWgqJgX8/c9bsYGeUcAda6Q4E37T4Q0qQVSPvSvW1mlbMRPz56MI9NXh6cUf/OnI0R+3x+5+kcFlab8dYaDxSXkZaSFMx3FD5ZMPwCIHyyab+OzUOaImfFWD3u2GznuLOO7MhDFx/DBQOd3+sTowfx5codvD17Y0gGgnDPTFsFwJod1U/E6UdRaRmnPDQ1JOmfV2qycOuZ/YJBo0vruk+BEi0oVUVdB9pwFjSaqOl3nRH3+SM6teCaE3sybfl2xn6wMGK2cLvmafRol8W7N53InHW76NEuk1S3E+6V64dx1b+cWeCBmsbjlw/i5Rlr+X79bh7+JHJIaWV9B9N/czp5uw8Gm328Hpy0jLuG949IsRIu2pDS7DbNeOfGE+nUKoODJWV8v24XP/ln9OAV3rQRzWmPTIvY9md3pFPXKuaEGtarLa2bpfHJ4i3MWJXP7z4I7R9Z+sAImkXpj3rmp0O4yW3W2V9SSitSmbRoS8R+udv2MW35Nto3T+Oru89k5dZ9HNk5tInj/VtOYtbanXEHSQC0yXJqQyLC5T/oHtzet2ML+nZswa7CYpbkFXDesZ2DtclYQ3mTk4Q2mWlxg3w8bbPSYg4fPlhSHjNgQEVt7Yph3flo3qaQJiXjsI5wAzjDCQMW3T+cD245KTg8ca7bJu1NiuYdBz60R5uQZqhT+nUIVpcDk5wuHNyVCwZGTiwMjBYpOBB7VNXNp/ehY4sMBnZrHfX5Z6evos9vJ1baXxD4srjq+B7Bbbed1S/Y5JSRmhw39Up44CwqLeP5L1YHZx+HuyusL8TbYXvJ0Oy4KWTm3/cjANq3cD7HK54PnZl9RKcWUQMGwMhjOvPE6EFOmQ+WxpwQd/Zj09l7sJQd+4rJSE3mmOxWEV/izdNTOKN/R54YPShuR2taJaN2bj/7cL78zRk8/ZMh/OLU3gCMmxh92HBZuXLl8d0jMrp6vXvTCVG392yXyen9nSv38OSCAAs8nfYB3swNgQufv1w0gEX3D4/Y11jQMK4bTukdvN88PYWM1OSIL6Wzj6yYcZ5cyXyG3u7QzHaettZoi07d+J85/OHDRYx4InLmbbA8Gf6u9irrL8jb7QzD7Na24oo//ErSO0ggvOISyAq8s7CY/3y7jqc/z2XcxKW8PCN02U9w2rzPOzb2fJRzj+nMi9cMC365hwtkPo6Vx6iyxXfaZDqf+3tzNwU7pmti1KCu9O3YglvPiD5fo7L5LUlJEhwZFPgyf+Hrtdz86hzWRmmSOrJzS57+yeDg434dm7Ny3Eh+f95R/O3ygQzt0TZqR3RJmTKgi/OZPXdlZIbX8DkkAHN+/0Peu/lEgGDGBBGxLMsxWN3LxBQ+EerG0/rw1mynXT01Jf4/1HNXDmVJXkHIl/Jp/cNWJDuuO699tz4iWWJ2m2Zs3FVxVV7V+Rjr8gt57/tN3H52v+A//oQFefz2fSfVQjfP3It4w41TkpJCmqTy3eaScROWhoyw+tPHkTmFSsvKY6a893ZAjxrUNTjsN5rBUTqlT+nXnrtHHBHzGKiYVBroK6gt5xzTmaen5nL4Yc2D+Zyq6poTe/LU506Or4kLtzB/w56Iffp0aE5KchJjzzmScROXkt2mGanJSSHrybx03TC+XZ3P4s0FPDttFXuLSikuK+eaE3sysFsrhvaoCCpXDOsWnDUdLiM1mSHd2zDj3jM5rAYTZ5sKq2mYmMIvtHp3aM70u07nhN7tePzyQXGPbeOmgPDyZuDt0yGLkQOiNz+ETz705s95YvSgStd5uO7FWTzx2Ury3KvKgyVl/PL1iqGb3hFQ6amx/wWeuXJIyOOXZqzj+S9WhwQMrxRP1SQrPYXWmWk8+ONjuO/8o+KW12vCr07m4YsrxvD3DZtMd/PpfXjl+uNiDl4IqKx/J6BZajIf3HKS7/IluR9XNaevANCueXrIOQaa7bzzOnq4k9Z+fmpv/nrpQP562aCI18lITeb0/h255Yy+zPrd2QBcOjSbpCQJBoy/XjqQc47pxP/9+FieumJwxGvcdla/4P3OrZpVOoPeWE3DeLx384kh62p08nzJBzpxe7TL4vUaZOYN+Otlg2Imc+zbsTlLPSNtvKODRg3qyn/nb2ZKnMype9z+kUDQW5pXEPIl523yiTZR68JBXfhk8RbOOvIwbjurH094sofGaocHp4a0Nn8/VwzrFvwCDCS0u2hwVwY9MDnqce/edAKTFm7h/IFdOLpLK47uUlG+pCRh2p2nc/qj0wD/81jiSUtO4vYf9iMrLYXLcrrF7BuJpme7LHq3z+K+84/iha/X8vmy6mWwHRalaenSodnMWbeLmWt2BvsWAC72kU0hIzWZZX8aEdG3cvHQ7ODx7aIMS/1/VUjfYhwWNExQ+Pj85ukprH3wXBZu3EPbWhoH3jIjhYKDpTRPT45IHvfNPWeSmZZMcVk5BQdKSE1OYsrSrcGrzoDzB3aJGTQ+XbyFHW7fQ0mpUl6u3PXOAgCuPaknHVtkkJKcxKq/nMOMVflR5yQ8PrriirQqQy77d2rB2vz9HNO1dUT/TevM2J/f0B5tQ5pSwvVol8nRXVqyeHNBpQn3/BjUvXXcXFLxZKQm8/mdpwNO09mA+/5Xrdfp2T4rYhh1z3ZZvHrDcdWehR+tz8wrPE+TqR4LGqZSx2THbwqpirtHHsHY9xfRoXkGzdKSnXHxZ/TjkpzskGVKX7puGPuLS1m+ZW/E8qWjBjkL3YgIZz82nbJy5afHdefV79Yz5pU5wf2Ky8pYkldArjvP4IZTegdrTMlJwsn9QpvPounRLv4Y+LtHHMEPerahV/us4CzpWP2nr95wXLVW1RMRPv7lybw1e0Nw/oMfpx7eIWTyWkBtfXk2T09h/FVDKa3ml/w395zJ6h2FzN+wm/zCIpKShCSESr77qy3wuz++d1sevXQgKUnWOl8d0sBXTa2xnJwcnT3b/+pnpu55s68GZktXV96eA2zfW8QHczfz769Dk+99cMtJvPzNWt6b64weWvanEZVejUYr6/NfrublGeui5nz66u4zgkkNX/tuPb99fyFvjjk+ZI34RNlWcDAkMd7p/Ttw6dBunHFEhyY7/2Db3oM0T09psudfFSIyR1VzwrfbJ2fqnbezsabDGju3akbnVs14f27ksNKnP89lytKKpIBVDRjglPUXp/Xh+N7tuPTZGRSXlZOSJMGra++XzxXDujGsV5t6ywFUmY5hM8VfvHZYgkrScNRkWQHjsPqZaRSizVvYUlBRM7h3ZPwhqpUZ2K01r1zvfOke4Zk17R2yKyINJmAYU1espmEahR8eeRivhWWA9Y7A+UUtJJgb0qMNd/zwcH48pCuFRWV8v35XtWovifCnUVVLmW5MLFbTMI3CGUd05IVrfxCyLZD+5LUbjquV90hNTuJXZ/Uju00m/Tu1CK7Mdii46oSeiS6CaSSspmEajfCJcAEn9q18lFRjlZGaFDI50piasqBhGo1oTUXRJpE1JYHEh8bUloQ0T4nIn0RkgYjME5FPRaSL57l7RSRXRJaLyHDP9qEistB97kmxbGImjHdm84ijO/HslUP5188iRgw2KekpySEZiY2pqUT1aTyiqseq6iDgY+APACJyFDAaOBoYAfxDRAJ/8c8AY4B+7m1EfRfaNGwZnpQgz141lBEDOkXksTLG1ExCgoaqepfwygICMwxHAW+oapGqrgFygWEi0hloqaoz1JmN+DJwYX2W2TR80Rb0McbUroT1aYjIOOBqYA8QWEauK+BdbWaju63EvR++PdZrj8GpldC9+6EzwsXU3H3nH8UPejbtfgxj6lKdXZqJyBQRWRTlNgpAVceqajfgVeDWwGFRXkrjbI9KVcerao6q5nToULP1d82h5dqTelWaNtwYU311VtNQ1bN97voaMAG4D6cG0c3zXDaw2d2eHWW7McaYepSo0VP9PA8vAJa59z8CRotIuoj0wunwnqmqecBeETneHTV1NfBhvRbaGGNMwvo0HhSR/kA5sA64EUBVF4vIW8ASoBS4RVUDCfdvAl4EmgGT3Jsxxph6ZKnRjTHGRIiVGt3GKBpjjPHNgoYxxhjfLGgYY4zxzYKGMcYY3xp9R7iIbMcZoVUd7YEdtVicQ4Gdc9PQ1M65qZ0v1Pyce6hqxOzoRh80akJEZkcbPdCY2Tk3DU3tnJva+ULdnbM1TxljjPHNgoYxxhjfLGjENz7RBUgAO+emoamdc1M7X6ijc7Y+DWOMMb5ZTcMYY4xvFjSMMcb4ZkEjChEZISLLRSRXRO5JdHlqi4h0E5GpIrJURBaLyG3u9rYiMllEVro/23iOudf9HJaLyPDElb5mRCRZROaKyMfu40Z9ziLSWkTeEZFl7u/7hMZ8ziLy/9y/6UUi8rqIZDTG8xWRf4vINhFZ5NlW5fMUkaEistB97kl3yQl/VNVunhuQDKwCegNpwHzgqESXq5bOrTMwxL3fAlgBHAU8DNzjbr8HeMi9f5R7/ulAL/dzSU70eVTz3O/AWfDrY/dxoz5n4CXgBvd+GtC6sZ4zztLPa4Bm7uO3gGsa4/kCpwJDgEWebVU+T2AmcALOqqiTgJF+y2A1jUjDgFxVXa2qxcAbwKgEl6lWqGqeqn7v3t8LLMX5hxuF8yWD+/NC9/4o4A1VLVLVNUAuzudzSBGRbOBc4J+ezY32nEWkJc6Xy78AVLVYVXfTiM8ZZ22gZiKSAmTirOzZ6M5XVb8AdoZtrtJ5ikhnoKWqzlAngrzsOaZSFjQidQU2eB5vdLc1KiLSExgMfAccps7qiLg/O7q7NZbP4nHgNziLfgU05nPuDWwHXnCb5P4pIlk00nNW1U3Ao8B6IA/Yo6qf0kjPN4qqnmdX9374dl8saESK1rbXqMYli0hz4F3gdlUtiLdrlG2H1GchIucB21R1jt9Domw7pM4Z56p7CPCMqg4GCnGaLWI5pM/ZbcMfhdME0wXIEpEr4x0SZdshc75VEOs8a3T+FjQibQS6eR5n41R1GwURScUJGK+q6nvu5q1ulRX35zZ3e2P4LE4CLhCRtThNjWeKyH9o3Oe8Edioqt+5j9/BCSKN9ZzPBtao6nZVLQHeA06k8Z5vuKqe50b3fvh2XyxoRJoF9BORXiKSBowGPkpwmWqFO0LiX8BSVX3M89RHwM/c+z8DPvRsHy0i6SLSC+iH04F2yFDVe1U1W1V74vwuP1fVK2nc57wF2CAi/d1NZwFLaLznvB44XkQy3b/xs3D66xrr+Yar0nm6TVh7ReR49/O62nNM5RI9GqAh3oBzcEYWrQLGJro8tXheJ+NUQxcA89zbOUA74DNgpfuzreeYse7nsJwqjLBoiDfgdCpGTzXqcwYGAbPd3/UHQJvGfM7A/cAyYBHwCs6IoUZ3vsDrOP02JTg1huurc55AjvtZrQKexs0O4udmaUSMMcb4Zs1TxhhjfLOgYYwxxjcLGsYYY3yzoGGMMcY3CxrGGGN8s6BhTC0RkbFuptUFIjJPRI4TkdtFJDPRZTOmttiQW2NqgYicADwGnK6qRSLSHie77DdAjqruSGgBjaklVtMwpnZ0BnaoahGAGyQuwcmFNFVEpgKIyI9EZIaIfC8ib7t5wBCRtSLykIjMdG993e2XumtEzBeRLxJzasZUsJqGMbXA/fL/Cict9xTgTVWd7ua8ylHVHW7t4z2cmbmFInI3kK6qD7j7Pa+q40TkauAyVT1PRBYCI1R1k4i0VifFuTEJYzUNY2qBqu4DhgJjcNKSvyki14TtdjzOwjhfi8g8nDxBPTzPv+75eYJ7/2vgRRH5Oc4CYcYkVEqiC2BMY6GqZcA0YJpbQ/hZ2C4CTFbVK2K9RPh9Vb1RRI7DWURqnogMUtX82i25Mf5ZTcOYWiAi/UWkn2fTIGAdsBdnaV2Ab4GTPP0VmSJyuOeYyz0/Z7j79FHV71T1D8AOQlNdG1PvrKZhTO1oDjwlIq2BUpylNccAVwCTRCRPVc9wm6xeF5F097jf4WRUBkgXke9wLuYCtZFH3GAkOBlM59fHyRgTi3WEG9MAeDvME10WY+Kx5iljjDG+WU3DGGOMb1bTMMYY45sFDWOMMb5Z0DDGGOObBQ1jjDG+WdAwxhjj2/8HPBIlbLBvaBQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEGCAYAAACUzrmNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABBy0lEQVR4nO2dd5wU5fnAv89eo/cqRdpZALFdUFGxIIoVYyyYGEnUoIkaNSYG1MTEhGhMfiYxxkSiiUSNhFgxwYIoYiwgoDSRjnCCVKkHV5/fHzOzO7s3W+5u9/a4e76fz31m9t0p78ztzPM+9RVVxTAMwzASEcp2BwzDMIyGjwkLwzAMIykmLAzDMIykmLAwDMMwkmLCwjAMw0hKbrY7kCk6deqkffr0yXY3DMMwDirmz5+/TVU7x7Y3WmHRp08f5s2bl+1uGIZhHFSIyGdB7WaGMgzDMJJiwsIwDMNISsaEhYj8TUS2iMgSX1sHEZkhIivdZXvfdxNEZJWILBeRc3ztx4vIYve7h0REMtVnwzAMI5hMahZPAKNi2sYDM1W1EJjpfkZEBgJjgEHuPo+ISI67z5+BcUCh+xd7TMMwDCPDZExYqOpsYEdM82hgsrs+GbjY1z5FVUtVdS2wChgqIt2BNqr6vjpFrP7h28cwDMOoJ+rbZ9FVVTcBuMsubnsPYINvu2K3rYe7HtseiIiME5F5IjJv69atae24YRhGU6ahOLiD/BCaoD0QVZ2kqkWqWtS5c7UwYcMwDKOW1Lew2OyalnCXW9z2YqCXb7uewEa3vWdAe8aoqlKmfriB8sqqTJ7GMAzjoKK+hcU0YKy7PhZ4ydc+RkQKRKQvjiN7rmuq2iMiJ7pRUFf79skIz3/0OXc8t4hJs9dk8jSGYRgHFRnL4BaRZ4DTgU4iUgzcA9wPTBWRa4H1wGUAqrpURKYCnwAVwI2qWuke6rs4kVXNgVfcv4yx50A5AFv3lGbyNIZhGAcVGRMWqnplnK9GxNl+IjAxoH0eMDiNXUtITshxk1RUmRnKMAzDo6E4uBsMITfnz1wWhmEYEUxYxJDrahZVVTY3uWEYhocJixgiZigTFoZhGB4mLGLIzfHMUGaHMgzD8DBhEYPnszDNwjAMI4IJixhyQ84tqVITFoZhGB4mLGII+ywqTVgYhmF4mLCIwRMWlWaGMgzDCGPCIoYc945UmhnKMAwjjAmLGCJJeSYsDMMwPExYxOCJCBMWhmEYEUxYxOLKCAudNQzDiGDCIgZ1pYVpFoZhGBFMWMTg+bVNWBiGYUQwYRGDCQvDMIzqmLCIwRMR5rMwDMOIkBVhISK3iMgSEVkqIre6bR1EZIaIrHSX7X3bTxCRVSKyXETOyWTf1FUtrES5YRhGhHoXFiIyGPgOMBQ4GrhARAqB8cBMVS0EZrqfEZGBwBhgEDAKeEREcjLVP09GWFKeYRhGhGxoFkcCH6hqiapWAG8DXwVGA5PdbSYDF7vro4EpqlqqqmuBVTiCJkNYNJRhGEYs2RAWS4DhItJRRFoA5wG9gK6qugnAXXZxt+8BbPDtX+y2VUNExonIPBGZt3Xr1lp1TsN5FjafhWEYhke9CwtVXQb8GpgBvAosBCoS7CJBh4lz7EmqWqSqRZ07d65d/9xlpVWdNQzDCJMVB7eqPq6qx6nqcGAHsBLYLCLdAdzlFnfzYhzNw6MnsDFzfXOW5rMwDMOIkK1oqC7usjdwCfAMMA0Y624yFnjJXZ8GjBGRAhHpCxQCczPVN8vgNgzDqE5uls77nIh0BMqBG1X1SxG5H5gqItcC64HLAFR1qYhMBT7BMVfdqKqVmeqYWm0owzCMamRFWKjqqQFt24ERcbafCEzMdL/Aqs4ahmEEYRncMVhSnmEYRnVMWMRgZijDMIzqmLCIwRzchmEY1TFhEYNpFoZhGNUxYRGDpVcYhmFUx4RFDCYrDMMwqmPCIgY11cIwDKMaJixiMFFhGIZRHRMWsZi0MAzDqIYJixiqzAxlGIZRDRMWMZioMAzDqI4Jixj8ioWV/DAMw3AwYRGD+nQLS8wzDMNwMGERg1+zsKlVDcMwHExYxODXJcptalXDMAzAhEV1fKqFFRM0DMNwyNa0qreJyFIRWSIiz4hIMxHpICIzRGSlu2zv236CiKwSkeUick4m++YXDxWVZoYyDMOALAgLEekBfB8oUtXBQA4wBhgPzFTVQmCm+xkRGeh+PwgYBTwiIjmZ6p8/AqrcNAvDMAwge2aoXKC5iOQCLYCNwGhgsvv9ZOBid300MEVVS1V1LbAKGJqpjplmYRiGUZ16Fxaq+jnwW2A9sAnYpaqvA11VdZO7zSagi7tLD2CD7xDFbls1RGSciMwTkXlbt26tZf8i6xY6axiG4ZANM1R7HG2hL3AI0FJErkq0S0Bb4FtcVSepapGqFnXu3LlW/YvWLExYGIZhQHbMUGcBa1V1q6qWA88Dw4DNItIdwF1ucbcvBnr59u+JY7bKCP4S5eVmhjIMwwCyIyzWAyeKSAsREWAEsAyYBox1txkLvOSuTwPGiEiBiPQFCoG59dHRkrLK+jiNYRhGgye3vk+oqnNE5FlgAVABfARMAloBU0XkWhyBcpm7/VIRmQp84m5/o6pm7C3u91nsK6vI1GkMwzAOKupdWACo6j3APTHNpThaRtD2E4GJme4XRNeGKik1zcIwDAMsg7saplkYhmFUx4RFDP5o2ZJSExaGYRhgwqIafjPUPnNwG4ZhACYsqhFlhjLNwjAMAzBhEZfmeTkWOmsYhuFiwiIGLymvZUGuaRaGYRguJixi8MxQrQpMszAMw/AwYRGD57Jo1SzXQmcNwzBc4iblichDKey/W1XvTmN/sk5EszAzlGEYhkeiDO7RwE+T7D8eaFTCosrzWeTn8sXuA1nujWEYRsMgkbD4napOTvA9/qlPGwueGSo/N2Qlyg3DMFwS+SyejfeFiFwIoKq/T3eHso6rWeTlhCizEuWGYRhAYmExU0T6xDaKyDXA7zPVoWyjgIgrLCpMWBiGYUBiYXEbMENECr0GEZngtp+W6Y5lC1Vnar78XLHJjwzDMFzi+ixUdbqIlAKviMjFwHXAV4DhqvplPfWv3lEUESEvJ2TCwjAMwyVhnoWqzgS+BcwC+gEjGrOggIhm4QgLc3AbhmFAAmEhIntEZDfwCtAGZ2KiLb72WiEih4vIx76/3SJyq4h0EJEZIrLSXbb37TNBRFaJyHIROae2504Fz2eRn2sObsMwDI+4wkJVW6tqG3eZr6otfZ/b1PaEqrpcVY9R1WOA44ES4AWcnI2ZqloIzHQ/IyIDgTHAIGAU8IiI5NT2/Mmo0ogZqqyiKlwryjAMoymT7XIfI4DVqvoZThKgl9cxGbjYXR8NTFHVUlVdC6wChmaqQ1VVSo4I+TkCQEWVCQvDMIxEZqgFyXZOZZskjAGecde7quomAHfZxW3vAWzw7VPstmWEyirICTmaBWBObsMwDBJncB8pIosSfC9A29qeWETygYuACck2DWgLHO6LyDhgHEDv3r1r1a8qVUJungVAeYVCfq0OZRiG0WhIJCyOSGH/utTwPhdYoKqb3c+bRaS7qm4Ske7AFre9GOjl268nsDHogKo6CZgEUFRUVCv7UZWqo1nkOsLCnNyGYRiJ8yw+y/C5ryRiggKYBowF7neXL/na/ykiDwKHAIXA3Ex1qrJKCfl8FmaGMgzDSKxZZAwRaQGMBK73Nd8PTBWRa4H1wGUAqrpURKYCnwAVwI2qmrFZiapUCZnPwjAMI4qsCAtVLQE6xrRtx4mOCtp+IjCxHrpGpRcN5ZmhrD6UYRhGaqGzInKoiJzlrjcXkdaZ7Vb2iI2GMp+FYRhGCsJCRL6DU678UbepJ/BiBvuUVRwzFOSHzVCWZ2EYhpGKZnEjcDKwG0BVVxLJgWh0VKljhjKfhWEYRoRUhEWpqpZ5H0Qklzh5Do0BLxoqz4uGMp+FYRhGSsLibRG5E2guIiOBfwMvZ7Zb2SMcDWV5FoZhGGFSERY/BrYCi3FCXacDd2eyU9kkHA2VY9FQhmEYHglDZ0UkBCxS1cHAX+unS9mlsoqYPItGa3EzDMNImWSTH1UBC0WkdoWWDkJUlZwQ4TwLc3AbhmGklpTXHVgqInOBfV6jql6UsV5lkUqNdnCbz8IwDCM1YfHzjPeiARGpDeVoFqXmszAMw0guLFT17froSEPBqzrbulkeAHsOlGe5R4ZhGNknqbAQkT1E8irygTxgX12mVm3IeNFQzfJC5OeE2L2/IttdMgzDyDqpaBZRdaBE5GIyOK1ptqlSCIVARGjTPJdd+02zMAzDqPEc3Kr6InBm+rvSMKhyfRYAbZrnsduEhWEYRkpmqEt8H0NAEY253IfrswBoXZDL3lIzQxmGYaQSDXWhb70CWAeMzkhvGgB+zSI3J0RFlUVDGYZhpCIsHlPVd/0NInIykTmyGxV+zSI3JJbBbRiGQWo+iz+m2JYyItJORJ4VkU9FZJmInCQiHURkhoisdJftfdtPEJFVIrJcRM6py7mTUVVFWLPIywlRYUl5hmEY8TULETkJGAZ0FpEf+L5qA+TU8bx/AF5V1UtFJB9oAdwJzFTV+0VkPDAe+LGIDATGAIOAQ4A3ROSwTM3DXaWKq1iQExIqq0yzMAzDSKRZ5AOtcARKa9/fbuDS2p5QRNoAw4HHAVS1TFV34vhBJrubTQYudtdHA1NUtVRV1wKryGDobmVVxAyVl2NmKMMwDEigWbiZ22+LyBOq+lkaz9kPp+T530XkaGA+cAvQVVU3uefeJCLebHw9gA98+xe7bdUQkXHAOIDevWtX+/CqEw+lQ8t8AHJD5uA2DMOA1BzcJSLyGxwzUDOvUVVrm2uRCxwH3Kyqc0TkDzgmp3hIQFvgcF9VJwGTAIqKimqlEowd1ifS0RyhwjQLwzCMlBzcTwOfAn1xigquAz6swzmLgWJVneN+fhZHeGwWke4A7nKLb/tevv17AhvrcP6UycsJUW6ahWEYRkrCoqOqPg6Uq+rbqnoNcGJtT6iqXwAbRORwt2kE8AkwDRjrto0FXnLXpwFjRKRARPoChcDc2p6/JuSEhErTLAzDMFIyQ3n1LjaJyPk4o/qedTzvzcDTbiTUGuDbOIJrqohcC6wHLgNQ1aUiMhVHoFQAN2YqEiqWvByh3KKhDMMwUhIWvxSRtsDtOPkVbYDb6nJSVf0Yp2xILCPibD8RmFiXc9aG3JDlWRiGYUDyObhzgEJV/Q+wCzijXnrVQKioUr4sKWfL7gN0adMs+Q6GYRiNlGRzcFcCjXL61FR4eaHjR3/ozZVZ7olhGEZ2ScUM9Z6IPAz8i+g5uBdkrFcNBK/irDdrnmEYRlMlFWExzF3e62tTGvGcFrG0a27CwjCMpk0qM+U1KT+Fn1MLO/HOym20NWFhGEYTJ2mehYh0FZHHReQV9/NAN7y10fPLiwcDkSq0hmEYTZVUkvKeAF7DqfgKsAK4NUP9aVAU5DrFdSvVci0Mw2japCIsOqnqVKAKQFUrgHpJiss2IffuVFhinmEYTZxUhMU+EemIW7xPRE7Eyblo9OS60qLKhIVhGE2cVKKhfoBTn6m/iLwLdKYO81kcTOS4vgrTLAzDaOqkEg21QEROAw7HKRe+XFXLk+zWKMjJcYSFaRaGYTR1kgoLEWkGfA84BccU9Y6I/EVVD2S6c9nG0ywmTl/GjpIyfjzqiCz3yDAMIzuk4rP4B87ER38EHgYGAk9mslMNBW96VYA/z1qdxZ4YhmFkl1R8Foer6tG+z2+JyMJMdagh4RcWhmEYTZlUNIuP3AgoAETkBODdzHWp4WCywjAMwyEVzeIE4GoRWe9+7g0sE5HFgKrqkIz1LsuIZW4bhmEAqQmLUek+qYisA/bgJPdVqGqRiHTAqWzbB2ee78tV9Ut3+wnAte7231fV19LdJ8MwDCM+Sc1QqvoZsBtoC3T0/lT1M/e72nKGqh6jqt6MeeOBmapaCMx0PyMiA4ExOE72UcAj7qRM9YqZpAzDaMqkEjr7C+BbwGrcLG4yU6J8NHC6uz4ZmAX82G2foqqlwFoRWQUMBd5P8/kT4mVzG4ZhNEVSMUNdDvRX1bI0nleB10VEgUdVdRLQVVU3AajqJhHp4m7bA/jAt2+x21YNERkHjAPo3bt3GrsbqRNlGIbRFElFWCwB2gFb0njek1V1oysQZojIpwm2DTIABaZUu0JnEkBRUVFa065NszAMoymTirC4Dyd8dglQ6jWqaq3n5lbVje5yi4i8gGNW2iwi3V2tojsR4VQM9PLt3hPYWNtz1xbLuTAMoymTirCYDPwaWIxbprwuiEhLIKSqe9z1s3GmbJ0GjAXud5cvubtMA/4pIg/izKlRCMytaz9qSq4JC8MwmjCpCIttqvpQGs/ZFXjBzWHIBf6pqq+KyIfAVHcWvvXAZQCqulREpgKfABXAjapa7/NpmGZhGEZTJhVhMV9E7sMZ4fvNUAtqc0JVXQMcHdC+HRgRZ5+JwMTanC9dmGZhGEZTJhVhcay7PNHXlonQ2QZNyISFYRhNmFTmszijPjrS0CmvrLO7xjAM46AlaTyoiHQVkcdF5BX380DXr9Ck2La3jAoTGIZhNFFSSR54AngNJxIJYAVwa4b60+BonudUFqmsUrbsKU2ytWEYRuMkrrAQEc9E1UlVp+KGzapqBU5BvybBG7efxjdPPBSAvaUVWe6NYRhGdkikWXi5DPtEpCNu1rQ7t8WuTHesodCjXXNOO6wzAKXlZoYyDKNpkkhYeOE/P8AJm+0vIu/iTLN6c6Y71pAoyHNu04UP/48D5U1GqTIMwwiTKBqqs4j8wF1/AZiOI0BKgbOARRnuW4OhIDdSEX1vaQXN8uq9QrphGEZWSSQscoBWVC/k1yJz3WmYFORGFDDLtjAMoymSSFhsUtV7660nDRi/JlFZldZitoZhGAcFqfgsmjx+zaLChIVhGE2QRMIisE5TU8RzcINpFoZhNE3iCgtV3VGfHWnI+B3cplkYhtEUsenfUqBlgfksDMNo2piwSIGC3Bwe/rpTfNeEhWEYTRETFimSl+Pcqooqy+I2DKPpkTVhISI5IvKRiPzH/dxBRGaIyEp32d637QQRWSUiy0XknGz015v8yDQLwzCaItnULG4Blvk+jwdmqmohMNP9jIgMBMYAg4BRwCMiUu8p1N60qubgNgyjKZIVYSEiPYHzgcd8zaOBye76ZOBiX/sUVS1V1bXAKmBoPXU1TI5pFoZhBKDaNN4J2dIsfg/cgVv23KWrqm4CcJdd3PYewAbfdsVuWzVEZJyIzBOReVu3bk1rh01YGIYRy9R5G+g7YTpbdh/IdlcyTr0LCxG5ANiiqvNT3SWgLfCNraqTVLVIVYs6d+5c6z4GkRtyHdyVJiwMw3B4dn4xAB9t2JndjtQD2dAsTgYuEpF1wBTgTBF5CtgsIt0B3OUWd/tioJdv/57AxvrrroOnWVz1+Jwmo3YahpGYvBznvXD9k6mOfQ9e6l1YqOoEVe2pqn1wHNdvqupVOHNmjHU3Gwu85K5PA8aISIGI9AUKiUzMVG940VAAO/aV1ffpDcNogHgh9U2BhnSl9wMjRWQlMNL9jKouBaYCnwCvAjeqar3PQJTjExbrtu/L6LnO/t3bPPnBZxk9R0Nm4879PPbOGtPgjAbPrOXp9Y02ZLIqLFR1lqpe4K5vV9URqlroLnf4tpuoqv1V9XBVfSUbfQ1JRFh87c/vZ/RcKzbv5ScvLsnoORoy33/mI37532V8tr0k210xjLhU+YJdOrbMz2JP6oeGpFk0aMoqozO3MzXqrUpTtNWyTbvZtrc0LceqL1SVJZ/vorTCude79pdnuUeZZ/aKrby8sN5dcEYaKPdVc2jdLNHUQI2Dxn+FaWJ/WbTla/f+Ctq2yEv7eWKFUm059w/v0LFlPvN/MjItx6sPXvz4c27718Lw5z0HKur1/LsPlLPnQAU92jWvl/OVVVRx9d8c99uFRx9SL+c00ke5LzKytKKKvaUVFOSGGq0fo3FeVQY4qmdb+nduyXWn9AVg857MxFWnS1gAbD/IHPGffrEn6vPO/fXb/7MfnM3J979Zb+e78q8f1Nu5jPRTXhF5Vg+UVzL4ntf43tMLstijzGLCIkVaFeQy8/bTOeMIJ1fwywy9iMsq6i4sfjZtaY22v2/6svAIN6vEWOBKSus3juGLek6smv/Zlxk79uwVW5mzZnvGjn+wUVJWkfaE2nJ3YNciPydsOp3xyea0nqMhYcKihrRt7piedmbAnj7h+UUU/fKNOh/niffW1Wj7R2evYfaK9ER1PPzmSvqM/y8VtdCQqmL8QAcq6j3oLWs8M3c9Vz02J23Hu/pvc7liUtPVXP7+7lpGP/w/Xl3yBdMWbmTgT19jwvOL0noOzwqQI0JJWeP/rZqwqCHtXD9FupNwtu0t5Zm5G5Jv2MD586zVAOyrxcMTGzNQWt50ysFPeH4x/1u1LdvdaDT8/OVPWFi8ixuems/3n/kIgBc/rnsgwV9nr+FjN1vb81nsKU2fb23HvjLGP7eIfe4xt+4p5b7py2o1+Eo3JixqiKdZpJstu+seubTnQHm9/qjeXbWNzTGmm7xc5ye1rxYPUKyVoLQeNQvL6Yhmy54DrN1W+3yitz7dwvVPzmtQ97Vzq4I6H2Pi9GVc/Kd3gYgZ6gcjD6vzcT3+9r+1TPlwA0+5eVa3/3shj/oEVDYxYVFDWhVkJoBs+766C4ujfvY635/yURp6kxrfeGwOJ/xqZvjz+u0l7CxxzHO1ERYa47Q4UE+aRUlZBe+tjtj30xW+7PH60i/497zMaI2qyjsrt0a9lA+U113IDn/gLc747axa7/+df8zjtaWbWbA+c36ZmlKQV7fXXazg8/yLh3drzeFdW9fp2B7t3XyNjTv3A7CoeCcAElQhr54xYVFDRITvnzkAkfRWoN2+t24Oc28UPn3xF+noTlLKfRrMp1/sBuCiP/0v3FYb1dwTNB71pVnc/eISvuHzF5TXYjbEXSXxfVjjnpzPj56NtpfHG3Gn+pv601urWLZpNy989DnffHwu/3YL2gGc99A7KR0jEXUV1H07tQRgSgMyrTbLrds0OKUxwSfe3Db5OaGoCg/3vbKs1hpVO9dysc0NoPGeiYbgEzFhUQvatshHlbTamPccqP6yea8Gx69tAtvugPOmwq9f+TS8vqh4FxD9sq+NZvG5O5ryiH0408WmXdHnWfr57qjP5TWsLPyfRRs5+t7Xw6PAVHhtaXDUTHkKZsTdB8r5zWvLOfcP74SF0Ccbd/PIrFU8+f461mxNXzma2mopHdwR8pufbuGoe15jby1+D8s27U7rgCE/t26vu9h78bZb6iMvJxQuKAjw6NtrWL+jdtUHvP9/rL9uXz1HBgZhwqIWnFrYCYCFabQjBr0Yv/7YHCbNXp3wgfVGovFGtqu2ROcuLCreGWWLvvwvtStd8tj/1obXgzLFayMsYv0f/3j/s/C1r9u2j+1pyEj/14frOem+N6Ne7C0LokecNfX7vOsK9aUbd1f7Ll7RyZKy4PuTirD4/MuIsAv///eX88Cry/nJS9Fh01v31O2eXfPEh6xPoezKpl376TP+v7y3ehu7D5QzZ61TrWf7vjL2lFawZuveGp130679nPuHd/jZtE9q1e8gUrm3iYh9Rn/3xgrAqRvn1ywAJHBmheR4EValFZVRYfT7y+s3QTUIExa14LCurenUqiDqoa0r8UbRv5r+KX+dvSbwu/dXb6f/ndP5/jMfMfJ3swO3OevB6PaLHn43yhbtT4TbVVIezlT/aP2XKY8qg3JO9saMhFSVnSWJTW1BOSaPvLUKgNN/O4vjf/kGfcb/t04vwNfdEf2mXY5gWrVlDy1dP9T3zxzg9KPGLxXnxRBreVi3bR/H/WJG4B5tmgUHSsxekVyb/GJX9XyQFz76PHDbr0wMDsW+64XFTJm7HoAVm/fENX+9t3o7E15IHnI6Z40jHL7+1zl896nqkYIFNTQBeddTmwFZvFyluvpy/Pv7/Vo7S8rIjcnaLqus+bkqKqvC5zhQXhlVNeK9VdnPmTFhUUu6t22W1izuA+WViMDT151QrdzE3jij0PmfOQ/otDTVFjr63tc59w+z+Xznfr76yHspFzP80tVq/E64WM3it68v55h7ZyQUGEHzmwf5Pv714fqU+hXEzE+daVJaFeSyZfcBznpwNu+s3EbRoe3p2b4FUHMzlHfdfgd9SVkF5yfwHYTiPHk3/jN5BvCaOkQpeTw9Zz3jn1/Mmq17Oft3s/ndDGeUvLe0opqZLpSCd9WfI/NuwIutogZ+oKoq5YFXlwNEmXdS5e4XFwe2f1lSntSXoKqBjuwrJ33A+74giHt8ia+De7SlIMbElYrPZ/aKrVHmuQF3vcKvpjvm3dKKKvb5nnu/TypbmLCoJS0Lcigpq2Tz7gN8+sVulm2qboKoCaUVVTTLzeHkAZ144NIhUd/lxXmzNMtLbbSmqimPqtZtLwn7T1IN1/M0C38/Y23Uf3rLyb9I5FsJGt3+/d111dp++/qKlPoVi3/EWVZZxdm/j2hdlark5TovpvKYkamq8tLHn7OrpJwD5ZVhh76HhLdzlgfKKxn409fi5poUf1nCQzNX1eoadpWU84v/1M00438ZeiVh/rNoI7dPXcjge17jpPuiS56kUusomWO+JpUJ/Fr2+h0lbNlzgD+9tSrlUN54/qAd+8pYsTmxOeyuF5fQd8J0wPEDVVRWsXrrXt5fs53xz0eEkDeFwMNfP5ZeHVpU0xQTPW+7Ssp5b/U2rv7bXH7074Us3LCzWub3ouJdXPV4dJJmbH26X01fxjsr669EuhUSrCXN83J4a/nWqNDRdfefX+vjlZZXhkP7msWE+MXaQz1io4eCKMgNccmf32Pb3lLeuePMlPrijSQrU4joOLWwEyu37EVVyc0RvN9zPIdm0Etj4Yad/PSlJezNcOFAv9OxtLwq6v5t31sWfilu3n2Az3aUcFzvdrRulsfC4l3cMuVjhh/WmeIvS1izdR8rJ54b3t4beN/94hKOP7Q9XVonjue//sn5gf6NVPjY52u5+JhDapVo9uPnImYlLwpv3fYS1sXxTcT7/fk5kEQYJBMWd7+4mPdWbefNH54e9aL9sqScoROdZ+w3ry1n9a/OIyck7CutCJsPYwnSRnp3aMH6HSVJJy775xxHa/1k427Oe+gdbjpjAKe4Psog2rdwHPltmscKi6rwcfp3aRk2w+0treDoe18Pb7dqy15Gu3kbsXiBClcU9eJf8zbwyaZdHH9oB8ARzpNmr2HS7DV1eu/UBNMsaknz/LqF4fkpKatg8vufkeO+dWI1hoffCh6FfpnEBwBO3ZqP1u9kw479XDd5XtR38V7o3jzja7buq+Z0hoi99tazCjllQCfW7yhh94GKqH57Zqhn5q7noZkrw+1BvpmJ05exsHhXXF9BKsX9nptfzKotiUeN/vsVOzqvUg2//K+Y9AFj/zaXo372OlPnbeANd9S3qHhn+AH2NIy7Xlgc5av47WvL44Y5frJxN1VVmrR0fKLr2OiLGLtj1BF8MGFEwmPFUvxlCVPnRUwaNwT4F2LJCTBD7S+rjHLS704SjZfMD/TUB+tZs20fqonvz4rNe/j9GysYdM9rcYM6+nduFV7Pd/+nXjJtvMCCWDyt+uG3VjEmQdkUz/x0eNdWUe1/fWcN2/aWct5D73D43a/yv5WOL2psTA22VATx6Yd3Bpx5dDzNPJVnP93Uu7AQkWYiMldEForIUhH5udveQURmiMhKd9net88EEVklIstF5Jz67nMQyUxAk2avZvriTezYV1YtiueRWavCLyCA5xY4zjzPJNA85tjxVPx4mkX3ts3C635TyBvLolXdeBUy/YUI7w7wW3h5CHk5IVq4o7vSisqo7HZPEE14fjEPzoiYjYJCIeM5ez1iQ2oh4q/xuP3fCxn5u7cTHsdvXvKOGfY3aHA/7nh2UVhY+0uml1VU8eNnF/H0nPUs9wUJDOjSiv1xTBDnPfQOf3t3bVgYxyPWzOVnk+9e5OWE6Na2WY0StmoSr5+oVPtxv5jBwJ++Fv4clCc0rH9HXrrxZACeXxDsgI9l1oqtcYM1wIm8+/0bzuBj7rodPL+gmN0HyqOekSO7twHg/kuO4obT+gGR+SaunTwvpUKbxV+mFvrqaQxXn9SHey4cGG5/e8VWrng0Eml41eNz2LL7QLXikdsC7pvfsnD+kO6cM6gbnV1t9b7py9z9Iu+U+kp8zIZmUQqcqapHA8cAo0TkRGA8MFNVC4GZ7mdEZCDOXN2DgFHAIyKSvmF9LdkY8AJ7e0Ukk/ZX0z/le08v4LhfzOD4mOKAD7y6nOv+ERnlF8TYhAt8wsIbuQRlFQeV8H7pxpO5+qQ+4c+J1P8Fcaqezl0XeRG/HTBtpOcAzs8JhfteXqmUV1ZxybE9OKxrq7ihs0H1nvxC5uyBXTnDHUklYvWWiP3au+fJrGalAaPbXq5Tu0o1/EDGw/9CenZ+MQvd/JJiX1RcQW6o2gvZP3hctmlPVOn4aTedXO08/hDPx/+3NhySCrDDN6L0zC3JrvvVJZFEzRWb9yTYMprJ1wwFoH3L6kI0ViBuDdAGVCO5DfGitSA6mu6NJFVbX/f5I77zj3n8YOpChvzsdX41fRmlFZVMmbue0ooqOrbMZ8zQ3uFnqUV+xGSVSqHNDSlGOnqm41BIuPT4nlHfrY7Jdwk6ZpAW5X8eHhpzLKGQ8OdvHAdEBoh+4XzJI++l1Ne6Uu/CQh08PTvP/VNgNDDZbZ8MXOyujwamqGqpqq4FVgFD66/HwQSFL47921wmx/khvrtqG68u2RRo1ol98Pxq//dHFAJORNSs5VuiHvYv91XXLA5p15zcFFRbiLxkD/FpIrFUj/Ko5JRfO2ahvBwJO4XLKqqoqHT8Fq0KcuMmEQWZoTq1jkxJeUS31vz928n/vX6zRqrRS0GC0/MvqJJUWPi5z5eU6C9tvnl3aZSpo1ubZtx2VqR20HMLoqNaDu3YMm4/l27cFTaXPTffednu2h85diLH828vOzq8fsNT83lm7np2lpRx0z+dcjCtUyhbM6BLK3q2b55SQcetew5QdGj7au3+RLhYDXnqvA08MmsVk96JhIYn8ynEm5v+pY838vCbqxj//GJe+vjz8L3xBGqLFMzGfsd/MpOmh//5aN0sjwe+NiTutl/7c2ov9Zb5udx13pH89IKBYTNVUZ8OnNC3Q/j+ZGMWzKz4LEQkR0Q+BrYAM1R1DtBVVTcBuMsu7uY9AH/NgGK3Lei440RknojM27o1s1EC8X7US90Iili+8dgcbnhqQZRDfMYnm5m+eFM4DG/+3WcB0KlVPif07cDT150QfoHtKinnW3//kLN9KnpQZFF+irN09Rn/37BT8sR+HePOIRyKETzFX+4Pj27yckPk5zgPYXllFRVVVeTmhGhZkMue0uD5A4LMUH7NKvZ88diwoyScWJdqXkSQsPBGhlWqtEnD1Jj/mreBpz9wnKSTrxnK+xPOTDijYtD/q6yiitKKSs5/KFI+xXvZrd8eGa3mui9CT+C9ftvw8HexTt4Jzy+O8k8ckmQ2wCE92wKOuTWeWc1jy54DLFi/M1yR2aNT64Ko6/NCQTfvPsCbn27mjmcX8cCry6O00NoWL8wNSTj/pqSsMjyI8YRGPGGhqlz+6Pv85e3VUQOZVKMbY/NHLv9Krxr3PZYWBTl8Z3g/rnEnWvPo1KqAbW4NuVjzlfdc/f3dtRT98o2UfTM1ISvCQlUrVfUYoCcwVEQGJ9g86O0ROJRU1UmqWqSqRZ07Jzdl1IWJXz2K3h1aMOmbx0e17y2tYPHnu1I6xnf+MY97X444Wr0SCbk5If51/UmcPKBTuFbMkoBjBjm5cnMkZRu29zKvUo3raIvVUvwvoYpKDX8uq6iivFLJCzmaxcINO+l/5/RqxwvSLPzmoVS1okdnr+Gih98NnzsVgrbzIr+q1Kn7deMZ/VM6ViL+u3gT4Nj8RSTKBBJLUAmKA+VVrIwJ8WyWF+Kyv7wXNn1BJFT5jlFHcHlRTwZ0jnayxvKJLwIryLTksfCnZ/PvG04CnIHLlpgkyNg8hB/+e5Er4KLv78SvDo4ykT385ioqKqs44VczueaJiBnWP0qOPVeq5OZI1Lk8IeHd39hnwjPrrti8l7lrd3D/K59yxE9erfF5YzVvgNvrWIV2087g/K12LfJYs3VfYDUD7/eyY18Z2/eV1rkOVhBZjYZS1Z3ALBxfxGYR6Q7gLre4mxUDfnHdE8j6DPcXHn0Is+84o9oI7bWlX/DVGtgQ/SYMCXjLt3ND874b44wuragMdFbWZP7f43q3c459+gDOHdwtcJvYCYn89vSte0rDD+OdLyxm1/7ysGYRj1umfBzlfykpq+DRtyNmiBz3BRg7So3HzGWbOdGnrQWxasse5q7dEe67PwDAC1TwIk6uHNo7pfMCXO86T4Po1qZZuJjeRQnm1w4S0nsOlEe92MEJVPhwXbSPydPCLj2+Jw9cenRSrWy3z0HfvW18zaJti7zwiLl3hxbM/+zLKLPrz30DnO17S8MTZ8X6qdo0y6NXh+b07+zch0mz13BpQHmZJb7aXJ7G/rXjou3/t7jm2HjkhiQqKTI/bIZylmUV0b9jT1tK1ZEdz9kfVMn25iR9TUa86ZDDv9XfzmL73jI6ty7gqB6OBnjBHx0tdM+BCloV5KasodeEbERDdRaRdu56c+As4FNgGjDW3Wws8JK7Pg0YIyIFItIXKAQawBygDrEly9Nc3Tru/Bnxwgb9I/9vDeuT8NgFuTkcf2h7Du/Wmp9cMJCrTox+UV5zcl++LCln4879bN9byiWPvBsVdlmpGn4ovWKCuTlC+yQv+v3llby8cCPrt5fw9b9GJx55msXLN53CDaclH+W/vHBjUjPUWQ/O5vJH3w9vN/maofzlKkcjvGBId2b/6AwmftVRbnu2b8HlRT3jHsuPf36Ekwd0jPqua5uCsCAI0h5ev204vxg9CHCEzg/PjoxGdx+oYFlMRJQX/5+M84d0B6LDR4Pw18P6+gnxBaQXAHDm/80Kt/kdxP5qun6zo6eZiAg/PPvwcHtQomdQ0b3fXhZt+082eMgJRWsWXjSeN/KPrQvl5XIEaTJB2u274yM5Sn5BlorZN5m26kU/Hh/g8/Hjv75te0vp3KogKoH3u0/NZ8+BiqTRhbUlG5pFd+AtEVkEfIjjs/gPcD8wUkRWAiPdz6jqUmAq8AnwKnCjqma/BKNLqiPgdB5/zprtvLY0uBS5iHDWkV0BqkVnxPJlSVnYlpubE6qmJR3rah7rtu1jxea9LFi/k0m+OlU3nTGg2oswLxQK/NF38PlEtuwp5eZnPmL4b96q9vLwRoe9OrRg/LlHBPoR/AIxJ17djAA8M1T3ts0YNbgbi352NqOP6UHvji2ibM8j3PvnJ3akC9EP701nRI8mYwcN5wyKPuZhXVvzTTdqbcK5R3LTmZH9n3hvXWDmeir86evHse7+8xncoy1PXhscKPDrrx1F19YR7SrRpECelhgv5Nafq+Mv1/KVPh3C681qkZMUq2VXVimPjy2Ku31uKBRlm/aEQX5Ys4gWFl45kaCAk5EDq////fzm0iFcUeQYO2JrQnkM7Ru5/k5JJl3q52peA92Q33j4tfxt+8ro2CqfI7pF5tF4ZckXbN9XmrE5d7IRDbVIVY9V1SGqOlhV73Xbt6vqCFUtdJc7fPtMVNX+qnq4qr5S331ORLsW+fz0goHJN0xCkO0TgjWLKyZ9EK4uevOZA6rZSPt0ahl+YQzr37Ha/h5b90T/sGJHVIe0c14oe0ormOcLpwV48tqhtCzIrWb2UpR+AaNaf+7IloAH1CM2sinI3u/fJpm2HTQpkNfneCOwUws7cdaRXRnqvvCuPulQfnHxoGrb+TPcO7TM59COLXx9jH45PfrN+C+6TNEnINIK4LTDuvCd4RET2ogju4TXf3TO4VHbHul7gT03v7haYIc/27qySvnPzafw9299JfqEcbTtVIMxvGN7Jtkglm/ew7O++kme3AqboSqrWPCTkdx53hGAE4gAwVFPd18wkEuO68GT1w7l0I4t+N7pjmbw84sGcXSvdoRCwn2XHMWnvxgVtz9PXXtCWAh0DBAWfu21dwfnd5ObI1w5tDd/vPLYwGP6hcW6bfs4pK3jE/OHmi8u3hXOKUk3lsGdBo5yI0diefHGk6PU13j89eoi3vrh6YHfNcvLqZak5+ecQd0S2kgnXV3EiCMiL4OzfaOm7fvKooTRGYd3idrXGxF9uHYH/+dLrIPIgx6rWTTPy4nyCYS3920XFJPvEeskjY1iOT0mByNZgTV/QTevbHayl1SL/FweG1vEYd0codemWR4t8nOrhYZWxQiLZ28YFvZ51KbC6T+vOyHqc68OzeP+LlIhXpWBZnmhqKTSIT3bhddvPGNA1LYn9e8Y1lRv//dCbomZiXHjzv2cdpjzPzl3cHcG92jLGUdE/47izUnRo33iiCw/lapJzZt+vEGC97srq6iiQ8t8eneICNBdJeX8Z9Gm6v1q15wHLz+GUws78/aPzuCOUY6AGTusTzjJMBSShIm5+bmhcHiyXzv2BoVdWjcLh6x71oM2zfK475KjuDCOj8tv5tu1v5xTD3PKkLT3CdHt+8pqFAJeE0xYpIF4L5+OLfOrqfhz7hwRdnx6jBzYNWEo47UxIXR+vB/aO3ecwT+uqW52aFWQy6BDnNHhrWcVcseo6JGjv6ZNYdfWUXVmWrsj748C7MzeXNuxmkWL/FxaN8ur9pL3a06JSozHmgu8ke0R3Vrz8k2n8GhM9Fky/ImL77glF1J1/nlC2hupfftk5//w4OVODsNIn7mqfYs8Orcu4LaRjuBOlPvh+UtiGTYgugbRgM6tqv1WAJ694STm3pW8zEe8QYb3kvNrQokY4hsMxZpttu0to1/nliz86dncfOaA2F0BOK53tJCdc+cIOrbM51dfPSrheTu1irwEq6o06qWYDM+MNdh1AHv+O3929PVPzau2Xzrx+uDX3pfdO4pfXDyYm0cMYPI1Q935u53tOrZKfH1XxITlehrJ987oH2Vd8J73dGOFBNNAvBm4WjfLjfpu3PB+dG3TjI4t82sUTx5bWNCP9wD16tCCXh2CH37vtSVItQcuUW6B9yMPihjxBGSHmPwMz1b9yb2jmL54U7ikyDmDunHRMYfwwKvLqwmLE/p2YMXmPXxZUl7NfPPApUO48OjujBrcPW4//RT98g0euvIYhvXvxK795XWaoCoiLByhef6Q7pw/xBGml7g+jMnXDGXzrgNh23WX1s2458KBUTb7WEbFiTwD57fkCcxFvjBZP+1b5tOldfxESo94I19PcL96y/Dwuf545bFxQ779gj+oXHmrgtyEuSRd2jTjlhGF/MGtEda1TTPm/2Rkte2m3XRyOBwa4NkbhvHCR5+zfkcJVw/rQ8sYk6RXHDAIr58dWuZHDYD8vqkP1uyott/LN50S9zpqijcmEXF+S9ed2pdQSPjmiYcCzuCssGtr/uSWkzksyTzeQ3q24w9jjuGWKR8D0M3VTAZ0ac3NI1qzsHgnbyzbwllJfC61xYRFGog3A1dsCOmd5x0JODka5/jKYycj0fSiqWWmOsuQVPeBJLID5+eGKMgNsXl3dU3AE4IdWuZzwZDuVFYpD1w6JPxiBRg1qBs/OudwerZvzoVDDnFfRst5ZNbq8DYnD+jI09edyOT31nHPtKXVNIuWBbkpCwpwokTuffkTXr11OBf/6d1qQjmZs9GP55ht0zz+Y+KZYPx4Gkgst488LNCf4+fS43uGI5+8EMrLi3pGRaHF82/FEi93xhvxNs/PCZuqLjz6kLjmD38CWFBYZ6JQaY9Rg7uFhUUQxx/aPmwO88wofTq15LaAnIWckFBZpVx/Wj/ueiF4zpV4lsZEA6+7zz8yrkm5Ntw28jDG/WM+hV1bsyyBf+P64f0Y2rdDwgGGh19b7NQy+rf8yDeO5/Od+wO10XRgwiINDDqkLdef1o8PVm+PTpqK84s9vFviEUQs8ZLOrh/eLzA3IxbPth4KSbXojYEBKuvvrjg6nKXdulkupQHFzvxnffjrxwWeNxSSKBt4kOPNi8rx7lVdp74EZ/a/hRt2Bmpv029JfeQYq1nUlVTi739+0aCwsHji246j+O4LBkYJi1QqlXqIJK8dlYyLjj6Ev7ztCPggE2Iqjuoju7dh/t1nxY2q+te4EwH46CcjwybOIF743jC6t20eHlXHFRZxnosgwTagSyv6dWrJZUV1z772c2ph54RCwiM3J5SSoIBIwEenVgXVzKn5uaGMCQown0VayAkJE849kp5xzEBB3HneEQzt04FZKTgwv3t6/6gQOY9U8hCgejDKGz84Lbx+aECfv3psz/DoON6LMpG2E48ubSKmkxW/PJdHvnEc91zkRBmFHZEpCIv/u+zoaklSscUHv+mbOMZfBykV842HJyzSUQYkVfwDjOGFzjXFZuPGy70J4oIhEW1hSC1HzUEDCj/xSt3H0rFVQTVT6a+/dhS/v+KY8CCmfcv8hKGfx/ZuHxYUfn5+UXS0Wjy/VGGXVpx5RGwgRz6Tri6q0X3NFs3znfuUaqWDdGKaRRpZFWcWrv+77OhqEQrjhvdn3PDUXvbtWuRz9Ul9uPOFxfTv3JIbTuvPWUd2pX2cek6xeFE8x/ZqBzgjKY9kJgTPzHVqYSd+ftEg3lq+lVnLt1DYNbE5JQj/SyA/N8R5R0XMSxcM6c77q7fz4xgHfBBfO74n3ds24+uPRQRC7P31ZysPG9CR15ZurnEpjy5tCgiJY2PPBt4Lz59XUtOJbvzZ8s9/d1jakkZ/MPIwKiqreOjNVXXSXK74SuoZ84mIFWiJTHAnD+jEm59uCbd9v44Z1/VJ8zznGaqJdpkuTFikkZP6d2S5WxXWP9L5WpLkuFQ4zR05P3Dp0UkzPWMZcWRXFvxkZDVnNCS3f3smsJEDu9Kvcyv6dW6VMDorGf07t+SUmKgfcJyx/3f50QF7BDNsQCfW3X8+fcb/F0gcffTrrw3hD2NyahTXD3D6YV2YefvpSYvuZZpUTI3x8MItTy3sFDeBrDZ87fiedHZNId8ZXvvfQ125fng/Pt6wkyO6taZdizxGDerGlA83xDVDAZw7uFu4mu/fv/0VhvWPPxNeQ8PzuZiwOMg5Z1A3nnhvHe1b5DE2SamNmtKjXfM6TZ8YJCgg+YvIMzEkqiVUE2befnpajhOL33zV3Fcpddzwfgmd+IkIhSSjNuCa8Nx3h6Uc6urH0wy/fXKfOp3/7IFded0310TzvBzyc0PcelZ1B3R9MsENGgHH37F22z6mfLiBVglMh4e0a86w/h15b/X2hEKlIeI9r6kEFaQb81mkEa+oWLLJ6w8merqJU52SxIBnG/8seP4X/LpalrxuaBx/aPsaRXJ53HPhIG4feRinH9Yl+cYJeCgmqzhRVFG2EHGE+w9GHsYj3wgOuvDwJ+sdTPTp2IKbzhjAo3FydTKJaRZpxDPp1DX6pCHxp68fx9/fWxeubtnQuO+So8Kj59cDZlk7/fC6vSQPdtq2yKtzFVRwzIQn9usQzk3IRAnsdCAiKfkgxnylF7OWb+WI7jWLTMw2IsIPz0nu18sEDW94cBDjJUFVNiJp0aVNM3486oi02rvTyZVDezP6mB6MPqZH2BHsWRZ+esFArhya3nDIpsyUcSeF1zNRArs+GTW4O+vuP5+e7Wtu2muqmGaRRsLzZR8EwuKDCSMOOhU8Gfk5IcorKzmud3uWbtzNwEPa1Mk5nC2+c2rftPmI0s0Zh3fmrYB52Y3GjwmLNOJpFgeDyyIoVv1gxxvtfvOkQxk77FAGdDm4TAwed51f9yrGmWLS1UW1yrExDn5MWKSRiM/iIJAWjRAvnDA3JEnLahi1Iy8nVKPZGI3Gg/3X04inWWQjrM2IZLXm1mBCJMMwUiMb06r2EpG3RGSZiCwVkVvc9g4iMkNEVrrL9r59JojIKhFZLiLn1HefUyUvJ8Td5x/JszeclHxjI+14lUZzcw4+P4VhNHSyMQSrAG5X1SOBE4EbRWQgMB6YqaqFwEz3M+53Y4BBwCjgERFpmHF7wHWn9jtobeUHO56wMCOgYaSfbEyruklVF7jre4BlQA9gNDDZ3WwycLG7PhqYoqqlqroWWAUETy5sNGm8qrbmMzKM9JNV47qI9AGOBeYAXVV1EzgCRUS8bKoewAe+3YrdtqDjjQPGAfTunZ4CZcbBw9++9RVe/OjzahVpDcOoO1nzBIpIK+A54FZV3Z1o04C2wKGjqk5S1SJVLercufqkNEbjpleHFtw8ovCgzK0wjIZOVoSFiOThCIqnVfV5t3mziHR3v+8OeDWEiwF/Gm5PYGN99dUwDMPITjSUAI8Dy1T1Qd9X04Cx7vpY4CVf+xgRKRCRvkAhMLe++msYhmFkx2dxMvBNYLGIfOy23QncD0wVkWuB9cBlAKq6VESmAp/gRFLdqKrBczMahmEYGaHehYWq/o9gPwTAiDj7TAQmZqxThmEYRkIs1dUwDMNIigkLwzAMIykmLAzDMIykmLAwDMMwkiKNtTSCiGwFPqvl7p2AbWnszsGAXXPToKldc1O7Xqj7NR+qqtWymhutsKgLIjJPVYuy3Y/6xK65adDUrrmpXS9k7prNDGUYhmEkxYSFYRiGkRQTFsFMynYHsoBdc9OgqV1zU7teyNA1m8/CMAzDSIppFoZhGEZSTFgYhmEYSTFh4UNERonIchFZJSLjs92fdCEivUTkLRFZJiJLReQWt72DiMwQkZXusr1vnwnufVguIudkr/d1Q0RyROQjEfmP+7lRX7OItBORZ0XkU/f/fVJjvmYRuc39TS8RkWdEpFljvF4R+ZuIbBGRJb62Gl+niBwvIovd7x6SmswUpqr25/htcoDVQD8gH1gIDMx2v9J0bd2B49z11sAKYCDwADDebR8P/NpdH+hefwHQ170vOdm+jlpe+w+AfwL/cT836mvGmb/+Onc9H2jXWK8ZZ3rltUBz9/NU4FuN8XqB4cBxwBJfW42vE2cuoJNwKn+/Apybah9Ms4gwFFilqmtUtQyYAozOcp/SgqpuUtUF7voeYBnOgzYa5+WCu7zYXR8NTFHVUlVdC6zCuT8HFSLSEzgfeMzX3GivWUTa4LxUHgdQ1TJV3UkjvmacaRaai0gu0AJnFs1Gd72qOhvYEdNco+t0ZyBto6rvqyM5/uHbJykmLCL0ADb4Phe7bY0KEekDHAvMAbqq6iZwBArQxd2ssdyL3wN3AFW+tsZ8zf2ArcDfXdPbYyLSkkZ6zar6OfBbnMnSNgG7VPV1Gun1BlDT6+zhrse2p4QJiwhBtrtGFVcsIq1w5j6/VVV3J9o0oO2guhcicgGwRVXnp7pLQNtBdc04o+zjgD+r6rHAPhzzRDwO6mt2bfSjcUwthwAtReSqRLsEtB0011sD4l1nna7fhEWEYqCX73NPHJW2USAieTiC4mlVfd5t3uyqprjLLW57Y7gXJwMXicg6HJPimSLyFI37mouBYlWd435+Fkd4NNZrPgtYq6pbVbUceB4YRuO93lhqep3F7npse0qYsIjwIVAoIn1FJB8YA0zLcp/Sghvx8DiwTFUf9H01DRjrro8FXvK1jxGRAhHpCxTiOMYOGlR1gqr2VNU+OP/LN1X1Khr3NX8BbBCRw92mEThz1zfWa14PnCgiLdzf+Agcf1xjvd5YanSdrqlqj4ic6N6vq337JCfbXv6G9AechxMptBq4K9v9SeN1nYKjbi4CPnb/zgM6AjOBle6yg2+fu9z7sJwaREw0xD/gdCLRUI36moFjgHnu//pFoH1jvmbg58CnwBLgSZwIoEZ3vcAzOH6ZchwN4draXCdQ5N6r1cDDuFU8Uvmzch+GYRhGUswMZRiGYSTFhIVhGIaRFBMWhmEYRlJMWBiGYRhJMWFhGIZhJMWEhWHUERG5y618ukhEPhaRE0TkVhFpke2+GUa6sNBZw6gDInIS8CBwuqqWikgnnGqv7wFFqrotqx00jDRhmoVh1I3uwDZVLQVwhcOlOLWK3hKRtwBE5GwReV9EFojIv906XYjIOhH5tYjMdf8GuO2XuXM0LBSR2dm5NMOIYJqFYdQB96X/P5zy2G8A/1LVt92aVEWqus3VNp7HyaTdJyI/BgpU9V53u7+q6kQRuRq4XFUvEJHFwChV/VxE2qlTatwwsoZpFoZRB1R1L3A8MA6nPPi/RORbMZudiDMhzbsi8jFOHZ9Dfd8/41ue5K6/CzwhIt/BmZjLMLJKbrY7YBgHO6paCcwCZrkawdiYTQSYoapXxjtE7Lqq3iAiJ+BM3vSxiByjqtvT23PDSB3TLAyjDojI4SJS6Gs6BvgM2IMzhS3AB8DJPn9ECxE5zLfPFb7l++42/VV1jqr+FNhGdMlpw6h3TLMwjLrRCvijiLQDKnCmsBwHXAm8IiKbVPUM1zT1jIgUuPvdjVPhGKBARObgDN487eM3rhASnIqiC+vjYgwjHubgNows4neEZ7svhpEIM0MZhmEYSTHNwjAMw0iKaRaGYRhGUkxYGIZhGEkxYWEYhmEkxYSFYRiGkRQTFoZhGEZS/h+Ml8Jy0ttFAgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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 final structure of the simulation. We take advantage of the fact that the TIP3P water model is a rigid water model which means the nearest neighbors, i.e. the bound H atoms, of each O atom never change. Therefore they need to be indexed only once."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "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": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEMCAYAAADDMN02AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAPR0lEQVR4nO3df7DldV3H8edLVuKnAu6VIXS7VEQiiuKlDDQVclJQCQf8kb9ymLbGIlIzaczslw7MUAPiD2YjAosBEzFFDFMUsBHRXSQQN4OUxU1yF23U1BlafPfHObiXy3LvWdjv93v2fp6PmZ2953u+d7+v/czu637P53zP55uqQpLUjkcMHUCS1C+LX5IaY/FLUmMsfklqjMUvSY2x+CWpMSuGDjCJlStX1uzs7NAxJGmnsm7dururambh9p2i+GdnZ1m7du3QMSRpp5Jkw7a2O9UjSY2x+CWpMRa/JDXG4pekxlj8ktSYzoo/yQVJNiX50rxt+yX5RJLbxr/v29XxJUnb1uUZ/4XA8xZsOx24uqoOBq4eP5Yk9aiz4q+q64BvL9h8AnDR+OuLgF/r6viSpG3r+wNc+1fVXQBVdVeSxz7YjklWA6sBVq1a1VM8tWj29Cvv9/iOM44fKInUj6l9c7eq1lTVXFXNzcw84BPHkqSHqO/i/2aSAwDGv2/q+fiS1Ly+i/8jwGvGX78G+HDPx5ek5nV5OeclwPXAIUk2JjkFOAN4bpLbgOeOH0uSetTZm7tV9fIHeerYro4pSVra1L65K0nqhsUvSY2x+CWpMRa/JDXG4pekxlj8ktQYi1+SGmPxS1JjLH5JaozFL0mNsfglqTEWvyQ1xuKXpMZY/JLUGItfkhrT983Wpak3/+bri914fdL9pGnjGb8kNcbil6TGWPyS1BiLX5IaY/FLUmMsfklqjMUvSY2x+CWpMRa/JDXG4pekxlj8ktQYi1+SGmPxS1JjLH5JaozFL0mNsfglqTGDFH+S1ye5NcmXklySZLchckhSi3ov/iQHAr8HzFXVYcAuwMv6ziFJrRpqqmcFsHuSFcAewDcGyiFJzem9+Kvqv4CzgDuBu4DvVNW/LNwvyeoka5Os3bx5c98xJWnZGmKqZ1/gBOAg4CeBPZO8cuF+VbWmquaqam5mZqbvmJK0bA0x1fMrwNeqanNV/R9wOXDUADkkqUlDFP+dwNOT7JEkwLHA+gFySFKThpjjvwG4DLgRuGWcYU3fOSSpVSuGOGhVvQ142xDHlqTW+cldSWqMxS9JjbH4JakxFr8kNcbil6TGWPyS1BiLX5IaY/FLUmMsfklqjMUvSY2x+CWpMRa/JDXG4pekxlj8ktQYi1+SGjPIevzSzmL29Ct//PUdZxw/YBJpx/GMX5IaY/FLUmMsfklqjMUvSY2x+CWpMRa/JDXG4pekxlj8ktQYi1+SGmPxS1JjLH5JaozFL0mNsfglqTEWvyQ1xuKXpMZY/JLUmEGKP8k+SS5L8u9J1if5pSFySFKLhroD1znAVVV1UpJdgT0GyiFJzem9+JM8Cvhl4DcAquoe4J6+c0hSq4aY6vlpYDPwd0m+mOT8JHsOkEOSmjTEVM8K4Ajg1Kq6Ick5wOnAW+fvlGQ1sBpg1apVvYfU8jL/pungjdPVtiHO+DcCG6vqhvHjyxj9ILifqlpTVXNVNTczM9NrQElaznov/qr6b+DrSQ4ZbzoW+HLfOSSpVRMVf5KjJ9m2HU4FLk5yM/AU4B0P48+SJG2HSef4z+WB0zHb2jaRqroJmHso3ytJengWLf7xB6uOAmaSvGHeU48CdukymCSpG0ud8e8K7DXeb+95278LnNRVKElSdxYt/qq6Frg2yYVVtaGnTJKkDk06x/8TSdYAs/O/p6qO6SKUJKk7kxb/B4DzgPOBe7uLI0nq2qTFv6Wq3ttpEklSLyb9ANcVSV6X5IAk+933q9NkkqROTHrG/5rx72+at60YLbgmSdqJTFT8VXVQ10EkSf2YqPiTvHpb26vqfTs2jiSpa5NO9Rw57+vdGC2sdiNg8UvSTmbSqZ5T5z9O8mjg7ztJJEnq1ENdlvkHwME7MogkqR+TzvFfwegqHhgtzvYE4B+7CiVJ6s6kc/xnzft6C7ChqjZ2kEeS1LFJ5/ivTbI/W9/kva27SNJ0WnjfXmlnNekduF4CfB44GXgJcEMSl2WWpJ3QpFM9bwGOrKpNAElmgE8yulG6JGknMulVPY+4r/THvrUd3ytJmiKTnvFfleTjwCXjxy8FPtZNJElSl5a65+7PAvtX1ZuSvBh4BhDgeuDiHvJJknawpaZrzga+B1BVl1fVG6rq9YzO9s/uNpokqQtLFf9sVd28cGNVrWV0G0ZJ0k5mqeLfbZHndt+RQSRJ/Viq+L+Q5DcXbkxyCrCum0iSpC4tdVXP7wMfSvIKthb9HLArcGKHuSRJHVm0+Kvqm8BRSZ4DHDbefGVVfarzZJKkTky6Vs+ngU93nEWS1AM/fStJjbH4JakxFr8kNcbil6TGWPyS1JjBij/JLkm+mOSjQ2WQpBYNecZ/GrB+wONLUpMGKf4kjwOOB84f4viS1LKhzvjPBv4Q+NGD7ZBkdZK1SdZu3ry5t2CStNz1XvxJXgBsqqpFF3mrqjVVNVdVczMzMz2lk6Tlb4gz/qOBFyW5A7gUOCbJPwyQQ5Ka1HvxV9UfVdXjqmoWeBnwqap6Zd85JKlVXscvSY2ZaHXOrlTVNcA1Q2aQpNZ4xi9JjbH4JakxFr8kNcbil6TGWPyS1BiLX5IaY/FLUmMsfklqjMUvSY2x+CWpMRa/JDXG4pekxlj8ktQYi1+SGmPxS1JjLH5JaozFL0mNsfglqTEWvyQ1xuKXpMZY/JLUGItfkhpj8UtSYyx+SWqMxS9JjbH4JakxFr8kNcbil6TGWPyS1BiLX5IaY/FLUmMsfklqTO/Fn+TxST6dZH2SW5Oc1ncGSWrZigGOuQV4Y1XdmGRvYF2ST1TVlwfIIknN6f2Mv6ruqqobx19/D1gPHNh3Dklq1aBz/ElmgacCNwyZQ5Jakqoa5sDJXsC1wNur6vJtPL8aWA2watWqp23YsKHnhNrZzZ5+5SDHveOM4wc5rrRQknVVNbdw+yBn/EkeCXwQuHhbpQ9QVWuqaq6q5mZmZvoNKEnL2BBX9QT4W2B9Vf1138eXpNYNccZ/NPAq4JgkN41/HTdADklqUu+Xc1bVvwLp+7iSpBE/uStJjbH4JakxFr8kNcbil6TGWPyS1BiLX5IaY/FLUmMsfklqjMUvSY2x+CWpMRa/JDXG4pekxlj8ktQYi1+SGmPxS1Jjel+PX1ruFrvXr/fj1TTwjF+SGmPxS1JjLH5JaozFL0mNsfglqTEWvyQ1xuKXpMZY/JLUGItfkhpj8UtSYyx+SWqMxS9JjbH4JakxFr8kNcbil6TGWPyS1BiLX5IaM0jxJ3lekq8kuT3J6UNkkKRW9V78SXYB3g08HzgUeHmSQ/vOIUmtGuKM/xeA26vqq1V1D3ApcMIAOSSpSUPcbP1A4OvzHm8EfnHhTklWA6vHD/83yVd6yLaYlcDdA2eYFo7FVts1FjmzwyTD89/FVtMyFj+1rY1DFH+2sa0esKFqDbCm+ziTSbK2quaGzjENHIutHIutHIutpn0shpjq2Qg8ft7jxwHfGCCHJDVpiOL/AnBwkoOS7Aq8DPjIADkkqUm9T/VU1ZYkvwt8HNgFuKCqbu07x0MwNdNOU8Cx2Mqx2Mqx2GqqxyJVD5helyQtY35yV5IaY/FLUmMs/gUmWU4iybOT3JTk1iTX9p2xL0uNRZJHJ7kiyb+Nx+K1Q+TsWpILkmxK8qUHeT5J3jkep5uTHNF3xr5MMBavGI/BzUk+m+TwvjP2ZamxmLffkUnuTXJSX9mWYvHPM8lyEkn2Ad4DvKiqngic3HfOPky4tMbvAF+uqsOBZwN/Nb5Sa7m5EHjeIs8/Hzh4/Gs18N4eMg3lQhYfi68Bz6qqJwN/wZS/yfkwXcjiY3Hf/6MzGV3MMjUs/vubZDmJXwcur6o7AapqU88Z+zLJWBSwd5IAewHfBrb0G7N7VXUdo7/bgzkBeF+NfA7YJ8kB/aTr11JjUVWfrar/GT/8HKPP6SxLE/y7ADgV+CAwVT1h8d/ftpaTOHDBPj8H7JvkmiTrkry6t3T9mmQs3gU8gdEH8G4BTquqH/UTb6pMMlYtOgX456FDDCXJgcCJwHlDZ1loiCUbptkky0msAJ4GHAvsDlyf5HNV9R9dh+vZJGPxq8BNwDHAzwCfSPKZqvpux9mmzUTLkLQkyXMYFf8zhs4yoLOBN1fVvaMXxdPD4r+/SZaT2AjcXVXfB76f5DrgcGC5Ff8kY/Fa4IwafRjk9iRfA34e+Hw/EaeGy5DMk+TJwPnA86vqW0PnGdAccOm49FcCxyXZUlX/NGgqnOpZaJLlJD4MPDPJiiR7MFpZdH3POfswyVjcyeiVD0n2Bw4BvtpryunwEeDV46t7ng58p6ruGjrUEJKsAi4HXrUMXwVvl6o6qKpmq2oWuAx43TSUPnjGfz8PtpxEkt8eP39eVa1PchVwM/Aj4PyqWvRyrp3RJGPB6KqNC5Pcwmi6481VNQ1L0e5QSS5hdNXSyiQbgbcBj4Qfj8PHgOOA24EfMHoltCxNMBZ/AjwGeM/4THfLNK9S+XBMMBZTyyUbJKkxTvVIUmMsfklqjMUvSY2x+CWpMRa/JDXG4pekxlj8Ug+SnJvkxiRHDp1FsviljiXZE3gs8FvACwaOI1n80n2S/GmSP3gY3z+b5IdJbpq/fbyu0wHANcA7x/vuPr6Zzz1JVj6M2NJ2s/ilHes/q+op8zckeQywB/A94F6AqvrheL9mF3PTcCx+NS3JW8a3l/wko0XmuvDHwFnArYzuZiYNyuJXs5I8jdGqo08FXgzs8Ddek8wCRwHvZ7SK6xN39DGk7WXxq2XPBD5UVT8Y3zzmx8tOJzk5yTlJ3pXk7eNtH573/AfG91Ndyl8Cfz6+Z4HFr6ngssxq3QOWp01yNDBXVaeNH5+X5FnA/DX2H1FV9y72Byd5CqNXEs9I8m5gN0a3qJQG5Rm/WnYdcOL4Cpu9gReOt58CnLtg3yOAQ8c/BC5isjdlzwReOO9mHIfjGb+mgGf8alZV3Zjk/YzuG7wB+Mz4qUcyfiWQ5CBG1+BvBt5YVV9Icjwws9ifneQYYM+qunre8b6ZZM8k+1XVt3f4X0iakDdikRZI8iTgLcAmRj8E3gr8DfDSqronyZ8Bl1XVLQu+bxb4aFUdth3HuoPRtNKyu3OZppfFL+0gSR4PfBb41sJr+bex7+7A9YxeOTzJVwDqk8UvSY3xzV1JaozFL0mNsfglqTEWvyQ1xuKXpMZY/JLUGItfkhpj8UtSY/4fGlivxiQOHnkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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(r\"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 using the number of neighbors in a specified cutoff distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_neighbors = final_struct.get_numbers_of_neighbors_in_sphere(cutoff_radius=9).max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "neighbors = final_struct.get_neighbors(num_neighbors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "neigh_indices = np.hstack(np.array(neighbors.indices)[O_indices])\n",
    "neigh_distances = np.hstack(np.array(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": 18,
   "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": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEcCAYAAADZQfNOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAXGElEQVR4nO3df7RlZX3f8fdHBEFAYWTAQcDRBq1oItrRKJhGIWSBgqBVojGWuGxGkzTV+CMlUlPNagx22Wq0RkuNOsYfgVQJiD8ijqBNVHQgCFIwIAEljMwAKqBWBL79Y++rh8u9z5yZe/c5Z2ber7VmnbOfvffZ37t5uJ/77H323qkqJElazP2mXYAkabYZFJKkJoNCktRkUEiSmgwKSVKTQSFJajIotNNLckeSR067ji1JsjpJJbn/Nq7/uiTvWe66tOMzKDSYJL+Z5PIkP0zynSTvSrLPGOsdlORDSW5J8oMkX0ly/FB1VtVeVXXtUJ8/DUmenuSG0baqelNV/btp1aTtl0GhQSR5NfBm4LXAg4GnAA8Hzk+yW2O9FcDfAXcCjwX2A94KfDjJ84aue4F6dhnoc+8zKtjWkYI0NINCyy7Jg4A3Ar9XVZ+uqp9U1XXAyXRh8RuN1X8fuAN4aVV9p6p+VFUfAf4E+G9JssD25g7JrE1yY5KNfVDNzX9yki8l+V4/73+MhlW/7s/179/fj3w+meQHwDMW2N6KJO/rt/XdJH8zMu+3klyT5NYk5yY5cN52fjfJ1cDVc3/1J/mPSb4DvC/J/ZKcmuSb/YjqrD48F9rPL0lyZZLbk1yb5GV9+57Ap4AD+8NqdyQ5MMkbknxwZP1nJ7mi3y8XJnnMyLzrkrwmyWVJvp/kzCS7N/67aQdmUGgIRwC7Ax8bbayqO+h+gR3TWPcY4KNVdc+89rOAQ4BHNdZ9BnAo8KvAqUl+pW+/my6A9gOeChwN/E7jc36dLpj2phvdzPeXwAPpRjz70414SHIU8Kd0gbgKuB74q3nrngT8InBYP/1QYAVdgK4F/kO/zC8DBwLfBd65SJ2bgOOBBwEvAd6a5IlV9QPgOODG/rDaXlV14+iKSR4FfAR4JbAS+CTw8XmjvZOBY4FHAL8A/OYidWgHZ1BoCPsBN1fVXQvM29jPb627cZH15uYv5o1V9YOquhx4H/BCgKq6uKq+XFV39SOb/0n3i3gx51TV31fVPVX1/0ZnJFlF90v45VX13X609Pl+9ouA91bVJVX1Y+APgacmWT3yEX9aVbdW1Y/66XuA/1xVP+7bXgacVlU39J/xBuB5Cx2WqqpPVNU3q/N54DPALzV+rlG/Bnyiqs6vqp8AbwH2oAv5OW+vqhur6lbg48DhY362djAGhYZwM7DfIsfcV/XzSfLukUMjrxtZd9Ui683NX8y3R95fT/cXOUkeleS8/oT6bcCbaAfOtxvzDgZurarvLjDvwH67wE9HULcAD2t89uZ5YfRw4Oz+cND3gCvpRkQHzN9YkuOSfLk/zPU94Jm0f65Wrff0tY3W+p2R9z8E9hrzs7WDMSg0hC8BPwaeO9rYHzs/DlgPUFUvHzk08qZ+sc8C/ybJ/L55Mt0vsn9sbPfgkfeHAHOHW94FXAUcWlUPAl4H3Odcx4jWLZW/DaxY5NtbN9L9ogd++vM+BPjnxmfPn/42cFxV7TPyb/eqGv0MkjwA+CjdSOCAqtqH7vDR3M+1pdtCz681dPvvnxddQzstg0LLrqq+T3cy+x1Jjk2ya3/45a+BG+iO8S/mrXTH3P8iyUOT7J7khcBpwGurfV/81yd5YJLH0h2zP7Nv3xu4Dbgjyb8EfnsJP9tGuvMsf55k3/5n+9f97A8DL0lyeP+L/E3ARf3hrnG9G/iTJA8HSLIyyYkLLLcb8ABgM3BXkuPozs3MuQl4SJIHL7Kds4BnJTk6ya7Aq+nC/YtbUat2EgaFBlFV/5XuL/e30P2Svojur+Wj+2Pvi613C/A0upPh/5fu0M2rgBdX1ZmLrdf7PHAN3YjlLVX1mb79NXQnqG8H/hc/C5Bt9WLgJ3SjlE10J4SpqvXA6+n+0t8I/AvgBVv52X8GnAt8JsntwJfpTn7fS1XdTnfi+yy6E96/3q83N/8qupPV1/aHsQ6ct/436L599g66w3knACdU1Z1bWa92AvHBRdre9aOVfwJ2XeQEuqQlcEQhSWoyKCRJTR56kiQ1OaKQJDUZFJKkpu3ibpX77bdfrV69etplSNJ25eKLL765qlYu9XO2i6BYvXo1GzZsmHYZkrRdSXL9lpfaMg89SZKaDApJUpNBIUlqMigkSU0GhSSpadBvPSW5ju6OnXcDd1XVmv75v2cCq4HrgJMXeQiMJGkGTGJE8YyqOryq1vTTpwLrq+pQuttBnzqBGiRJ22gah55OBNb179fRPUhekjSjhg6KonsAy8VJ1vZtB/RPCZt7Wtj+C62YZG2SDUk2bN68eeAyJUmLGfrK7COr6sYk+wPnJ7lq3BWr6gzgDIA1a9Z4i1tJmpJBRxRVdWP/ugk4G3gycFOSVQD966Yha5AkLc1gQZFkzyR7z72ne/D71+me63tKv9gpwDlD1SBJWrohDz0dAJydZG47H66qTyf5KnBWkpcC3wKeP2ANkqQlGiwoqupa4PELtN8CHD3UdiVJy8srsyVJTQaFJKnJoJAkNRkUkqQmg0KS1GRQSJKaDApJUpNBIUlqMigkSU0GhSSpyaCQJDUZFJKkJoNCktRkUEiSmgwKSVKTQSFJajIoJElNBoUkqcmgkCQ1GRSSpCaDQpLUZFBIkpoMCklSk0EhSWoyKCRJTQaFJKnJoJAkNRkUkqQmg0KS1GRQSJKaDApJUpNBIUlqGjwokuyS5B+SnNdPr0hyfpKr+9d9h65BkrTtJjGieAVw5cj0qcD6qjoUWN9PS5Jm1KBBkeQg4FnAe0aaTwTW9e/XAScNWYMkaWmGHlG8DfgD4J6RtgOqaiNA/7r/QismWZtkQ5INmzdvHrhMSdJiBguKJMcDm6rq4m1Zv6rOqKo1VbVm5cqVy1ydJGlc9x/ws48Enp3kmcDuwIOSfBC4KcmqqtqYZBWwacAaJElLNNiIoqr+sKoOqqrVwAuAz1XVbwDnAqf0i50CnDNUDZKkpZvGdRSnA8ckuRo4pp+WJM2oIQ89/VRVXQhc2L+/BTh6EtuVJC2dV2ZLkpoMCklSk0EhSWoyKCRJTQaFJKnJoJAkNRkUkqQmg0KS1GRQSJKaDApJUpNBIUlqMigkSU0GhSSpyaCQJDUZFJKkJoNCktRkUEiSmgwKSVKTQSFJajIoJElNBoUkqcmgkCQ1GRSSpCaDQpLUZFBIkpoMCklSk0EhSWoyKCRJTQaFJKnJoJAkNRkUkqQmg0KS1DRYUCTZPclXknwtyRVJ3ti3r0hyfpKr+9d9h6pBkrR0Q44ofgwcVVWPBw4Hjk3yFOBUYH1VHQqs76clSTNqsKCozh395K79vwJOBNb17euAk4aqQZK0dIOeo0iyS5JLgU3A+VV1EXBAVW0E6F/3H7IGSdLSDBoUVXV3VR0OHAQ8Ocnjxl03ydokG5Js2Lx582A1SpLaJvKtp6r6HnAhcCxwU5JVAP3rpkXWOaOq1lTVmpUrV06iTEnSAob81tPKJPv07/cAfgW4CjgXOKVf7BTgnKFqkCQt3VhBkeTIcdrmWQVckOQy4Kt05yjOA04HjklyNXBMPy1JmlH3H3O5dwBPHKPtp6rqMuAJC7TfAhw9boGSpOlqBkWSpwJHACuTvGpk1oOAXYYsTJI0G7Y0otgN2Ktfbu+R9tuA5w1VlCRpdjSDoqo+D3w+yfur6voJ1SRJmiHjnqN4QJIzgNWj61TVUUMUJUmaHeMGxV8D7wbeA9w9XDmSpFkzblDcVVXvGrQSSdJMGveCu48n+Z0kq/rbhK9IsmLQyiRJM2HcEcXcldSvHWkr4JHLW44kadaMFRRV9YihC5EkzaaxgiLJv12ovao+sLzlSJJmzbiHnp408n53ultwXAIYFJK0gxv30NPvjU4neTDwl4NUJEmaKdt6m/EfAocuZyGSpNk07jmKj9N9ywm6mwE+BjhrqKIkSbNj3HMUbxl5fxdwfVXdMEA9kqQZM9ahp/7mgFfR3UF2X+DOIYuSJM2OcZ9wdzLwFeD5wMnARUm8zbgk7QTGPfR0GvCkqtoE3fOwgc8C/3uowiRJs2Hcbz3dby4kerdsxbqSpO3YuCOKTyf5W+Aj/fSvAZ8cpiRJ0izZ0jOzfw44oKpem+S5wNOAAF8CPjSB+iRJU7alw0dvA24HqKqPVdWrqur36UYTbxu2NEnSLNhSUKyuqsvmN1bVBrrHokqSdnBbCordG/P2WM5CJEmzaUtB8dUkvzW/MclLgYuHKUmSNEu29K2nVwJnJ3kRPwuGNcBuwHMGrEuSNCOaQVFVNwFHJHkG8Li++RNV9bnBK5MkzYRxn0dxAXDBwLVIkmaQV1dLkpoMCklSk0EhSWoyKCRJTQaFJKlpsKBIcnCSC5JcmeSKJK/o21ckOT/J1f3rvkPVIElauiFHFHcBr66qxwBPAX43yWHAqcD6qjoUWN9PS5Jm1GBBUVUbq+qS/v3twJXAw4ATgXX9YuuAk4aqQZK0dBM5R5FkNfAE4CK651tshC5MgP0XWWdtkg1JNmzevHkSZUqSFjB4UCTZC/go8Mqqum3c9arqjKpaU1VrVq5cOVyBkqSmQYMiya50IfGhqvpY33xTklX9/FXApsXWlyRN35DfegrwF8CVVfXfR2adC5zSvz8FOGeoGiRJSzfWTQG30ZHAi4HLk1zat70OOB04q3+mxbeA5w9YgyRpiQYLiqr6OyCLzD56qO1KkpaXV2ZLkpoMCklSk0EhSWoyKCRJTQaFJKnJoJAkNRkUkqQmg0KS1GRQSJKaDApJUpNBIUlqGvKmgNpJrT71E2Mtd93pz5rK50naOo4oJElNBoUkqcmgkCQ1GRSSpCaDQpLUZFBIkpoMCklSk9dRaGrGvT5CmoZx+ufOcu2OIwpJUpNBIUlqMigkSU0GhSSpyaCQJDUZFJKkJoNCktTkdRTaYfjcCk3azvLsFUcUkqQmg0KS1GRQSJKaDApJUtNgQZHkvUk2Jfn6SNuKJOcnubp/3Xeo7UuSlseQI4r3A8fOazsVWF9VhwLr+2lJ0gwbLCiq6gvArfOaTwTW9e/XAScNtX1J0vKY9DmKA6pqI0D/uv+Ety9J2koze8FdkrXAWoBDDjlkytVISzOrF1LNWc6HSO0IP4PubdIjipuSrALoXzcttmBVnVFVa6pqzcqVKydWoCTp3iYdFOcCp/TvTwHOmfD2JUlbacivx34E+BLw6CQ3JHkpcDpwTJKrgWP6aUnSDBvsHEVVvXCRWUcPtU1J0vLzymxJUpNBIUlqMigkSU0zex2FNJQd4fv2s35dhu5tufvcpPuwIwpJUpNBIUlqMigkSU0GhSSpyaCQJDUZFJKkJoNCktTkdRQa245w/YFmh9eCbD8cUUiSmgwKSVKTQSFJajIoJElNBoUkqcmgkCQ1GRSSpCavo5BmyI5wrcr2/uwF3ZcjCklSk0EhSWoyKCRJTQaFJKnJoJAkNRkUkqQmg0KS1GRQSJKavOBuB+fDYXZuXqym5eCIQpLUZFBIkpoMCklSk0EhSWqaSlAkOTbJN5Jck+TUadQgSRrPxIMiyS7AO4HjgMOAFyY5bNJ1SJLGM40RxZOBa6rq2qq6E/gr4MQp1CFJGsM0rqN4GPDtkekbgF+cv1CStcDafvLHSb4+gdqWaj/g5mkXMYb71Jk3T6mStu1hf24PNYJ1Lrftpc5HL8eHTCMoskBb3aeh6gzgDIAkG6pqzdCFLZV1Lq/toc7toUawzuW2PdW5HJ8zjUNPNwAHj0wfBNw4hTokSWOYRlB8FTg0ySOS7Aa8ADh3CnVIksYw8UNPVXVXkn8P/C2wC/DeqrpiC6udMXxly8I6l9f2UOf2UCNY53LbqepM1X1OD0iS9FNemS1JajIoJElNUw2KJAcnuSDJlUmuSPKKBZZJkrf3t/u4LMkTR+YNfiuQMWt8UV/bZUm+mOTxI/OuS3J5kkuX66tqS6jz6Um+39dyaZI/Gpk3kduqjFnna0dq/HqSu5Os6OdNan/unuQrSb7W1/nGBZaZat/cijpnoX+OU+dU++eYNU69b47UskuSf0hy3gLzlrdvVtXU/gGrgCf27/cG/hE4bN4yzwQ+RXf9xVOAi/r2XYBvAo8EdgO+Nn/dCdZ4BLBv//64uRr76euA/WZkXz4dOG+BdSeyL8etc97yJwCfm8L+DLBX/35X4CLgKbPUN7eizlnon+PUOdX+OU6Ns9A3R7b3KuDDi+yzZe2bUx1RVNXGqrqkf387cCXdldujTgQ+UJ0vA/skWcWEbgUyTo1V9cWq+m4/+WW6a0Mmasx9uZiJ3VZlG+p8IfCRIWpp6fvbHf3krv2/+d/8mGrfHLfOGemf4+zPxUzq//WtrXEqfRMgyUHAs4D3LLLIsvbNmTlHkWQ18AS6FB+10C0/HtZoH0yjxlEvpUvyOQV8JsnF6W5LMrgt1PnUfmj9qSSP7dsmvi9hy/szyQOBY4GPjjRPbH/2Q/tLgU3A+VU1k31zjDpHTa1/jlnnVPvnuPty2n0TeBvwB8A9i8xf1r45E8/MTrIX3Q5/ZVXdNn/2AqtUo30QW6hxbpln0P2P+LSR5iOr6sYk+wPnJ7mqqr4wpTovAR5eVXckeSbwN8ChTHhfjlHnnBOAv6+qW0faJrY/q+pu4PAk+wBnJ3lcVY3ec2wm+uYYdQLT759j1Dn1/jnuvmSKfTPJ8cCmqro4ydMXW2yBtm3um1MfUSTZle4Xxoeq6mMLLLLYLT8mdiuQMWokyS/QDQNPrKpb5tqr6sb+dRNwNt3QbxBbqrOqbpsbWlfVJ4Fdk+zHhG+rMs7+7L2AeUP7Se7PkW1+D7iQ7i/IUVPvm6Madc5E/9xSnbPSP1s1jphm3zwSeHaS6+gOHR2V5IPzllnevrmlkxhD/qNLtw8Ab2ss8yzufVLmK337/YFrgUfws5Myj51SjYcA1wBHzGvfE9h75P0XgWOnuC8fys8usnwy8K1+vYnsy3Hr7Jd7MHArsOeU9udKYJ/+/R7A/wGOn6W+uRV1zkL/HKfOqfbPcWqchb45r5ans/DJ7GXtm9M+9HQk8GLg8v64IMDr6Do2VfVu4JN0Z/CvAX4IvKSfty23Ahmqxj8CHgL8eRKAu6q7s+QBdMNX6P4DfbiqPj1AjePW+Tzgt5PcBfwIeEF1vWdS+3LcOgGeA3ymqn4wsu4k9+cqYF26B23dDzirqs5L8vKROqfdN8etcxb65zh1Trt/jlMjTL9vLmjIvuktPCRJTVM/RyFJmm0GhSSpyaCQJDUZFJKkJoNCktRkUEiSmgwKaQKSvCPJJUmeNO1apK1lUEgDS7InsD/wMuD4KZcjbTWDQuoleUOS1yxh/dVJfjRyxTkA/RW8q+juHfT2ftk9+gfc3Nnfz0iaWQaFtLy+WVWHjzYkeQjwQOB24G6AqvpRv9zgNwuUlsqg0E4tyWn9YyE/Czx6oM38J+AtwBXAYQNtQxqMQaGdVpJ/RXe76CcAzwWW/URz/3CmI4Az6Z7m99jmCtIMmvbdY6Vp+iXg7Kr6IUCSc+dmJHk+3QN+dgG+X1WnLdQ2xjb+C/DHVVVJDAptlwwK7ezuc/vkJEcCa6rqFf30u5P88gJtj66qbyz2wUkOpxupPC3JO4HdgcsH+BmkQXnoSTuzLwDP6b+BtDfd4y2he1zoO+Yt+4YF2u7cwue/GTihqlZX1Wrg8Tii0HbIEYV2WlV1SZIzgUuB6+meaAawK/1II8kj6K6BuGF+W1X902KfneQouiegrR/Z3k1J9kyyou79rGVppvngImmeJD8PnAZsoguN19NdB3Gvtqq6ed56q+keS/m4rdjWdXSHtG7e0rLStBgU0jJJcjDds5JvmX8txQLL7gF8ie45zT/vCEOzzKCQJDV5MluS1GRQSJKaDApJUpNBIUlqMigkSU0GhSSpyaCQJDUZFJKkpv8PEVjgKsZCzi0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "bins = np.linspace(1, 8, 120)\n",
    "count = plt.hist(O_neigh_distances, bins=bins)\n",
    "plt.xlim(2, 4)\n",
    "plt.title(\"O-O pair correlation\")\n",
    "plt.xlabel(r\"d$_{OO}$ [$\\AA$]\")\n",
    "plt.ylabel(\"Count\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now extent our analysis to include statistically independent snapshots along the trajectory. This allows to obtain the radial pair distribution function of O-O distances in the NVT ensamble."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "traj=ham_eq[\"output/generic/positions\"]\n",
    "nsteps=len(traj)\n",
    "stepincrement=int(nsteps/10)\n",
    "# Start sampling snaphots after some equilibration time; do not double count last step.\n",
    "snapshots=range(stepincrement,nsteps-stepincrement,stepincrement)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in snapshots:\n",
    "    struct.positions = traj[i]\n",
    "    neighbors = struct.get_neighbors(num_neighbors)\n",
    "    neigh_indices = np.hstack(np.array(neighbors.indices)[O_indices])\n",
    "    neigh_distances = np.hstack(np.array(neighbors.distances)[O_indices])\n",
    "    O_neigh_indices  = np.in1d(neigh_indices, O_indices) \n",
    "    #collect all distances in the same array\n",
    "    O_neigh_distances = np.concatenate((O_neigh_distances,neigh_distances[O_neigh_indices]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To obtain a radial pair distribution function ($g(r)$), one has to normalize by the volume of the surfce increment of the sphere ($4\\pi r^2\\Delta r$) and by the number of species, samples, and the number density."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "O_gr=np.histogram(O_neigh_distances,bins=bins)\n",
    "dr=O_gr[1][1]-O_gr[1][0]\n",
    "normfac=(n/a)**3*n**3*4*np.pi*dr*(len(snapshots)+1)\n",
    "# (n/a)**3            number density\n",
    "# n**3                number of species\n",
    "# (len(snapshots)+1)  number of samples (we also use final_struct)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEcCAYAAAAydkhNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYR0lEQVR4nO3de7hddX3n8fdHDIKAppioARKOnUFnRC0wKUhRi6h9CKJ0HMeBVh2ZzkSoOjrVziCO9X5pH0dHRUlT7xVBFLVUYitab0wFDWkAMVgj4kNMkAByCVAQ+M4fe8Vsds46l5yctc/l/Xqe/WTvtX5r7e/ZJ2d/9u/3W2vtVBWSJI3mIcMuQJI0cxkSkqRWhoQkqZUhIUlqZUhIkloZEpKkVoaENIYk25L85rDrGE+SkSSV5KG7uP2ZST6yu+vS7GdIaCiSvCzJVUnuSnJDkrOTLJzAdgclOSfJzUnuTPK9JCdOV51VtW9VXTtd+x+GJMcm2dS/rKreWVX/dVg1aeYyJNS5JK8F/hz4U+CRwFOBg4GLk+w5xnb7A5cA9wKHAouA9wGfSfLC6a57lHr2mKb97tQb2NUegjRVhoQ6leQRwFuAV1XV31XVr6rqOuBF9ILixWNs/j+AbcAfVdUNVXV3VZ0LvAP4P0kyyvNtH4ZZmWRzki1NSG1ff2SS7ya5tVl3Vn9QNdv+6+b+J5oez5okdwLPHOX59k/y8ea5fpnkS33r/luSjUluSXJhkgMGnucVSX4M/Hj7p/0k/yvJDcDHkzwkyRlJftL0pM5vgnO01/nUJBuS3JHk2iQvb5bvA3wFOKAZStuW5IAkb07y6b7tn5/k6uZ1+WaSf9u37rokr0tyZZLbknw2yV5j/N40ixkS6trvAHsBX+hfWFXb6L15PWeMbZ8DXFBVDwwsPx9YBjx+jG2fCRwC/B5wRpJnN8vvpxc+i4CjgWcBfzzGfv6AXijtR69XM+ivgYfT6+k8ml5PhyTHAe+iF4ZLgJ8B5w1s+/vAUcATm8ePBfanF54rgf/etPld4ADgl8CHWuq8ETgReARwKvC+JEdU1Z3ACmBzM5S2b1Vt7t8wyeOBc4HXAIuBNcDfDvTyXgQcDzwOeArwspY6NMsZEuraIuCmqrpvlHVbmvVjbbulZbvt69u8parurKqrgI8DpwBU1eVVdWlV3df0aP6S3ptwm7+pqv9XVQ9U1b/0r0iyhN4b8GlV9cuml/StZvUfAh+rqnVVdQ/weuDoJCN9u3hXVd1SVXc3jx8A3lRV9zTLXg68oao2Nft4M/DC0YaiquqiqvpJ9XwL+Crw9DF+rn7/Cbioqi6uql8B7wH2phfw232gqjZX1S3A3wKHTXDfmmUMCXXtJmBRyxj7kmY9SVb1DYec2bftkpbttq9vc33f/Z/R+yROkscn+XIzeX478E7GDpvrx1i3FLilqn45yroDmucFft1zuhk4cIx9bx0IooOBLzZDQLcCG+j1hB4z+GRJViS5tBnauhU4gbF/rrFqfaCprb/WG/ru3wXsO8F9a5YxJNS17wL3AC/oX9iMla8Avg5QVaf1DYe8s2n2NeA/JBn8f/siem9i/zzG8y7tu78M2D7EcjZwDXBIVT0COBPYaW6jz1iXTb4e2L/lKK3N9N7kgV//vI8Cfj7GvgcfXw+sqKqFfbe9qqp/HyR5GHABvR7AY6pqIb0ho+0/13iXfh6sNfRev5+3bqE5y5BQp6rqNnoT1x9McnySBc2Qy+eATfTG9Nu8j94Y+0eTPDbJXklOAd4A/GmNfd37NyZ5eJJD6Y3Rf7ZZvh9wO7Atyb8BTp/Cz7aF3rzKh5P8RvOzPaNZ/Rng1CSHNW/i7wQua4a4JmoV8I4kBwMkWZzkpFHa7Qk8DNgK3JdkBb25mO1+ATwqySNbnud84LlJnpVkAfBaesH+j5OoVXOEIaHOVdVf0PvE/h56b9CX0fuU/KxmrL1tu5uBp9Gb+P4hveGaPwFeUlWfbduu8S1gI72eynuq6qvN8tfRm4y+A/grdoTHrnoJ8Ct6vZMb6U3+UlVfB95I7xP+FuBfASdPct/vBy4EvprkDuBSehPdD1JVd9Cb5D6f3uT2HzTbbV9/Db2J6WuboasDBrb/Eb2jzD5IbwjvecDzqureSdarOSB+6ZDmsqaX8lNgQctkuaQx2JOQJLUyJCRJrRxukiS1sichSWplSEiSWs36K0suWrSoRkZGhl2GJM0ql19++U1VtXi8drM+JEZGRli7du2wy5CkWSXJz8Zv5XCTJGkMhoQkqZUhIUlqZUhIkloZEpKkVoaEJKmVISFJamVISJJazfqT6WajkTMuetDj69793CFVIkljsychSWplSEiSWhkSkqRWhoQkqZUhIUlqZUhIkloZEpKkVoaEJKmVISFJamVISJJaGRKSpFaGhCSplSEhSWplSEiSWhkSkqRWnYVEkr2SfC/JFUmuTvKWUdokyQeSbExyZZIjuqpPkrSzLr906B7guKralmQBcEmSr1TVpX1tVgCHNLejgLObfyVJQ9BZT6J6tjUPFzS3Gmh2EvCppu2lwMIkS7qqUZL0YJ3OSSTZI8l64Ebg4qq6bKDJgcD1fY83NcskSUPQaUhU1f1VdRhwEHBkkicNNMlomw0uSLIyydoka7du3ToNlUqSYEhHN1XVrcA3geMHVm0ClvY9PgjYPMr2q6tqeVUtX7x48XSVKUnzXpdHNy1OsrC5vzfwbOCagWYXAi9tjnJ6KnBbVW3pqkZJ0oN1eXTTEuCTSfagF07nV9WXk5wGUFWrgDXACcBG4C7g1A7rkyQN6CwkqupK4PBRlq/qu1/AK7qqSZI0Ns+4liS1MiQkSa0MCUlSK0NCktTKkJAktTIkJEmtujxPYl4ZOeOiYZcgSVNmT0KS1MqQkCS1MiQkSa0MCUlSK0NCktTKkJAktTIkJEmtDAlJUitDQpLUypCQJLUyJCRJrQwJSVIrQ0KS1MqQkCS1MiQkSa0MCUlSq85CIsnSJN9IsiHJ1UlePUqbY5PclmR9c/uzruqTJO2sy2+muw94bVWtS7IfcHmSi6vqhwPtvlNVJ3ZYlySpRWc9iaraUlXrmvt3ABuAA7t6fknS5A1lTiLJCHA4cNkoq49OckWSryQ5tGX7lUnWJlm7devW6SxVkua1zkMiyb7ABcBrqur2gdXrgIOr6reADwJfGm0fVbW6qpZX1fLFixdPa72SNJ91GhJJFtALiHOq6guD66vq9qra1txfAyxIsqjLGiVJO3R5dFOAjwIbquq9LW0e27QjyZFNfTd3VaMk6cG6PLrpGOAlwFVJ1jfLzgSWAVTVKuCFwOlJ7gPuBk6uquqwRklSn85CoqouATJOm7OAs7qpSJI0Hs+4liS1MiQkSa0MCUlSK0NCktTKkJAktTIkJEmtDAlJUitDQpLUqsszrtVi5IyLHvT4unc/d0iVSNKD2ZOQJLUyJCRJrQwJSVIrQ0KS1MqQkCS1MiQkSa0MCUlSK0NCktTKkJAktTIkJEmtDAlJUitDQpLUypCQJLXqLCSSLE3yjSQbklyd5NWjtEmSDyTZmOTKJEd0VZ8kaWddXir8PuC1VbUuyX7A5Ukurqof9rVZARzS3I4Czm7+lSQNQWc9iaraUlXrmvt3ABuAAweanQR8qnouBRYmWdJVjZKkBxvKnESSEeBw4LKBVQcC1/c93sTOQUKSlUnWJlm7devWaatTkua7zkMiyb7ABcBrqur2wdWjbFI7LahaXVXLq2r54sWLp6NMSRIdh0SSBfQC4pyq+sIoTTYBS/seHwRs7qI2SdLOJh0SSfZJsscubBfgo8CGqnpvS7MLgZc2Rzk9FbitqrZM9rkkSbvHuEc3JXkIcDLwh8BvA/cAD0uyFVgDrK6qH0/guY4BXgJclWR9s+xMYBlAVa1q9ncCsBG4Czh1Mj+MJGn3msghsN8Avga8HvhBVT0AkGR/4JnAu5N8sao+PdZOquoSRp9z6G9TwCsmUrgkafpNJCSeDdwPnFFVV25fWFW30JtfuKCZa5AkzTHjzklU1a+a3sOzx2qzW6uSJM0Ik5m4/qckb2rmKCRJ88BkLsuxFHgycHqSy4ArgSur6nPTUpkkaegmHBJV9SKAJA8DDqUXGEcBhoQkzVETOQQ2zVFHAFTVPcC65jZqG0nS3DCR+YVvJHlVkmX9C5PsmeS4JJ8E/vP0lCdJGqaJDDcdD/wX4NwkjwNuBfamFzBfBd5XVeunq0BJ0vCMGxJV9S/Ah4EPN+dDLALurqpbp7m2WWXkjIuGXYIk7XaTOpy1OWdiC7A4ycHTVJMkaYaY9DfTJXk78Kjm/mJgZXP2tSRpjtmVry99ZFWdDtDMUbwNr7ckSXPSrpw9ff/2O1X1U8BLckjSHLUrIXFJknclOSDJAcBjdndRkqSZYdLDTVX1+SQ3An9N75vk3rHbq5IkzQi7MnH9NnqHwf4z8Gj8elFJmrN2ZeJ6Yd/E9QhOXEvSnDXVievrcOJakuYsJ64lSa12deJ6M71hpofixLUkzVkTDokkPwauAq4A1gNva4abJElz1GSGm/4SuAG4GVgBXJXkqiRvbS78J0maYyYTEi+uqj+uqrOq6jTg6cA/ALcD7x1v4yQfS3Jjkh+0rD82yW1J1je3P5tEbZKkaTCZkLgtyVO2P2i+Q+KpVfUe4JgJbP8Jet9NMZbvVNVhze2tk6hNkjQNJjNx/XLgnCTr6c1JPAF4oFm353gbV9W3m/MqJEmzxIRDoqquSXIk8ALgKcBG4E1J9gHO2031HJ3kCnpncb+uqq4erVGSlcBKgGXLlo3WZFYb/AKj69793CFVImm+m9QhsFV1P/C55tbv7buhlnXAwVW1LckJwJeAQ1rqWA2sBli+fHnthueWJI1iVy7LMS2q6va++2uSfDjJoqq6aZh1zQT2LCQNy66ccT0tkjw2SZr7R9Kr7ebhViVJ81tnPYkk5wLHAouSbALeBCwAqKpVwAuB05PcB9wNnFxVDiVJ0hB1FhJVdco4688CzuqoHEnSBMyY4SZJ0sxjSEiSWhkSkqRWhoQkqZUhIUlqZUhIkloZEpKkVoaEJKmVISFJamVISJJaGRKSpFYz5lLhmruGealzL7MuTY0hIU0jA1KznSEhSfPI4IeH8RgS0gwx3h+vPQENgxPXkqRWhoQkqZXDTbOQE5Ld8bXWfGdISNIsNt0fZAwJzSmTPXKja/31dd0rsVekXWFIzEG+GUgz22Q+LAz779mQkKQpGvYb+XTqLCSSfAw4Ebixqp40yvoA7wdOAO4CXlZV67qqT9oVs/mM6vG2n8tvfJq4LnsSnwDOAj7Vsn4FcEhzOwo4u/lX2m0m+8ao6THZAJptgTWd9Xb9f7SzkKiqbycZGaPJScCnqqqAS5MsTLKkqrZ0U6GGYbb98e9Ok/1jN8Bmjvn0u5hJcxIHAtf3Pd7ULDMkpAmYT29c45nqUNru/vAy1u9mup9rqvubSSGRUZbVqA2TlcBKgGXLlk1nTdK8YcjMDDPt9zCTQmITsLTv8UHA5tEaVtVqYDXA8uXLRw0SaRhm2h+4esb7vUx1/Vw2k0LiQuCVSc6jN2F9m/MRmm7z+Y9fmoguD4E9FzgWWJRkE/AmYAFAVa0C1tA7/HUjvUNgT+2qNk1Nl0dyzKeJ7ZlmqpcyN5Bnpy6PbjplnPUFvKKjcuaU+fTHZ2hoNPPpb6BrM2m4SZrxfDPSbDPV/7OGhKTdYncGqGE8cxgSmtV8M5Gml99MJ0lqZU9CnfPTvzR7GBLzkEcIaabxg8PM5XCTJKmVPYl5YKqf0ux5aKax59EdexKSpFaGhCSplcNNcjhJUit7EpKkVvYkNKM4ISnNLPYkJEmt7Ensorn8iXcu/2ySJseQ0KQZItL84XCTJKmVISFJamVISJJaGRKSpFaGhCSplSEhSWrV6SGwSY4H3g/sAXykqt49sP5Y4G+AnzaLvlBVb+2yRk2dh8hKc0dnIZFkD+BDwHOATcD3k1xYVT8caPqdqjqxq7okSe26HG46EthYVddW1b3AecBJHT6/JGmSugyJA4Hr+x5vapYNOjrJFUm+kuTQbkqTJI2myzmJjLKsBh6vAw6uqm1JTgC+BByy046SlcBKgGXLlu3mMiVJ23XZk9gELO17fBCwub9BVd1eVdua+2uABUkWDe6oqlZX1fKqWr548eLprFmS5rUuQ+L7wCFJHpdkT+Bk4ML+BkkemyTN/SOb+m7usEZJUp/Ohpuq6r4krwT+nt4hsB+rqquTnNasXwW8EDg9yX3A3cDJVTU4JCVJ6kin50k0Q0hrBpat6rt/FnBWlzVJktp5xrUkqZUhIUlqZUhIkloZEpKkVoaEJKmVISFJamVISJJaGRKSpFaGhCSplSEhSWrV6WU5ZjO/klPSfGRPQpLUypCQJLUyJCRJrQwJSVIrQ0KS1MqQkCS1MiQkSa0MCUlSK0NCktTKkJAktTIkJEmtDAlJUqtOQyLJ8Ul+lGRjkjNGWZ8kH2jWX5nkiC7rkyQ9WGchkWQP4EPACuCJwClJnjjQbAVwSHNbCZzdVX2SpJ112ZM4EthYVddW1b3AecBJA21OAj5VPZcCC5Ms6bBGSVKfLr9P4kDg+r7Hm4CjJtDmQGBLf6MkK+n1NADuSfKD3VvqrLUIuGnYRcwQvhY7+Frs4GuxwxMm0qjLkMgoy2oX2lBVq4HVAEnWVtXyqZc3+/la7OBrsYOvxQ6+FjskWTuRdl0ON20ClvY9PgjYvAttJEkd6TIkvg8ckuRxSfYETgYuHGhzIfDS5iinpwK3VdWWwR1JkrrR2XBTVd2X5JXA3wN7AB+rqquTnNasXwWsAU4ANgJ3AadOYNerp6nk2cjXYgdfix18LXbwtdhhQq9FqnYa8pckCfCMa0nSGAwJSVKrWRsSSZYm+UaSDUmuTvLqYdc0LEn2SvK9JFc0r8Vbhl3TMCXZI8k/JfnysGsZtiTXJbkqyfqJHvI4VyVZmOTzSa5p3jeOHnZNw5DkCc3/h+2325O8prX9bJ2TaM7EXlJV65LsB1wO/H5V/XDIpXUuSYB9qmpbkgXAJcCrm7PW550kfwIsBx5RVScOu55hSnIdsLyq5v0JZEk+CXynqj7SHGH58Kq6dchlDVVzuaSfA0dV1c9GazNrexJVtaWq1jX37wA20Ds7e95pLmOyrXm4oLnNzvSfoiQHAc8FPjLsWjRzJHkE8AzgowBVde98D4jGs4CftAUEzOKQ6JdkBDgcuGzIpQxNM8SyHrgRuLiq5utr8X+B/wk8MOQ6ZooCvprk8uZyNvPVbwJbgY83Q5EfSbLPsIuaAU4Gzh2rwawPiST7AhcAr6mq24ddz7BU1f1VdRi9s9SPTPKkIZfUuSQnAjdW1eXDrmUGOaaqjqB3heVXJHnGsAsakocCRwBnV9XhwJ3ATl9XMJ80Q27PBz43VrtZHRLN+PsFwDlV9YVh1zMTNF3obwLHD7eSoTgGeH4zDn8ecFySTw+3pOGqqs3NvzcCX6R3Neb5aBOwqa+H/Xl6oTGfrQDWVdUvxmo0a0Oimaz9KLChqt477HqGKcniJAub+3sDzwauGWpRQ1BVr6+qg6pqhF43+h+q6sVDLmtokuzTHNRBM7Tye8C8vGJyVd0AXJ9k+5VPnwXMu4NcBpzCOENN0O1VYHe3Y4CXAFc1Y/EAZ1bVmuGVNDRLgE82Ryo8BDi/qub94Z/iMcAXe5+neCjwmar6u+GWNFSvAs5phlmuZWKX/ZmTkjwceA7w8nHbztZDYCVJ02/WDjdJkqafISFJamVISJJaGRKSpFaGhCSplSEhSWplSEhDluSDSdYl+e1h1yINMiSkIWrOhH40vZOa5vVlzTUzGRLSBCR5c5LXTWH7kSR3910dAICqupPeGfPfBD7QtN27+TKYe5MsmkLZ0pQZElJ3ftJcqffXkjwKeDhwB3A/QFXd3bTb3HWB0iBDQmqR5A1JfpTka8ATxt1g1/xv4D3A1cATp+k5pF1mSEijSPLv6F1J9nDgBcBun1Ruvizrd4DP0vtmxUN393NIUzWbrwIrTaenA1+sqrsAkly4fUWS/wg8DdgDuK2q3jDasgk8x9uBt1ZVJTEkNCMZElK7nS6RnOQYYHlVvbp5vCrJ746y7AlV9aO2HSc5jF4P5WlJPgTsBVw1DT+DNCUON0mj+zbw75sjjfYDntcs/yPggwNt3zzKsnvH2f+fA8+rqpHmS5J+C3sSmoHsSUijqKp1ST4LrAd+BnynWbWApoeR5HH0znHYNLisqn7atu8kxwH7VNXX+57vF803ye1fVbdMw48k7RK/dEiahCRPBt4A3EgvMN5I7zyHBy2rqpsGthsBvlxVT5rEc11HbxjrpvHaStPFkJA6kGQp8I/AzYPnSozSdm/gu8Bi4Mn2LDRMhoQkqZUT15KkVoaEJKmVISFJamVISJJaGRKSpFaGhCSplSEhSWplSEiSWv1/2RdxTk6IE8EAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.bar(O_gr[1][0:-1],O_gr[0]/(normfac*O_gr[1][0:-1]**2),dr)\n",
    "plt.xlim(2, 7)\n",
    "plt.title(\"O-O pair correlation\")\n",
    "plt.xlabel(r\"d$_{OO}$ [$\\AA$]\")\n",
    "plt.ylabel(\"$g_{OO}(r)$\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  },
  "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
}