{ "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\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": "9fd59fa15e0b46258ce4366dc1731c2d", "version_major": 2, "version_minor": 0 }, "text/plain": [ "_ColormakerRegistry()" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "98eb3f129f784853af47319a07a5a6d2", "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)\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": "0a6e30ec113840c3beb0ba9a20a1eae2", "version_major": 2, "version_minor": 0 }, "text/plain": [ "NGLWidget(max_frame=100)" ] }, "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": "24602902924748579a0c4c78a56ef052", "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": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEGCAYAAACZ0MnKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO2dZ5gb1dWA37Pavq644YIxxRCDMW0B0wkGYnonQAgkhFBCKilAHBIIkIQSkg9ICA4JCTV0SOg2vRiMKS7YFHcvuHd711uk8/3QjHYkzUizu9Jqy3mfR89Kd2akO9rRnHu6qCqGYRiGEYaiQk/AMAzD6DyY0DAMwzBCY0LDMAzDCI0JDcMwDCM0JjQMwzCM0BQXegL5pn///jpixIhCT8MwDKNT8f77769S1QGp411eaIwYMYJp06YVehqGYRidChFZ5Ddu5inDMAwjNCY0DMMwjNCY0DAMwzBCY0LDMAzDCI0JDcMwDCM0JjQMwzCM0JjQMAzDMEJjQiOAf721gP9N/7LQ0zAMw+hQmNAI4IGpi3lmxtJCT8MwDKNDYUIjgPKSCFuaooWehmEYRofChEYA5cURtjSa0DAMw/BSEKEhIqeLyMciEhORap/tw0Vkk4j8zDO2t4jMFJG5InKriEg+51hWUsSWxlg+P8IwDKPTUShNYxZwCvB6wPY/Ac+ljN0BXAiMdB7j8zY7HPOUaRqGYRhJFERoqOocVf3Ub5uInATMBz72jA0GeqnqFFVV4B7gpHzOsbwkQn2TaRqGYRheOpRPQ0SqgMuBa1I2DQVqPK9rnLGg97lQRKaJyLSVK1e2ai7lxUWmaRiGYaSQN6EhIpNFZJbP48QMh10D/ElVN6W+nc++GvQmqjpRVatVtXrAgLQeIqEw85RhGEY6eWvCpKpHtOKw/YDTRORGoA8QE5EtwGPAMM9+w4C8Zt6VlxRRZ0LDMAwjiQ7VuU9VD3afi8jVwCZVvd15vVFExgLvAucCt+VzLiWRIpqigcqMYRhGt6RQIbcni0gNsD/wjIi8EOKwS4C7gLnAPNKjq3JKcZEQVRMahmEYXgqiaajqE8ATWfa5OuX1NGB0HqeVRFGRoAqxmFJUlNeUEMMwjE5Dh4qe6khEnNxB0zYMwzCaMaERQCTiCI2YCQ3DMAwXExoBJDQNExqGYRgJTGgEECky85RhGEYqJjQCcIVGg5USMQzDSGBCIwBXaFRfN7nAMzEMw+g4mNAIIGJhtoZhGGmY0Agg4mnXYc5wwzCMOCY0AvAm9G3c0ljAmRiGYXQcTGgEUOwRGuvrTGgYhmGACY1AIkmaRlMBZ2IYhtFxMKERQJHHp2F9NQzDMOKY0AjAa56ytq+GYRhxTGgE4HWEm6ZhGIYRx4RGAMVJQsM0DcMwDDChEYhpGoZhGOmY0AjAm9y3pcmEhmEYBhSu3evpIvKxiMREpDpl2xgRmeJsnyki5c743s7ruSJyq4jktc6H1zy1eHVtPj/KMAyj01AoTWMWcArwundQRIqB+4CLVXVX4DDAzay7A7gQGOk8xudzgg3RZj/Gna/Pz+dHGYZhdBoKIjRUdY6qfuqz6ShghqpOd/ZbrapRERkM9FLVKaqqwD3ASfmc44qN9fl8e8MwjE5JR/Np7ASoiLwgIh+IyC+c8aFAjWe/GmfMFxG5UESmici0lStXtmoix+w2mBN2H8LJew6lsjTSqvcwDMPoauRNaIjIZBGZ5fM4McNhxcBBwDecvyeLyDjAz38RWHpWVSeqarWqVg8YMKBV8+9RVsytZ+3JDgOqqG2IUm/OcMMwDIrz9caqekQrDqsBXlPVVQAi8iywF3E/xzDPfsOAL9s8yRD0qSwFYH1tIwN7mcZhGEb3pqOZp14AxohIpeMUPxSYrapLgY0iMtaJmjoXeKo9JtSjLC5XN9Vb0ULDMIxChdyeLCI1wP7AMyLyAoCqrgVuAd4DPgI+UNVnnMMuAe4C5gLzgOfaY65VjtDYXG/mKcMwjLyZpzKhqk8ATwRsu4+4OSp1fBowOs9TS6PKcYJvbjBNwzAMo6OZpzoczZqGCQ3DMAwTGlmoKnM1DTNPGYZhmNDIgmkahmEYzZjQyEJZcVzTsEq3hmEYJjSyUlYc/4oarHufYRiGCY1slJrQMAzDSGBCIwvFRYJIctVbwzCM7ooJjSyICGXFRdSbpmEYhmFCIwylkSIzTxmGYWBCIxSlxRHTNAzDMDChEYq4ecpCbg3DMAJrT4nIhizHCrBUVXfK7ZQ6HmXFZp4yDMOAzAUL56nqnpkOFpEPczyfDkmpCQ3DMAwgs3nq1BDHh9mn01NRGrEqt4ZhGGQWGpeJyIGZDlbV+TmeT4dk+FaVLFpdW+hpGIZhFJxMQuNz4GYRWSgiN4jIHu01qY7Gdv2r+GJdndWfMgyj2xMoNFT1/1R1f+ItV9cAd4vIHBH5tYh0eee3l+36V6EKi9eYtmEYRvcma8itqi5S1Rscp/jZwMnAnLZ8qIicLiIfi0hMRKo94yUi8m8RmekIqCs92/Z2xueKyK1Or/B2Ybv+VQDMX7m5vT7SMAyjQ5JVaDg38uNF5H7ifbk/o+0O8FnAKcDrKeOnA2WquhuwN3CRiIxwtt0BXAiMdB7j2ziH0IxwhMaCVSY0DMPo3mTK0zgSOAs4FpgK/Ae4UFXbfOdU1TnOZ6RtAqpEpBioABqADSIyGOilqlOc4+4BTiIuxPJOr/IS+vcoZdFqExqGYXRvMmkavwSmAKNU9XhVvT8XAiMLjwKbgaXAYuBmVV0DDAVqPPvVOGO+iMiFIjJNRKatXLkyJxPrV1XG2tqGnLyXYRhGZyVQ01DVr7rPReQgYKSq3i0iA4Aeqrog0xuLyGRga59NE1T1qYDD9gWiwBCgL/CG8z5+/gvNMPeJwESA6urqwP1aQq+KYjbUWa6GYRjdm0wZ4QCIyG+AamBn4G6gBLgPyJbDcUQr5nM28LyqNgIrROQt57PfAIZ59hsGfNmK9281vcpLWLZhS3t+pGEYRocjTMHCk4ETiJuNUNUvgZ55ms9i4HCJUwWMBT5R1aXARhEZ60RNnQsEaSt5oVdFCRu2NLbnRxqGYXQ4wgiNBlVVHHOQczNvEyJysojUAPsDz4jIC86mvwA9iEdXvQfcraoznG2XAHcBc4F5tJMT3KVnuZmnDMMwspqngIdF5E6gj4h8Fzgf+HtbPlRVnwCe8BnfRDzs1u+YacDotnxuWygrLqLRWr4ahtHNySo0VPVmJ/x2A3G/xq9VdVLeZ9bBKImY0DAMwwijaeAIiW4nKLzEhYYSiylFRe2WjG4YhtGhCPRpiMjT2Q4Os09XobQ4/lU1xkzbMAyj+5JJ0zhIRP6bYbsAu+R4Ph2Wkkhcu2iMKmWh9DPDMIyuR6bb34khju82KdKlEUfTaIpBWYEnYxiGUSAyZYS/1p4T6eiUuOYpc4YbhtGNCZOnYRB3hAM0mNAwDKMbY0IjJK55qqHJhIZhGN2XMP00jhORbi9cXE3DKt0ahtGdCSMMzgQ+F5EbRWRUvifUUXFbf3zn39MKOxHDMIwCEqbd6znAnsTrPd0tIlOcfhX5KlrYIflyXR0A62qtaKFhGN2XUGYnVd0APEa8e99g4pVvPxCRH+Rxbh2KE/eI93w6ZKcBBZ6JYRhG4Qjj0zheRJ4AXibeS2NfVT0a2B34WZ7n12EY0LOM7QdU0bPcMvsMw+i+hLkDng78SVVf9w6qaq2InJ+faXVMKkoi1DdGCz0NwzCMghGmyu25Gba9lNvpdGwqSiLUmdAwDKMbE8Y8tVFENqQ8lojIEyKyfXtMsqNQURqhrsGEhmEY3Zcw5qlbiPfjfoB4kcIzga2BT4F/Aofla3IdjfKSCKs2WZ6GYRjdlzDRU+NV9U5V3aiqG1R1InCMqj4E9M3z/DoUZcVF5tMwDKNbE0ZoxETkDBEpch5neLZpaz5URG4SkU9EZIZj5urj2XaliMwVkU9F5Gue8b1FZKaz7VYRafdOSEUirTthwzCMLkIYofEN4JvACmC58/wcEakAvt/Kz50EjFbVMcBnwJUAIrILcfPXrsB44K8iEnGOuQO4EBjpPMa38rNbTZFANGZiwzCM7ktGn4Zzw75EVY8P2OXN1nyoqr7oefkOcJrz/ETgP6paDywQkbnAviKyEOilqlOced0DnAQ815rPby1FIsTUhIZhGN2XjJqGqkaBvfM8h/NpvvkPBZZ4ttU4Y0Od56njvjhlTqaJyLSVK1fmbKJFRYLJDMMwujNhoqc+dNq+PgJsdgdV9fFMB4nIZOJRVqlMUNWnnH0mAE3A/e5hPvtrhnFfHGf9RIDq6uqc3eaLBNM0DMPo1oQRGlsBq4HDPWMKZBQaqnpEpu0ich5wHDBONXEnrgG28ew2jHi4b43zPHW8XTHzlGEY3Z0wGeHfzvWHish44HLgUFWt9Wz6L/CAiNwCDCHu8J6qqlEnyXAs8C5wLnBbrucVYt6YH9wwjO5MmIzwnUTkJRGZ5bweIyK/auPn3g70BCaJyEci8jcAVf0YeBiYDTwPXOr4VQAuAe4C5hIv096uTnBwzFMmNQzD6MaEMU/9Hfg5cCeAqs4QkQeA61r7oaq6Y4Zt1wPX+4xPA0a39jNzQaTIzFOGYXRvwuRpVKrq1JSxpnxMpqNTZOYpwzC6OWGExioR2QEnWklETgOW5nVWHRQRWF/XyPcf+KDQUzEMwygIYYTGpcRNU18RkS+AHxP3L3Q7ipzKJU/P6JYy0zAMI1T01HzgCBGpAopUdWP+p9UxKWr3aleGYRgdi6xCQ0TKgFOBEUCxWydQVX+b15l1QIrav0aiYRhGhyJM9NRTwHrgfaA+v9Pp2BSZqmEYRjcnjNAYpqrtXlG2I2IywzCM7k4YR/jbIrJb3mfSCTDzlGEY3Z0wmsZBwLdEZAFx85QA6vTC6FYUoO+TYRhGhyKM0Dg677PoJJjIMAyju5PVPKWqi4hXnj3ceV4b5riuSFMsVugpGIZhFJQwBQt/Q7wi7ZXOUAlwXz4n1VFpshoihmF0c8JoDCcDJ+A0YFLVL4lXqO12WIVbwzC6O2GERoPTJMmtPVWV3yl1XLyahlq1W8MwuiFhhMbDInIn0EdEvgtMJl4uvdvh1TRemrOigDMxDMMoDGEc4TcDjwKPATsDv1bVdu+a1xHwWqeemt7u3WYNwzAKTpiQW1R1EjApz3Pp8EQ9JqmvbN0t3TqGYXRzChI6KyI3icgnIjJDRJ4QkT7O+JEi8r6IzHT+Hu45Zm9nfK6I3CoFyLTzmqcqSyPt/fGGYRgFp1D5FpOA0U5W+Wc0h/OuAo5X1d2A84B7PcfcAVwIjHQe7V4Py+sIb4xazoZhGN2PgggNVX1RVd2Wse8Aw5zxD52QXoCPgXIRKRORwUAvVZ3iRHLdA5zU3vP2ahoNTSY0DMPofgT6NERkJk6Ybeomclt76nzgIZ/xU4EPVbVeRIYCNZ5tNcDQoDcUkQuJayUMHz48R9NM9mlMfH0+3z98ZM7e2zAMozOQyRF+XFveWEQmA1v7bJqgqk85+0wAmoD7U47dFbgBOMod8nmfwEQJVZ0ITASorq7OWUKFN3pqw5YmVNWKGBqG0a0IFBpOnalWo6pHZNouIucRF0zj1JMpJyLDgCeAc1V1njNcg2PCchgGtHvMa2pGeH1TjPISc4gbhtF9CFN7aqyIvCcim0SkQUSiIrKhLR8qIuOJ17M6QVVrPeN9gGeAK1X1LXdcVZcCG525CHAu8Y6C7UpqwcK6hmh7T8EwDKOghHGE3w6cBXwOVAAXAG1N7rudeP2qSSLykYj8zRn/PrAjcJUz/pGIDHS2XQLcBcwF5gHPtXEOLSY1YKqu0YSGYRjdi7DJfXNFJKKqUeBuEXm7LR+qqjsGjF8HXBewbRowui2f21ZiKfWmak3TMAyjmxFGaNSKSCnwkYjcCCwFumXRwiuP/gobtzTy3sK1gJmnDMPofoQxT30TiBA3HW0m3pDp1HxOqqMyclBPHrn4ACoc57eZpwzD6G6E6tynqnWqukFVr1HVy1R1bntMrqNy3wX7AVDb0JRlT8MwjK5FpuS+h1X1jKAkvxwm93U6EpqGmacMw+hmZPJp/Mj526Ykv66IW6zQzFOGYXQ3As1TTm4EwPccE1XiAXyvfabXMXGFhkVPGYbR3QjjCD/SZ+zoXE+kM1FeauYpwzC6J5l8GpcQ1yi2F5EZnk09gbf8j+oeVFr0lGEY3ZRMPo0HiGdd/x64wjO+UVXX5HVWHZziSBGlkSIzTxmG0e3IVLBwPbAeOEtEIsAgZ/8eItJDVRe30xw7JOUlRdRZyK1hGN2MrBnhIvJ94GpgOeBWX1Kg24bcAlSWFpt5yjCMbkeYMiI/BnZW1dX5nkxnorI0YuYpwzC6HWGip5YQN1MZHspLIhY9ZRhGtyOMpjEfeFVEngHq3UFVvSVvs+oEVJZGzDxlGEa3I4zQWOw8Sp2HAVSURti4xRzhhmF0L7IKDVW9BkBEqlR1c/6n1DmoKImwYkN99h0NwzC6EGHave4vIrOBOc7r3UXkr3mfWQensjTCpvomfvPULBasMlnaUv7yyly+e880ABas2sxj79cUeEbtw/q6RqYu6NZpTkYnJ4wj/M/A14DVAKo6HTikLR8qIjeJyCciMkNEnnB6g3u3D3d6kv/MM7a3iMwUkbkicqvTK7xgVJQW88W6Ov49ZRHfu/+DQk6lU3LTC58yafZyAE7561v89JHpNKX20+1izKxZz+7XvMgZd05hi/nDjE5KGKGBqi5JGWrrFT8JGO2UV/8MuDJl+59I7wF+B3AhMNJ5jG/jHNqEW7QQoL7JbgAtQT1tc6MxZW1tIxBfhXdlbnzhk8Tzhi4uII2uS6iQWxE5AFARKXVW/3Pa8qGq+qKqul7kd4Bh7jYROYl4xNbHnrHBQC9VnaLxO849wEltmUNbcXtqADQ02Q0gLC98vIztrnw28bq2oYniorjS6AqPrsrQPhWJ53bNGJ2VMELjYuBSYChQA+xBbkujn4+jVYhIFXA5cE3KPu5nu9Q4Y76IyIUiMk1Epq1cuTKHU22mwqNpNNqqMTR3vDov6XVdQ5TS4vhluK62IfT7rK9rZH0nEzL9e5QlnpvQMFxmf7mhU1krwgiNnVX1G6o6SFUHquo5wKhsB4nIZBGZ5fM40bPPBKAJuN8Zugb4k6puSn07n49I6yaY2KA6UVWrVbV6wIABWU+wNVQmCY3AqXRqtjRGk0xJuaAo5T9Z2xBNaBotMU/tfs2L7P7bF3M5tbxTVtz8c6s3oWEAy9Zv4Zhb3+Ca/80u9FRCE0Zo3BZyLAlVPUJVR/s8ngIQkfOIdwX8hjbfmfYDbhSRhcTLl/zSqX1Vg8eE5Tz/MsTc84bXPNXYBW8Aqso+103mu/e8n7P3XLFhS5qk39zQlBiburDlUUVL19e1eV755p35q/nt/2bTGGs+e9M0DIAp81cB8MC7ixMLtLqGKHe+No/N9R0zDyxTP439gQOAASJymWdTLyDif1Q4RGQ8cTPUoapa646r6sGefa4GNqnq7c7rjSIyFngXOJcQgiufeM1TXdGp2RCNsbG+iclzlvPVm1/l0J0GcPUJu7b6/eat3MS4P76WNl7XEE3ojHe+Np8Dd+jPITuF1w4vf2wm95y/b6vn1R6cOfEdAC4+dIfEmAkNQ1X5yUPTE6831TfRs7yE656Zzf3vLmabrSo5ZrfBBZyhP5k0jVKgB3HB0tPz2ACc1sbPvd15r0ki8pGI/C3EMZcAdwFzgXmkR1e1KzsO7JF43hV9Gl6T24JVm/nX2wvb9H5BiZCbG6IUeWxW5/5zKvNWplong+lMYbq1nlL6ncmGbeSHVBOl+4tbun4L4G+T7whk6qfxGvCaiPxLVReJSM/4cJq/ocWo6o4h9rk65fU0YHRbPztX7Dqkd+J5rAu6NA658ZWcvM/6ukYeeHdxkg/IS11DEyWR5LVLbX34G2osxz6XfLLJU3bGNI3cMH3JOkoiRewypFehp9JiUgueumZudw3VUS0YYWpP9RSRD4GtAERkFXCeqs7K68w6AUN6l/OlsyroSqgqazaHj2TKxO+emcND01LTfJrZXB+lJJK8poqkessz0IlkBhs9Nur6DnpD6Gyc+Jd45+mFfzgWgN8/O4chfSo474ARBZxVOFILnq7e3EC/HmW4ecsddWERxhE+EbhMVbdV1W2Bnzpj3Z6S4lC5kZ2OLY3+F+vi1bW+40G8/tlKPlyyNuM+tY1RoimqWkty/Tu6zPBmfs/+ckPi+UqrW5YX7nx9Pr/578fZd+wApAqNo/70OtCsaXTUqMwwd70qVU3YKlT1VaAqbzPqRBS3YEXcmQgKfT32tjda9D7n/nMqny3PbM2sa2hKExpNLfix5DokONds8HyXX6yLR3r1rijhzbmrCjUlo4MQ1I+nKKFpdEy/VxihMV9ErhKREc7jV8CCfE+sM5Bqi28ry9ZvYdYX69m4pZHT//Z2ixzCuSRIaOSiFPy7vxyX9HpzfTQtkKAx1vx6U31Txu8hVWbM+mI9r366os3zzBV+3+XOg3qybEPXM2saLSOo/pgrNK7+32we7YCFPMPc9c4HBgCPO4/+wLfzOanOwg4DmiOofv/snDZHxBz+x1c57rY3eePzVby3cC03PPdJ9oPyQEsys1tKWYpJr7ahiaYUTcOb9/Ltu6f6huq6pDrCj7vtTb5193tMeGImsVZGKLz8yXLm50hg+5VGGdCrjFUbzTzVFmZ9sZ4RVzxT6Gm0Cb8mbu8vWpNknn125tJ2nFE4AoWGiJSLyI+Ba4nXgdpPVfdS1R+ramZDdTfhhtPGcPDI/kDclnrfO4vb9H5uz3E30qhQTZ5yUTjQazbyWvHKipujqCqdRlZpQsNjnnpvYfxSC3IKBkWY3P/uYpa2cjV//r+mcbgjqFZs3JJkPrv5hU/Z5/rJod/LTwAP6FHGqk3dU2h8tnwjI654hrkrNrbpfV74eFmOZlQ4/MxTp94xhZlfNHfXTg0S6Qhk0jT+DVQDM4GjgZvaZUadiB5lxZy5z/DE69ocZXCuc1anmxs6l9Dwlh2Zs7T5puCtuVRaXMTg3uUA9OtRyoYtjWm5Fl7zlEttwHeRyf8RaUX1fK+w+/jL9ex7/Uv8/JHmBKzbX5nLyo31oTPR1/l8l5WlEWobmr+rT5Zt6BI3wTA89dEXADw/q23nu03fylxMp12YtnANB9/4Mp8uSxaUQe2iF3kCTnJtAs8FmWa0i6qeo6p3Ek/ma1MPja5Kn8qSxPNojpyyNWvjF02hWoYECY2BPcuSXt/20ueMuOIZojGlZm0tX7nqeR6cGg+vPebWZqd5RWmEvYbHW6ZEioRHLt6fP5yyG/WNMZ6duYyYwvkHbpfY368sy+YAp6FXIKUSlHT5weK1gdu8kWPH3vomAM/NWsaazQ1J5pCVIc1LfkUVK0sjNMWUFz6O9xMZ/+c3uOje3JVr6ci4X3tRG4NISorTj891ouei1ZvZlIOF4NSFa1iypo7HP0j2T4TpqVLaCqGxqb6J/07PX5WlTDNKXO2eMuZGCl6hkaskPzfKplDRWas2pZtUDtt5QNqN9v9e+hyI35znr4x3L/SzwVaURLj3O/vx0k8PBWBY30rO3Hc4Kzw33n49Snnhx/F1SVNMefyDGuauaPYrLFmTHO7r+kbUE3Sbqo34CYZPl23klL++zY3P+/uLNtb7C8w5SzckvU41qfkxdcEa7n1nUdp4uVO37OL73u/SzZjenruKEVc8k+jQCM2aXFEbFkTXPzM7qfyGS66LQB5606t88x/vtugYv2i+kqL4teo1c6oqT30Uv7H/6tjg+q+tuaVc/d+P+eGDH3Lp/R/k5frKJDR2F5ENzmMjMMZ9LiIbMhzXrdiuf3P0cdTHrNIa3JIbLUlyyyVTF6xOG+tdUZL2o3RvnDFtvnX73QvKSyJUlRUnBQ6k7hspkoT9tjEa47KHp3PELc0O8DMnvsMfX/w0oYW5P8DGpvjfuSs28dzMZJOHX5y760uY9UXzJdwUjXHVk7N4e+6qpKztoLl6Pz8TZ9w5hcVraunfozQxNmpwryTn/fcfaO762NHDh1vK847JbdLs5Tz5YdwsFUsIjda/79/fSA/erGuIJnyCucD1oX24eF3oY1SV6usmc9WTyXnP7rXTFFM21TexalM9z89axtvz4r+zgb3Kg+fRCu3JXYw9M3NpXu4hgUJDVSOq2st59FTVYs/zzpeznycqS5uT6oOS4lrKKicbO8yNKSwLV21m72snpdlV/Vi8po6z9t2GqROaw2N7lBVT2xDlysdn8vIny5P2z7bq3rDFf/X+yEX7J56v3lSfsN/+6D8f+e5/28tzufT+D/jNU7MSnzl14Roe/6CGI255jZ8+krz69NM03J+QV0O5+L4PuPedRZx917vB5oiUU2xJLknvimZt9IEL9ksSZpPnNIcHN0RjNDTF+MebC7pEPbPFHu1wek385uteKm3RNLyH7uqUDxn16+dbFKCQjS/Xtbx6cs3aOlZvbuDedxZx75SFqCqxmCZ6wsdUGffHV6m+bnKivhRAaQZnd0uzwtfXNvL6Z809hPLhE+l4XpZOyA/HjQTI2UpntbMazqVqeefr81i9uSFrDkM0pqzZXM+AHmX0r2r2F7iRHg9OXcz5/5qWfIzPDdS7wAn6Afby3Ew31UdDXeBbGmP8e0qyyeeyh9NNFZDZXOFd1E+e0ywEF/lkvdc1Rnl2VrLZbX1dQ2Dp6k31TUlmDa/Q6FtVGngjqGuIcu87i7j26dncOyXdrNWZqGuI8uqnyTev1Zvq+cebcS2hLf4673VSszb3pfHnrdzEYTe/CrTMROw1t1711Me89tlK7pmykBdnx6+vppiy3LEi/Pbp5v4ZQQrm7sN6t1hoTJqzPPtObcSERg647Mid2GarCupyFO20qg1C46+vzmXy7PQLxzV59cvgOIa4wIopDOhZluSsLM1QMiWqmmaa26qqjL237QsEr5Z6lDVrafVN0ST/UBBNLTABPj8r3b/iBisE6QneVZqX1HDqi+/7gP1+95LvvtOXrOONz5szvqvKkku87RpQXK+uMZq4hla2MCR3yZpaDrrh5RSPWI0AACAASURBVFatkHPN6k31fOOud5LGiouEKfObzZ5tsZp4ncNBQRv3TlkY6r1UlRUbk0Ozl65rfu3tm5ONjSkadW1DlHmOrw/8F1cA26eYbV1Ki4taLDSqPIVBjx2Tn7LqJjRyRGVJcc40DdfM1Rpz143Pf8oF90xLG3dDP7NdhDNq4jHiOwxMvpAzHReNaWKuIsLi1bVsqm9M2PKDrFfem+nm+qaEgziVIb2bbb7eH2E2/v7GgjQ/Qb37nQbMqSW5MUGmrNTVaaoGddSuW/OVrXumHVfXEE18By1dMDz03hJq1tbxWAfIIP7Hmwv4IMUXUBwpQjzFvttiefXmLhy5yyDffa56Klz9qfveXcy+17+UZLb19spJFfiZ2JBy7cRUk3wKQYU7d966JzeeOiZtvLS4qMU+DW+jr1vP3LNFx4bFhEaOqCiNBMZdh8Evbj2X5ik3ySyonk1TNMYj05YkEov2Gh7XEm47a0/+fm51RlNPNKaJbHhV5ZCbXmFLYyxjOCwkaxreUvOpDO5TkfF9MlHfFGPhqs0JB3q2CJsla5PNU61JrqpN+b/5vUP1iL7pxzVEE7b+sAuGJWtq+Wz5xsTKvSOU6ffzxRUXSZJ20ZbwWK9p6/en7Oa7T1VAKf5U3nI0Qm+pGq8/afhW4fNBUjWNaExDO6KH9U2/xksjLdc0vBp/vgJpTGjkCDdhqyW8v2htwsY74YmZadtzESPu4qrxQSuXe6Ys4uePzuD/XvqcQb3KEive43cfwpG7DEq6wacSjWliBe/9wbmfedCO/X2P817Ul341uMXKjhnU92xsrm/isJtf5aAb4jU3XUHsOsLdqB6Xj79MDgx8yOOsD0tqpq+f+f6Xx4xKmO9cjrvtzYRpMqz/8uAbX+GoP73OrS/PBTpGfxG/wIjiiCR9D21x9A/1LCKqSv2vy+0GZK+puq62gc+dzHTv9+a9UVeVtcQ8lfx7Vc3sExk5sAev/fwwwH9hVBIpYu6KTYl7RBjaozKuCY0c0Rqhceodb3Ot4xD72uitk7b1rYyHuLZE2wgK2VTVRJZ50MrFW+5iiM8FPOG44FjyppgmCvB537+hKcaUKw/n7+dWZ527K0D8Vkc7+5hygMDGTl5S/yeupqEa70j444f8I7VchvWp4P4L9sv6OZk+szGqvHXF4Tz+vQMSY5WlxXz34O3TjnXNJJ8t39QqTbPwIsNfi/hibR2/8oSirqttZForesIDDOqVXGHAd5+ewWGsLqf9bUrC3OmVc+41XBKRFt2EN6T4V3779Owk53gq912wH9v2iwu3oSm/uX99ex9Ki4uoa4xy7dOzAysieFHVhC/tz1/fI/S8W4oJjRwxsFc5c5Zu4K2QJa+9xQ1VFVXoV9Ucz+86x+av3MyIK55hxBXPZK16G3SBb6pvru8UJDS8Tm8/raJXeUmgqr65vok/T44n+nlNB/VNMQb3rkiyEWfDT2gcvdvWPntC38pS3/Fqzwre63yFZk1DJLv579gxgxnYqzwp+smPtZsbkppWpQZELF1fx9A+FQmTX/P809/3g8XxWltTF6zhJ1kEmh+5yPWYvmRdm+qPeTWNiw6JC8b7312clDR615sLOO1vU0KXY3H5cl1dUphy6vUy/3fHsPuw3gnb/pI1tbzyiX/EoDd51Pu9uVpQVVlxKI0o/vtVNm5pSvqfrtncwBMpmqzLV3cewCBPfkZpcRG/Ozluavvm2G05bOeBSQKxPoS58r2Fa/mfkwl+0Eh/7T4XFERoiMhNIvKJiMwQkSdEpI9n2xgRmSIiH4vITBEpd8b3dl7PFZFbpVA1NgK40Fk1hu2TcNtLcxPP65tiNEZjSY7grR3nrzcc9Ppn5vi+V31TlA1bGtkS4K/wrnaCOsZ51ejUSrQuQeYFb/E9r2mmNVV/bzsr3Xk3oEcZo4cmRxwN7VPBIQE/jEcvOYBT9hoKwOWPzUja5obJZjJtJTQY5z6SSWjcO2Uhe147ib2unZQIs3V9W29e/lUguUSKl62q0oWetyrum5+v4q+vzuXcf04N/PxU2mqeeuqjLzjxL29xsaesyTvzV6fZ6zPhzWHZd7utMu4b1Dvejy2NUQ74w8uJ137+pqIioaw4kvDdjbvlNb79r/eyvneSecoVGqXhhMYhN73CpQ98wIYtjfQszx4BCOC3vjt7v+FM/81RXH3CrkDy7zDot+3FG1noZqHng0JpGpOA0ao6BvgMuBJARIqB+4CLVXVX4DCay5ncAVwIjHQe49t5zhkZ0b+KYX0rQoc8frq8OVpjw5ZG5izdQElEEqr3YGcV4n2/+qYoj75fk1hNuJx6x9uMufrFwNWIt05SkKYR8VxkZQFRTN4CgF5BtNZj2vL6YVoT/fW1XZu1iicvPZB5vzuG4kgRp++9TdJ+5x+0XcZY/7P3jReS9N5Do05GLsBbc1fz0RL/bF/XVOBqD70zhAJ7o3Rc08DS9VuoKIkwpHcFC/9wLGfuO9z32F5ZNJiGaIwbn/80LQzYvZH5VUltqyPcTax0NbTahibOnPgOF/w7PSIvCG/BSW9VYz+CEj9TmfXFej5JSUz9/PpjfPctKW42K4V1JHtlg2vCrCqL8MmyjazYsIV356dXSXBZsqaOZ2cuY+OWJnqWZ462chczqaYsl94VJQntyRta7P62H5m2hA1bGhMaqRevsI7ksTpuQYSGqr7oqWf1DjDMeX4UMENVpzv7rVbVqIgMBnqp6hSN65H3ACe1+8SzMLh3OctC9gz3qrH7Xv8SH3+5gcaoJsIzXU1j4ermENOq0mJ+9sh0fvDgh0nv5ZbECFrZe2/w9/nUQoJkx2t5wA/duxo70LPiW+0xO3hXpNuHcEZ+eNWRfHDVkUljJ+4xhF8dO4o9tumT+AGdu/+2PPjdsQnTk6qmORkP3WkAZ+0bFy6VPg7SDXWNSX26g8JT99+hH9DsBO3pMdf9/Gs7ZzyfaEz5bPlGdhrUI2tRvmzhnN5Irw+dm8STH37ByAnPsWj1Zk654+20Y9pqnXJNk64W5N50312Q3f/wl1fmcsPznyS1tc0WfXbPlEWhVvPH3fYmJzn9wLMRNurIq23GfMxTDU0xahui7Pu7l/j6xHd8Fxles9a62oasQmPU1r0S+7ZkfvVNMWZ9sZ6fPzqDMVe/yCl/fTttMeHVSPNZt64j+DTOB55znu8EqIi8ICIfiMgvnPGhgPcXXuOM+SIiF4rINBGZtnKlf7JWPuhTWRraFux3U4vGmoXG4N7x1a5XCHlLL/g5xoJW9t4L1F2BxWKaFBqZrGn4XxZBq9jrHLNZ/x6liQiSbx84gutOGu1/gIe+VaVpZpr/O3NPLkhxEosI++/Qj923iVsyVeOx/17OqN6G358Sj3f3i3pZV9eYFOEybVHyau2M6vjapbwkwr3f2Zc/OLHzXo3m7ACtwaUxGmPt5kYG9MwcbgxQ2YLEsZP/+jZPz/iSp2fEExZf/mRFWhFFaLtPYw/n+y0uEmIxbVERwJte+JQ7Xp2XKLgJ8aipTEyavZw7X5uXcZ9s5/SN/ZL/J6XFRdQ3RZOOW75hCyOueIZnZjQnfHoT96IxZX1dIz/6z4cJoZcavei3IPRGIy5ZW+f7u/ayo5P/dGBARKEXb37PlsZoWlTaHyd9FurYXBM+c6WFiMhkwM+DOUFVn3L2mQA0Afd75nMQsA9QC7wkIu8DfgUSA68kVZ0ITASorq5ut4CS3hUlfJxFaDRGY4G9ExRNrMwGOmaq5R6br9euuba2Me0Cne25iVz79GyuOm4XoPniL43Ek4Vm1KzjhNvjq7b9ttuKv59XnbQyCaNp+DGgZ3nC2XnITgOy/oBag7d2VOpqqtwj7Pyc758u28imAHPIfy4cy5fr6nh4Wg1frKvj4JEDkra/8YuvxpO1stwEm2LKlqZoVrMMJAcf7D6sN9Nr1mfYG77/QLOGec3/Zvvu0xafxpR5qxM+uRUb67n0gQ94LiV/aEtjlGv+N5tvHziCnQb5R7V5BXNDU/b5LM/i18gWlXj9ybtx5C6DWOKUFCmOFPHZ8k2M+vXziX3cwpAPTl2cyJSuKIkkFnmfLd/I9CXrEpVnAS46ZAeuf7bZj+inyW9paBYaKzfWZ80gH9irjClXHk6/quyLilRNI1WLmb5kHbGY+mq0+Sx2mjdxpKpHqOpon4crMM4DjgO+oc1LghrgNVVdpaq1wLPAXs74MM/bDwPyVzC+lfSuKGF9XSO/enImYwNKTNz20ud8/4EPA5vQuCuE8uIIkSKhrjFKkcCew/skOUnXbk5Xb2d41Gc3tnt9XSOfLdtIpEi47KidABICA+JmhynzVietyoI0jWyhB97VdVmeVjruDySm6T8MbyCBN37/ptPiGsPF973PZ8v9I9CG9K5IrP5O2TNdid1mq0q27VeVVe1visaob4wFfodBPPX9g0LlnWSjJT6N+Ss3JZWcP+vvyaU/UgUGwH+nf8mDUxfzr7cXJo0HRaJtaYry7A8P5jsHNQcDnJzy/WYSdPNWbuJbd6cHAvztnL2TXh+280C+OXZbgIQ24dW83Q6Q7qIsFtOkz737rYVJGhLEtWUvfj7D2sZkbSTb/71nWQmDe1eE+l97fYgPT1vi+/trbXfKtlCo6KnxwOXACY5wcHmBeAn2SscpfigwW1WXAhtFZKwTNXUu8FS7TzwLvcpL2NwQ5b53FrNsw5a0pisAXzh1bZYF/LNvPWtPTt97GKMG90xE7fQsL6GyNMJ0j1BY62MTXeqjPh9/25s8+dGXRFN+JF5WbKznH281JxAF3RgnfrOaQ3Ya4LsNSCoBXpKDG6AfCU3DJ3HKG+XkXfF5ww9Tbwwug3qXMahXOQv/cCzjRvmXpoDsK7gmJzs+qCRKJlqa/etHSxSNw//4Ggff+EraeKbVcr0jHFJLyAc5dncd0otdhvTi9OrmNd+fUnIIMk35vH9OTdzwvYwf7R+GDfCtA0YEbnNv1qk5FL3Ki5PqhUG6+fMXj83giQ+Tf9OpWpDf//3UvYYl6qpVtiBZ0NvU7dH3a9JK/0OzsL7o3vCBCm2lUD6N24GewCQR+UhE/gbg9B6/BXgP+Aj4QFXddmmXAHcBc4F5NPtBOgw9UtTHoOqrmdhhQA9uOn13iiNF9HFugmXFRWk/Ulfr8PolZvvYuL3lqYPKeV/15CyWrGm+mQaVOh81uBf3nL9v4Ny9N5vWdBwLg2vyKi0uSvyoz9p3ODedNobRQ5tLkXhV9n5VZbxz5TgyEcacBFDsE8q4s8dM0xSN1+EKMvGl8sfTd+df394n4z5+oblB/POtBfz4Px9mDchIzRPy0q9H8Oc1ONdQqu9uSUq12UsO24FZ13yNgU6SnXs9+MncB95dHJiD1JoqtlefsCsj+vnnFLnXzP3vJgeEpNaNCiK1+VNqBJvf//26k0Zz0h5x7SpTZYVUUvvIu03PvNQ3xthc35ToAtkeFCp6akdV3UZV93AeF3u23aequzqmrF94xqc5Yzuo6vc19UrvAIRROTVlXXXmPs2hpKln5IZklpdE0npNu87tdzyhgAtWJRfze/Hj1KZE2Veyp+09jAs8poSW4L3xtiShryVcdOj2/HDcSM4ZOzwRx967ooTTq7cJPMYbytxWiiRuXjnDs3L+xfjmiKrGaDyLP6x56tS9h3HYzgMz7tM/w03cjyc/+jJRaSAIbzJiqsM3UyVk9xpanqIpn5oSyTV6SO+kG6Rrdu0TkJB5y4txp+4dr87jhynRgan8O8PCxeWIAG3RFV5hM7138IkAjMaUx96vIRrTtJDh8pIizkoJligrLuKq43Zh+m+OapEG6gbDZKK+KdqmmnetoSNET3UZ/Oz42W7U3lV9asSDG+K3dnMDuw+LR7U896ODnbH4xZrJz3ChJ0FrQM+ywBpQXm46bUzW8uku3uqzfzl7rySh2ZKS0i2hvCTCZUfuRFlxhLP3G863DhjB9766Q8ZjRAQR4aiAiqgtQUT409f3SPg/jhszmHGjBiVMIvVN8SiXsJqGlyCzSiZNY/yu/maabGa0zfXNN5p1tY1JZrtBAZFfB/7h5UStrnkrN2VsEta3KjkHxb029h3hn+zn3nxveP6TjP2tRw7swaEZTKQuQebRloaiPv2Dg/n4mq8ljd07ZSE/fWQ6D05dnGaWKy+JpBVRLCoSIkWStbJAKhcesj33fiezgKxvinHKX9NDr/OJCY0c4qdp3PTCp8kDKb+zmOeHV56yOt19WNzcUt8U46bTxzDpJ4cwanAvepYXJ3waqXbwoJvFKz87jP2275e1amdLEu3dzNUDd+zHsWMGJ3UgC1MXqq2Ul0S4+oRd6RUyCzefqqlbfPCIW14H0v+XYRg5yL8wYyYh3r+nv0DJ1JukKRpLShBdV9uYlAEeVJ7li3V1iQS7xqgmhX3vNjS5SnGqoBvUq5wHLtiPW76+u+97h63bFvbGG1RF1++6zGRKrSiNpOXTuBGCqzc1MHVBsr+lNf/3IEoiRWnfayr1TbEkE/SwvvGE0nxiQiOH+AmNbElRTUlCI/mC/v7h8Y6AInHTz0jHdt67oiSxwkkVGpeP908+c00FfjLBryxzGKpHbMVWVaX87Kj4Z3o1pZb0IcgX/VJuXKnmv3+cFy+k2Bb/iytkU5PYWuMID5rHgAxCo0iEXx2bXkwyaPHw++fmsOOE55Ls4+vqGpKSMouKhJd+emjW+dY2RHl77iqOu+2NtEZGW/kIngN27B8Yhv3+orWsr82e45Qti97Fm3DqpawkkiRQykuKOHXvYb77ZmPZhi380wkgcYV0a/7vmcj2fs/MSNbKzshgps0Vhf9ldyH8fvTe1beq8vrnycmGTUklF5KPd4VQ6o2+oiTeu2PO0g1pDryv7jyQTfVRbvVxmgFJIZYAuwzuxVn7bhO6aY2XrapKk7K5vSaBoPpV7clLPz00xWEblxo/PmIkFSURxo0axIyrj2pTv2rXtZbqIG/NijOoTlOq8PMiwAUHb59IsHQJaiZ152vz08YemVaT9B0UFwk7DOjB1cfvwtUB+SBAYOdCCPZdZGLV5mbH79X//ZhfO3lGXsJqsI0BprOGpljSd9UYVeZ7nPB7bNMnsLwMJDvyazy9V7bdqpJ1tetbZZbMRLYFzcPTmqO5rj95NGfukzn5NBcU/pfdhfCzo7qr25q1tWx35bNJlT4hOaIptWlRc8RJeie452Yt4+j/e4MrH08uyNe7oiRjo/rU39KAnmUtzh49LqCNpPd9OkI9yT6VpYnS09AshHcb2puLDo37QXqVl7QoosXloB37s1VVKRc775MaVBU2GsvLtv2qWPiHY5n3u2P49LrxiXyGsKtrL6nRdkCaNuDy3+lfJmV+u1nh3lXugTv2C/3ZPxo3slU5J7d4Mpz/9fbCJLOLS9j3vSqglH9jNJaUIxWNaZJp7MlLD8z4viKSCGZxQ3Qfvmj/xMquPMdm2aBSNAf7FOsc3Ls8r0l9LiY0ckiJzz/MzY0ISuZzt+85vA83pLR8LC2Ov1+q0PCG1qZGgQzoWZYWX56JfbfbKlGZM1N8u5fbz97L127ami537ck1J4zmOwdtF8qRmo1+Pcr44KojE2G+dQ3JZsK22LYjTqXWm04bw8yrj0rT2r62a7pD/8qjv5L02hsR9ffX5zP+z6/z8Rd+hRXiuI7wEf0qE0X1vBFpPcvCC64dB/r7Zrw888ODuOHUZIext8QH+BczDKvBDuxZzqWeAIn7L9iPQb3KaIpqWo6UG+30w3Fxc/C7vxzHexOOCHzvT1MKJ+48qGcif6g8ZX6pddVyxXb906O6ghpS5RozT+UQP4XY7YsR1Izl2DGDmTxnBbeeuSd9U8wQpZH4qiXsrfiY3bZGREJrDn8/t5rDvxIP97zi6K9wjpNR21rceaY2lOkoDOhZliitkms2N6RmBrd9xVkcKaJnpCgtfPm2s/Zip1/F05TcGl0XHboD5x0wgslzlvPwtBrW1zbw04ens2RNLVOdZkdhStVv268qoSV6V62Zcjdczj9wOxRNqlQcxK5DerPrkN5c/lh6x0qXp1OECLTM/+T6C3921E4cuGN/SiJFbKxPF0Rn7zecsz31q7x9LlKJxjSpn0dZcRG9KooTJuRUH0RLcmzCcuXRX+Gb+2/LPVOac02OGDUwaxn6XGFCI4f4ZY48+n5NUvKXl9d+fhjb9qvi+DFDfLWDIJ9GENXbxi+asCv+Iz0hqK6ZpS24p3/Yzm1fyXc2UrXBXNq2vc7j43cfkmSi2cYTDVdeEuG4MUN4btYyatbW8lhKRYI1m7M7moOqBmTr9w6wVVVJIngjF0x8Pd3/0hKzl+uMdxdspZGiRBLrmftsw3/eW9LmOTbFNMkUm2tHuB9H7jKIytJiDtqxf6JW2Ml7Dms3k7CZp3JI0A/u+mfnJMoveHHt7UHmJPfmn63EdmL/4uTEpTClyfNBB3BntDsn7TEkKeEvl6a6IOevn4kC4qXc/XwafqVnUgm6hr1Vg/tVlXLVcbswKyV/Icj5nktaIjSO2S3ueztxjyGJY90E2ON3H5KT+aTmGOfDpfC3c/ZKynB3gy7+cvZeibH2DDwxTSOHuKuxHQf2SGolCfBvjyoZloSmEXJ/d2W12OnBceY+21DbEE1rM5ov3N+PhJ5x16E4UsSVR49KRLPkUnD6CY1Z13wtMFGtZ3lxWgUBSM/89iMWkIu6jycprzgiiQKE54wdzrL19aytbeC8kD6xlnLe/tuyzVaVXPfMnFD+EpdttqpM8r3tMLBHIs9k697Ze4iH4VfH5sfc6WX86MG8/MkKFq6OBwa4lZZ7V5YkKle3h4bjYkIjh+wypBePf+8A6htjaRVDW0NQ9NRXdx7AK58mh+7+/pTdONop4rbH8D78e8oiDt1pIDtv7W8am/STQ9o8v1TcVWo7BHB0SJKjnHL3JYzwRIC5K9tMEV89ykp8ix/6aR9PfO8ATvZkFJ+xT3LOwj/Oq2bxmlr2HN6XbftVsmh1bZIp87qTkp3Z+eDbB27Htv0q2XVIb8Zu33q7/aiteyac7Vv3Kue0vYcl+rO0FlcDcxNMXT/Qg98dG6psT1gmHLNL84LEM+4WNfTrH5MvTGjkmL2G9+Wz5Ruz7xgG5+pItVVOPLeakRPijtDJlx3CtIVrk1qKnrTHUMaNGuSbKf349w5g5cb6RKJgLkloGt3RPkX8hnHEqEFMnrM8qTNjW+lbVcrUCeM49tY3Q/meUgtnuvj1+d4zRQs9ec9koeGt+OtqNu2duFleEkk04WoLAz0O7qqyYm4+3T87PRPXnjSaq56clXjtmopuOn0MD01dkghXbutcU+ldWcKwvhXUrK1LCtN3S7mEafqVK8ynkQf6hIirDxPe2reylKF9Krj2xF2TxksiRZy5zzbsNrQ3Ow7smdaDWkQCS2vsNbxvqOiW1tDhKkgWgL98Y08e/94Biai5XDGwZznvTTgiqZJvED0DbuquPf80JwP6uwfHTUwn7RHOvu92qctlaKcbJJLpBp6rOmaZoqKycda+w6ksjaSF1bv5OQN7lvODcSPzumD6weE7Av6RbGECFXKFaRp5wDVTDOhZxsqAUNvfHJ/dFloSKeKtKw733faHlJyOjkAHLDzc7pQVR9rNhxREUJ/q6TXr6V1Rws2n787Np++e+H+5ta1Oy1JOwzV55dIU8uSlB1Lb0ES/HmVc+/Rs33bJLW1oFUSmzPps/P6U3fj9KbvxaEpf+bZUE2gpX99nOF9Pyfg+aY8hPPnRl+3q0zBNIw+Ul0T42zl78/QPDvLdftjOA7qkCeeY3QbTp7KkzfkeRtsIMk+Bf9b+z47amT99ffdEh8Mg6hNCI3drzYrSSEJo3XKGv7aRq8igXJTrT42Ka0+h4ccfz9iDT64d366faZpGnsjUWey3J4xux5m0H0P6VPDRr48q9DS6PZlu6n4lZipKI2m+DD/WOQUFg6rgtpUhTlLosL4VrNpUn2jXmqsFVi4qL7s+i+IioSmmBRcakSIhUtR+WgaYplEQfJq/GUbOyJRYmIs2vNsGdMVrK6MG92L2b7/Gm5cf7tshsa1UlrR9jezWBxs3Kl5JoTtGChaqR/hNIvKJiMwQkSdEpI8zXiIi/xaRmSIyR0Su9ByztzM+V0RulU5s38nHD8IwXDL5ABatTi8CGJZ/nFfND8eNbHGBy5bgZr8/dNFYth9QlaiDlQty2U3ye4ftSN/KEvbbPrdRUp2BQt29JgGjVXUM8BngCofTgTJV3Q3YG7hIREY42+4ALgRGOo/2NeS1khd98iFMZhj5JF/ZweNGDeKyI3fKy3unsuuQ3rz808O45Yw9cvaeram8G8Tu2/Thw18flZfaUh2dQvUIf1FV3UyjdwDXoKpAlYgUAxVAA7BBRAYDvVR1itMb/B7gpPaed2vYyScfItJ5lSSjE9CasuyGEZaO4Ag/H3jIef4ocCKwFKgEfqKqa0SkGvDGutUAgXqriFxIXCth+PD8NyVpKWaeMvJJrkJUuyIXHrJ9ojWv0TryJjREZDLgF0I0QVWfcvaZADQB9zvb9gWiwBCgL/CG8z5+S/PApABVnQhMBKiuru5wyQMmM4x84jVP/erYUcRU+d2znwBwYwfM72lPfnmMf3MmIzx5ExqqGtzFBBCR84DjgHHanBV2NvC8qjYCK0TkLaAaeINmExbO8+TmuJ2I9uiuZXRfvD0nLjh4e7Y0RhNCo6TYrj2jbRQqemo8cDlwgqp6wzkWA4dLnCpgLPCJqi4FNorIWCdq6lzgqXafeI4woWHkk9TAQq/mUei8AqPzUyhDye1AT2CSiHwkIn9zxv8C9ABmAe8Bd6uq2wT7EuAuYC4wD3iufafcelJbf5oj3GhPvELEhIbRVgriCFfVHQPGNxEPu/XbNg3olKnU71w5ji2NMcb+/iXANA0j/wzrW5FoPuTFrj2jrXSE6KkuT5+UsgudOC/RNtiQmwAABotJREFU6CS8eXlyocsjdxnEpNnLTdMw2owJDcPoBlx17C4IcOhO3a9/u5FbTGi0I8/+8GDemb+60NMwuiHD+1Uy8dzqQk/D6AKY0GhHdhnSi12G9Cr0NAzDMFqNpZkZhmEYoTGhYRiGYYTGhIZhGIYRGhMahmEYRmhMaBiGYRihMaFhGIZhhMaEhmEYhhEaExqGYRhGaKS5lUXXRERWAotaeXh/YFUOp9MZsHPuHnS3c+5u5wttP+dtVTWt7kyXFxptQUSmqWq3qr1g59w96G7n3N3OF/J3zmaeMgzDMEJjQsMwDMMIjQmNzEws9AQKgJ1z96C7nXN3O1/I0zmbT8MwDMMIjWkahmEYRmhMaBiGYRihMaHhg4iMF5FPRWSuiFxR6PnkChHZRkReEZE5IvKxiPzIGd9KRCaJyOfO376eY650vodPReRrhZt92xCRiIh8KCJPO6+79DmLSB8ReVREPnH+3/t35XMWkZ841/QsEXlQRMq74vmKyD9FZIWIzPKMtfg8RWRvEZnpbLtVpAXN41XVHp4HEAHmAdsDpcB0YJdCzytH5zYY2Mt53hP4DNgFuBG4whm/ArjBeb6Lc/5lwHbO9xIp9Hm08twvAx4AnnZed+lzBv4NXOA8LwX6dNVzBoYCC4AK5/XDwLe64vkChwB7AbM8Yy0+T2AqsD8gwHPA0WHnYJpGOvsCc1V1vqo2AP8BTizwnHKCqi5V1Q+c5xuBOcR/cCcSv8ng/D3JeX4i8B9VrVfVBcBc4t9Pp0JEhgHHAnd5hrvsOYtIL+I3l38AqGqDqq6jC58z8dbVFSJSDFQCX9IFz1dVXwfWpAy36DxFZDDQS1WnaFyC3OM5JismNNIZCizxvK5xxroUIjIC2BN4FxikqkshLliAgc5uXeW7+DPwCyDmGevK57w9sBK42zHJ3SUiVXTRc1bVL4CbgcXAUmC9qr5IFz1fH1p6nkOd56njoTChkY6fba9LxSWLSA/gMeDHqroh064+Y53quxCR44AVqvp+2EN8xjrVORNfde8F3KGqewKbiZstgujU5+zY8E8kboIZAlSJyDmZDvEZ6zTn2wKCzrNN529CI50aYBvP62HEVd0ugYiUEBcY96vq487wckdlxfm7whnvCt/FgcAJIrKQuKnxcBG5j659zjVAjaq+67x+lLgQ6arnfASwQFVXqmoj8DhwAF33fFNp6XnWOM9Tx0NhQiOd94CRIrKdiJQCZwL/LfCccoITIfEPYI6q3uLZ9F/gPOf5ecBTnvEzRaRMRLYDRhJ3oHUaVPVKVR2mqiOI/y9fVtVz6NrnvAxYIiI7O0PjgNl03XNeDIwVkUrnGh9H3F/XVc83lRadp2PC2igiY53v61zPMdkpdDRAR3wAxxCPLJoHTCj0fHJ4XgcRV0NnAB85j2OAfsBLwOfO3608x0xwvodPaUGERUd8AIfRHD3Vpc8Z2AOY5vyvnwT6duVzBq4BPgFmAfcSjxjqcucLPEjcb9NIXGP4TmvOE6h2vqt5wO041UHCPKyMiGEYhhEaM08ZhmEYoTGhYRiGYYTGhIZhGIYRGhMahmEYRmhMaBiGYRihMaFhGDlCRCY4lVZniMhHIrKfiPxYRCoLPTfDyBUWcmsYOUBE9gduAQ5T1XoR6U+8uuzbQLWqriroBA0jR5imYRi5YTCwSlXrARwhcRrxWkiviMgrACJylIhMEZEPROQRpw4YIrJQRG4QkanOY0dn/HSnR8R0EXm9MKdmGM2YpmEYOcC5+b9JvCz3ZOAhVX3NqXlVraqrHO3jceKZuZtF5HKgTFV/6+z3d1W9XkTOBc5Q1eNEZCYwXlW/EJE+Gi9xbhgFwzQNw8gBqroJ2Bu4kHhZ8odE5Fspu40l3hjnLRH5iHidoG092x/0/N3fef4W8C8R+S7xBmGGUVCKCz0Bw+gqqGoUeBV41dEQzkvZRYBJqnpW0FukPlfVi0VkP+JNpD4SkT1UdXVuZ24Y4TFNwzBygIjsLCIjPUN7AIuAjcRb6wK8Axzo8VdUishOnmO+7vk7xdlnB1V9V1V/DawiudS1YbQ7pmkYRm7oAdwmIn2AJuKtNS8EzgKeE5GlqvpVx2T1oIiUOcf9inhFZYAyEXmX+GLO1UZucoSREK9gOr1dzsYwAjBHuGF0ALwO80LPxTAyYeYpwzAMIzSmaRiGYRihMU3DMAzDCI0JDcMwDCM0JjQMwzCM0JjQMAzDMEJjQsMwDMMIzf8D4ObJoNY2aCsAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "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": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO2dd5wU5f3HP9/dvU49OBCpghQBBfVEEAuCBUvEFA0aE4xGUoi/RE0MaOIvRaMxMbFFjeWnpKHYYjcKigVQpIoISO/tQNrBld39/v6YeWafmZ2Znd3bvTtuv+/X6163++zM7DOzM8/3+daHmBmCIAiCAAChpu6AIAiC0HwQoSAIgiBYiFAQBEEQLEQoCIIgCBYiFARBEAQLEQqCIAiCRc6EAhH1J6LF2t9+IvopEZUT0dtEtMr8317bZwoRrSailUR0fq76JgiCILhDjZGnQERhAFsAnApgEoA9zHwXEU0G0J6Zf0FEAwFMAzAMwNEAZgDox8yxnHdQEARBANB45qMxANYw8wYA4wBMNdunArjUfD0OwNPMXMvM6wCshiEgBEEQhEYi0kjfMx6GFgAAnZl5GwAw8zYi6mS2dwXwkbbPZrPNk44dO3KvXr2y3FVBEISWzYIFC6qYucLts5wLBSIqBHAJgCmpNnVpS7JtEdFEABMBoEePHpg/f36D+ygIgpBPENEGr88aw3x0AYCFzLzDfL+DiLqYHesCYKfZvhlAd22/bgC2Og/GzI8ycyUzV1ZUuAo6QRAEIUMaQyhcgYTpCABeBjDBfD0BwEta+3giKiKiYwD0BTCvEfonCIIgmOTUfEREpQDOBfB9rfkuANOJ6FoAGwFcBgDMvIyIpgP4HEAUwCSJPBIEQWhccioUmPkQgA6Ott0wopHctr8DwB257JMgCILgjWQ0C4IgCBYiFARBEAQLEQqCIAiCRd4KhfVV1fhwVVVTd0MQBKFZ0VgZzc2OUX+aBQBYf9dFTdsRQRCEZkTeagqCIAhCMiIUBEEQBAsRCoIgCIKFCAVBEATBIu+FQmMsMiQIgnCkkPdCIRoXoSAIgqAQoRAToSAIgqDIe6FQH483dRcEQRCaDXkvFERTEARBSCBCQTQFQRAECxEKoikIgiBYiFAQoSAIgmCR90JBHM2CIAgJ8l4oiKYgCIKQQISCaAqCIAgWIhREUxAEQbDIe6EQk9pHgiAIFnkvFKQgniAIQoK8FwoxcSkIgiBY5L1QiIumIAiCYCFCQUpnC4IgWIhQEJkgCIJgkfdCQaKPBEEQEuS9UBDzkSAIQgIRCqIpCIIgWOSlUIhqcagx0RQEQRAs8lIo1EYTQkFkgiAIQgIRCmI+EgRBsMipUCCidkT0HBGtIKLlRDSCiMqJ6G0iWmX+b69tP4WIVhPRSiI6P1f9qqmPWa9FKAiCICTItaZwH4A3mXkAgCEAlgOYDGAmM/cFMNN8DyIaCGA8gEEAxgJ4iIjCueiUrimIT0EQBCFBzoQCEbUBcCaAJwCAmeuYeS+AcQCmmptNBXCp+XocgKeZuZaZ1wFYDWBYLvpWGxVNQRAEwY1cagq9AewC8CQRLSKix4moDEBnZt4GAOb/Tub2XQFs0vbfbLZlnZp6zacgBfEEQRAscikUIgBOAvAwM58IoBqmqcgDcmlLmsYT0UQimk9E83ft2pVRx/pUlOFPlw0BIBnNgiAIOrkUCpsBbGbmj833z8EQEjuIqAsAmP93att31/bvBmCr86DM/CgzVzJzZUVFRUYda11cgNP6dFDHy+gYgiAILZGcCQVm3g5gExH1N5vGAPgcwMsAJphtEwC8ZL5+GcB4IioiomMA9AUwL1f9C5GhmMh6CoIgCAkiOT7+9QD+RUSFANYC+C4MQTSdiK4FsBHAZQDAzMuIaDoMwREFMImZY+6HbTghUxyKo1kQBCFBToUCMy8GUOny0RiP7e8AcEcu+6RQmoIIBUEQhAR5mdEMAGHLfCRCQRAEQZG3QiGhKTRxRwRBEJoR+SsUlE9BpIIgCIJF3gqFcEh8CoIgCE7yVihYIakiFARBECzyXiiITBAEQUiQx0LB+C/RR4IgCAnyVigon4IIBUEQhAR5KxSICETiaBYEQdDJW6EAAMWRsG0VNkEQhHwnr4VCWVEEB2tFKAiCICjyXCiEcagu2tTdEARBaDbkt1AojKC6VoSCIAiCIq+FQmEkhPdXVTV1NwRBEJoNeS0UFm/ai7poHDsP1DR1VwRBEJoFeS0ULjq+CwDgYI2YkARBEIB8FwonGEKhPia5CoIgCECeC4XCsHH6dVFZqFkQBAHId6EQMYVCTHIVBEEQABEKAIBa0RQEQRAAiFAAIOYjQRAERX4LBfEpCIIg2Ih4fUBE9wfYfz8z/zKL/WlUEj4FEQqCIAiAj1AAMA7AbSn2nwzgyBUKpqZQL0JBEAQBgL9Q+AszT/XbmYjaZ7k/jYr4FARBEOz4+RSe8/qAiL4CAMx8b9Z71IhI9JEgCIIdP6Ewk4h6ORuJ6BoAR7QwULQqMhSlg1IpVRAEAYC/ULgBwNtE1Fc1ENEUs/2sXHesMSiKhFAYDmH/YUMobNl7GJv2HGriXgmCIDQdnj4FZn6diGoBvEFElwL4HoBTAJzJzF82VgdzCRGhTUkB9h2uBwCMvOsdAMD6uy5qym4JgiA0Gb55Csw8E8DVAGYB6A1gTEsRCIqiSAjT5m3E1r2Hm7orgiAITY5fnsIBAAyAABQBGANgJxERAGbmNo3TxdyyxRQG97z1RRP3RBAEoenxMx+1bsyONDWti/2icwVBEPKDnJa5IKL1RLSUiBYT0XyzrZyI3iaiVeb/9tr2U4hoNRGtJKLzc9k3RUGYAABtRCgIgiB4CwUiWphq5yDbADibmYcyc6X5fjKAmczcF8BM8z2IaCCA8QAGARgL4CEiCgc4foN49genAQBKi0QoCIIg+I2ExxHRpz6fE4C2GXznOACjzNdTYTixf2G2P83MtQDWEdFqAMMAzM3gOwJzfFfjFGrqZU0FQRAEP6EwIMD+qUZSBvAWETGAvzHzowA6M/M2AGDmbUTUydy2K4CPtH03m205JRwiFIQJNfWS1SwIguDnaN6QheOPZOat5sD/NhGt8NmW3LqRtBHRRAATAaBHjx5Z6CJQHAmLpiAIgoAcO5qZeav5fyeAF2GYg3YQURcAMP/vNDffDKC7tns3AFtdjvkoM1cyc2VFRUVW+llUEEZtVISCIAhCzoQCEZURUWv1GsB5AD4D8DKACeZmEwC8ZL5+GcB4IioiomMA9AUwL1f90ykpDIn5SBAEAf4+BQsi6gmgLzPPIKISABFmPpBit84AXjRy3RAB8G9mfpOIPgEwnYiuBbARwGUAwMzLiGg6gM8BRAFMYuZGmb6L+UgQBMEgpVAgoutg2PDLAfSBYdZ5BEaGsyfMvBbAEJf23V77MvMdAO5I2essU1wgQkEQBAEIZj6aBGAkgP0AwMyrAHTy3eMIo7hAzEeCIAhAMKFQy8x16g0RReASFXQkU1wQRo04mgVBEAIJhfeI6BYAJUR0LoBnAbyS2241LkWRsGgKgiAICCYUfgFgF4ClAL4P4HUAv8xlpxqb4oIQasWnIAiC4O9oJqIQgE+ZeTCAxxqnS42POJoFQRAMUi2yEwewhIiykzrcTCkuCKEmKuYjQRCEIHkKXQAsI6J5AKpVIzNfkrNeNTKSpyAIgmAQRCj8Jue9aGKc5iNmhpl0JwiCkFekFArM/F5jdKQpKS4IIa4F2cYZCItMEAQhDwmS0azWagaAQgAFAKpbyhrNgKEp6MTijHBIpIIgCPlHEE3BtlYzEV0Ko9ppi6HIIRTi3KJy8wRBEAKTdpVUZv4PgNE56EuTURyxX4ZYXISCIAj5SRDz0de0tyEAlWiBZS50RFMQBCFfCRJ99BXtdRTAehjrKbcYkoSCpCwIgpCnBBEKjzPzbL2BiEYisWLaEU9xgcN8JJqCIAh5ShCfwgMB245YWhcX2N6LT0EQhHzFU1MgohEATgNQQUQ3ah+1ARB23+vIpH2pXSiIT0EQhHzFz3xUCKCVuY0elrofwDdy2anGpn1Zoe29aAqCIOQrnkLBzGR+j4ieYuYNjdinRqd1kf0yiFAQBCFfCeJoPkREfwQwCECxamTmFpOrQEToXl6CTXsOAxDzkSAI+UsQR/O/AKwAcAyM4njrAXySwz41Cd8Z3st6LZqCIAj5ShCh0IGZnwBQz8zvMfM1AIbnuF+NTkSrgCcyQRCEfCWI+aje/L+NiC4CsBVAt9x1qWmIhHShIFJBEIT8JIhQuJ2I2gK4CUZ+QhsAN+S0V01AXSwhCMR8JAhCvpJqjeYwgL7M/CqAfQDObpReNQGtixOXQoSCIAj5Sqo1mmMAWsyym358/aRu+Pn5/QGI+UgQhPwliPloDhE9COAZ2NdoXpizXjUB4RBhYBdj3SDRFARByFeCCIXTzP+/1doYLWxNBQAImc5m0RQEQchXgqy81mL9CE7CZAiFmJTOFgQhT0mZp0BEnYnoCSJ6w3w/kIiuzX3XGp+QeTVEUxAEIV8Jkrz2FID/AjjafP8FgJ/mqkNNScjUFOasrmringiCIDQNQYRCR2aeDiAOAMwcBRDLaa+aiLDpU7j/ndVN3BNBEISmIYhQqCaiDjDXZSai4TByFlocSlMQBEHIV4IIhRsBvAygDxHNBvB3ANcH/QIiChPRIiJ61XxfTkRvE9Eq8397bdspRLSaiFYS0flpnkuDCWulLqpro4399YIgCE1OSqFg5iOcBSM09fsABjHzp2l8x08ALNfeTwYwk5n7AphpvgcRDQQwHkaJ7rEAHjIzqhuNsKYp3P/Oqsb8akEQhGZBkOijYgD/A+B3MEpnTzLbUkJE3QBcBOBxrXkcgKnm66kALtXan2bmWmZeB2A1gGFBvidbhLSrUR+VCCRBEPKPIOajv8OYvT8A4EEAAwH8I+Dx7wVwM0wntUlnZt4GAOb/TmZ7VwCbtO02m202iGgiEc0novm7du0K2I1g6OajgrD4FwRByD+CZDT3Z+Yh2vt3iWhJqp2I6GIAO5l5ARGNCvA9bqNw0nSdmR8F8CgAVFZW5mw6HxGhIAhCHhJEKCwiouHM/BEAENGpAGYH2G8kgEuI6EIYy3i2IaJ/AthBRF2YeRsRdQGw09x+M4Du2v7dYKzd0GhEtfLZkVAQJUoQBKFlEWTkOxVGUbz1RLQewFwAZxHRUiLydDgz8xRm7sbMvWA4kN9h5qtgRDJNMDebAOAl8/XLAMYTURERHQOgL4B5mZxUptRp9S3EfCQIQj4SRFMYm+XvvAvAdLNUxkYAlwEAMy8joukAPgcQBTDJLN3daOiaQlg0BUEQ8hDiAHV+zFyC7tCESHMonV1ZWcnz58/P2vEO18Vw3G1vWu/PG9gZj36nMmvHF9LjQE09DtREcXS7kqbuiiC0KIhoATO7Dm4pNQUi+h2AqwGsQcLx2yJLZ5cUhnH310/Azc8bVrG3Pt/RxD3Kby55cDbWVVVj/V0XNXVXBCFvCGIjuRxAH2Yexcxnm38tTiAoLhl6dOqNWihXPvYRpryQTl5ibllXZazp9OKizXhlSaPGHAhC3hJEKHwGoF2uO9JciITy18E8Z81uTJu3KfWGjcwNzyzB9dMWNXU3BCEvCOJovhNGWOpnAGpVIzO3yLWbw3ksFITMuerxj/GtU3vgguO7NHVXBKFBBBEKUwH8AcBS2DOTWyTkqJTKzElt2eJwXQyfb9uHk3uW5+T4QuMQjzM+XF2FD1dXif9DOOIJIhSqmPn+nPekmVIbjaO4IDd1+W5+/lO8smQrPpoyBke1DVROKmfUyxqkGVMblWsntByC+BQWENGdRDSCiE5SfznvWTOhpj53qRJLN+8FAFTXNX2ZbikVnjm10Ra55pTgwcHaKPre+jreWRE8OrE2GsNLi7cgSApAUxNEUzjR/D9ca2uRIaluHK6P5czLrm6P5rC4z0ERChlTU29oCvkcpJBPrNtVjfoY4563vsDoAZ0D7fP715Zj6twNOKpNMU7t3SHHPWwYKYUCM5/dGB1pThzbqRV27K/BgZooDtflbhbYnCYNYgLJHKVNSpBCflAQMX7ndEyun6z/EgAQOgLukSDrKXQmoieI6A3z/UCzREWLZcaNZ+GP3zAKw6pZYC5gU1eIxTP7jj3Vdfj5s0uyIrjqciQUXly0GQdq6rNyrFz1saEogSqaQn6gimXWx4LP6r48VAfAXkqnuRLEp/AUgP8CUFldXwD4aa461FwoKTScy4dz6FNQmkI6N5fOPW+txLMLNuOFRZsb3JfPtiSW3Y7Fs3PjfrZlH254Zgl++Z/PsnK8umbqDFeaQiQs9bJyxSPvrcGaXQebuhsAYPkFgk5S6qJxSyis3nkgZ/3KFp53MREp01JHZp4OMxyVmaMAWrxnrcSMOMqlo1mRSeTPpj2H8MLCLQDsy4hmys+fS2QyZysSad9hQ0PYdaA2xZbBiDZToaA0hX2H6/FldV1WjvnBql1ZE85HOofrYrjrjRW4/JG5DT5WLM74w5srsPtg5vdk1Pxdgj4nP/rXAsvi8KuXlmHbvsMZf3dj4De1UWWrq4moA0y/KBENB7DPc68WghIK1zz1Sc6+oyGawlcfmmNpMdm2Uw741ZtYsOHLBh9HPTQFWZpBZ6pRpSIWZ9zz1kps31eT0f76xGH59v0N7s/bn+/At5+Yhydnr2vwsVoCMfNByUYwxPurduHhWWtw28vLMu9PmkJhxvKdtvfbMrzPGgu/p1WNNDfCWOugDxHNhrE85/W57lhTU1xgXJp0HLD3zvgC767cmdS+Ze9hvPnZNs/9MpkBV2kznWxoCk6+/vCcBh9DDeLZWpsimqHvJRWzV1fhgXdW43evfp7R/gdqEoNVNoIH1lUZZpL3vsjucrOp2HmgJuPw2tpoDBt2V2e5RwafrN8DwGUZxgxQNv3aBlgAEppCZj3SfYAbdlfj2Ftex+qdzcM0BvgLhQoiuhHAKAAvArgbwBsAHgNwTu671rRkkrB274xV+O6TyZrFFY9+hB/8cyGisThG3zML4/46G8yMLXsNNTKappngkCOvIRdLh7pF0ryzYgf+9N+VgY8RzbamEM2NprD5S+N3UH6kdNEFdDZMPtW1xqDxwaoq18+/rK7DTdOXJN0HDWXYHTPxw39mVhH/1hc/w1l/nJW1oAId65nKws/vvKvV7zVz+Y7AZh0VGJKpj0sXCi8v3oponPFiFvyC2cLvaQ0DaAWgNYAyGOGrYQClZluLpjCSPafhzgOGunigJoq1u6qxZNNe7NFsz+na8K947GPb+1zkObgNbtc8NR8Pvrs68DGsqJxsCYUcaQrKCdixVVHKbeNxRty8NvsO1ePVT7fa7NOxLKgKbsENVz72EZ6etxEA8Oe3v8DzCzfj+QXBBpLnFmzGP+au991GOU/fWZGs6QZh9mpDgO2vya6g+tRM8AQS0XrZYk91Hfrc8jr+MXc9rp06H197KJh2rLSNumg8o0mA+n2nvLAU97z9BYDmkauk8MtT2MbMv220njQzOrcxyk707dSqwccqK4ygpr7OcrwC9ocn3TC1JZv22t43o/vJxiFzRpQ181GWfAq7DtTinx9twE/G9EUoRFYUSZCJQOUdM1BSEMbsyaNx07NLMGP5DpzRt6P1eUPMEgq37PI5a3ZjzprdGD+shzWJCLo64M+eXQIAuKyyu6cGnKmC82V1HW5/bbml7R7OsvYyb90e63U2THPqWWFOaHgPz1oDILitXxcEj7y3BpPOPjatPmzccwgAMM0U8kDzEgpBfAp5y1n9KlBaFCTp258y8xi6UNBfp6spOAev5hqloswbby3LzmJF2YqK+vG/F+K+massp7ByFAfJF9lTXWeZ/baa/3Wt71AWckZ084KzLMJnW/ZZtux0zYb3vOVt+svUX/Pwe2vw/MLNVoRZtjWFMu35S9fM6oY6BiNxP23bn57jV+/HFzvSDzFdtDE5iKM5JT76CYUxjdaLZko4RBknlumUmrZqXRDsPaSZj9K82TuWFdreZ2MG3bakwPb+7P4Vae2/aOOX1iB5+h/ewc3PLbEGyEyiRtxqxGRjUABgxbsXRYzfRfWvNs1ERaUB7dfs6IfrY9i05xA+3bwXsTjjqdnr0g5r1mthrdp5EL0mv2a9v/iBD60BPF0NzG8mnOnEwtmHg1kWCqUOP4+f3f9gbRS/f325r7Nczy1QwjddDUS/VpkkVLoFrxwRQoGZ93h9li+EiBB0cur3UClHqz446AIi3eij9k6hEFBwMbPnAB0JEdqXJgRDaaG3hhR3OdevPjQHo/40C4DhuJ0+f7NtcHPbxw9984rWhq0/W5pC1UEzu9S8bspc8/iH6/Bn08YbBOUr2XMwIeBr6uM44+53ccmDs/Ha0m349Suf4y8zgh8TsGsbbn4DNQmIBDAfLduaiB73W+s6U6FQHLEP2geyLBQKHf4o5YR344GZq/Do+2sx/RP3haIO18VsK/hlGt0VbaBQcLuPm5H1KFBGc94SDgUbzL7x8Bz0ueV1188+WLULS81sYX1AtguF9B5Ip/0xaGjcPW99gcH/+1+rHzX1MUz690Js2F2ddKP6RVZ4OXydD8jf3ltrvU53lq8G7GtGHoN7LjNKjmS7vHc0xqiLxvGfxYmB4v6Zqzy3dw7QapZcrQ3iumZZY7ZXHUgvoU0XCmpJUp1E/kfqkeTvczZYryt8HOnpCoXPtuzDwo1fJvkoDtZmN/rIqUXH4owf/nMB3vxse9K2agbu9Tz89tXPrXXX62NxPPCOPWgi6MCs549kEoHk1r9chJVniggFHyKhUMpokuraKOb7JHp9+4l51mvd9rz3kOZTSNNE5Rxgg2YMq8ghpeIv27oPr326Dd//xwLUx+yLCfkNwM6bWhecXvHWce06BgmlVINU5zZFlgnhf6YtsmLWs8HFD3yY1iz+6U822t47TW6A/dook4ASFPsO1WOFR3Lbpj2H8NYyY6DTfQpukUjq9w/iaO7QKqFV3vH6cs/t0hXaFz/wIb720Bwrn0eRbU3BqUXvqa7DG59txw/+uSDtY23+8pD12u3+DuLsPVwXw5w1u633mRSSXLDhyyQtZV1VNT7f2vDEx2wgQsGHUIhSzqDiaRgkVTw8YBcK6WoKsXgcZ/arwHs/HwUAuM9ndutGnWNGtfNAbdJD4icUnCUC9IHrnD+/57qPGnQm/WshBt72XyuE0cnCjV/ilSVbtYGPLDNN1cE6XOZT6mDXgVrsSNNp+NKiLYG3dQ6c+iU6oVtbsy2xzUumqULt99WHZmPsvR+4Hvv8e9/HxH8swLx1e/D5tsTg4JaroH6/VPfmZ1v24SEzsibRZ/d99PZ0av47gx6yLxTsfbnisY88t03V71SDfpDZunOSmMp85HU/OpNZn/5kEy683/3eaGxEKPgQptQPXjqTfBWKBhiO5tZmZEW6ZpFYnNG6OIJOrTNbra0uZgziypZeECZE44zhvRPLgvolil3pyJMIskiQuo6vLTUeBi+h8LWH5uD6aYsQiyWEQlAf3Cl3zMCpv5/p+flj769Fv1vfsLVtdXG+ejkznfeC7tAsioRAZJ/Zvm/OBtXAttbFFKRQJqPL/5a6vo9ybKeakLjZzHX/UyzOOO5Xb2LavI242Vb7yv+4um/MeetmWyj4adH700yU038/t8faTyYs3rTXCB5wXJtUz65eaFInaELnX97+wvNZyRUiFHxwagprdh3Evz+2mxDSCeXTZw1b9x1GR8uBmq6mwIiEKOOIBaXyKt/Cjv3GzH/Q0W2x/q6LcGa/Cl9bqQrJVPiV7lY1pJwDaiphq2bXkRDZBpqiBiQV3vH68kA24BF3vuNassE5a9VniZFQCJEQuUaSVddFfatjplviQCXbpbqGbveHfg6H62M4XB/D/760zCZAUpW6+IvmjH9t6VbbZ1n3KfjMxE/49Vuu7W6D+3tf7MKH2uCqBJtuAvR7ni7962xc8uDsJCGV6jfwqkvmFSTgrL9138xV+NbjH7tumytEKPgQcQiFSx74ELe8uNSmpqbjoNMdiIs27kVH096bbvRRNM4IE9kcjalU59P/8I71Wg1mzkiOMtN2XximtKIq/CJCxg4+CkDydUolCNeaYaPhUAjDjkloMG52/Fygm/oUThOLPkuMhMnwQbncD1UH67BpT+J4zuCF8Y96m0TcUI7rlELBZXTUTWBqwHUKylRriOimz9mrd9s+U0737ftq0o44cyNbYcjvO7QmJRS6aGujB/EpOK95Km3NyyTlFSQw/E5vTbexEKHgQzhE2L6/xqr3o254/cZI56bVM1Vro3F0KDM1hXTDNeOMcIhsjmE/ZzdgH+QSmoJ9Vteq2BhwiwvCgWLrZ63ciV6TX8Mqn1mwGsRjcbYJmlT5H980B8pIiFAQDuGmc/sBAE7ollgc9cVFm7NShG3yBQOSQh/dflddK4zG2baN4fsg15Df5dv24w9vrvA8drpmkMP1yfehG+6aQuIcvLSmVPWL/GpERWNxbP7yEIbfORP3v5Oer8uNTIrOxeKGwH78g7Woi8ZRG43hiQ/tFWfVNezXOVGxJ4ji7TQXpfMbfHxLIvUriPmoqUrFi1DwQc0cHnx3tS1iRn+oM9UUAKB9WQEiIcpIU1DZrOcONNaI/e0r3hU+nVpEnSUU7P1RUT6tiwsCZaY+Y8aDz12z23ObVqbfJMZsE4pT526wFZLzQj1UPx59LLq1L0FRgVr1Ko4bnlniaoPftOeQLeQ3FUWRUNKMz01oOZOWdCF3qC6GSIiSzIuKFdsTgvO/y7bbnPVuWlmn1qnrMGUiFIIkXs1f7z3BePvzHb6TgPoYWyYQ5+w8E9QgPO264bj29GN8t1VnVheL4+9zN+D215bj73PXu5aBP1xnHLdHeanVFqQEvdOEuLaq2ldL17UPVTon6Hc11RK5IhR80JdX1O2+mQsF+0DbqiiCiOnkTYdYnK2b7afn9AUAjDy2o+f2zuMnzEf2/qizbVMcCVTtUj0L1T4+BVWmIBZjzFhuL3dRefuMlN+hhB8RobQwbDn6lAnDbfA/4+53Mfbe91MeW1EUCbsIBfs2Nz6zGOt3JwIFamM9L18AACAASURBVKNx23U9WBMNXPjv+mmLkpz1TsoClFdJFS7tZtrQtVKvQcdPc7nu7/OTTEYAsOCX56BPRRnqY3FrcM5GPR+VTDi8d7ll3vRCnW5tfdwK066Nxl0z1VWNKj1kN0j0kfNZYkaSFqLz2Adr3T8IEOElQqEZov/+WzTziz6zT2dAd25aFAmjIBRKP/qI2RJYg45ui0iKCB3njFCZDZxCQXWvdXEEtY6ZMAAcaxYH1G38gH8RtFZFYavPCzfuTfp8VYraMfpsNxIKWSYc5WxtVeTuY3CWdPCzbxdGQkn1N52awguO0NW6aBz1sTg6mNnlB2ujtnyR687wn9WuTHHezvIObuj33gX3fYBTf28XsmpQufLUHlabHj3jpSlkYrIpLytE+9JCfLCqyppNZ7PIGxG5mlz0SYHSxKtro5b5s7ggjIdm2ZPURg/ohAPmvV+uVQc4WBvFo++vwSUPfujZD7fAEi9Nefu+moyrzgL2KK93G3CcdBGh4IP+kG/XIof0h6YhxegKIyEURJKFglp3wS1rEzAebD1xyfABeAuW6fPtaf9qMDjgEApqsG9t+hac2oJSk9UsS5Uy9isCZ2kK8TiOapMcQvvykq1JbTq6tlYQJuvaq0RAJXS8WG+q937mpKJIKGnilkpO10ZjqI/G0c00PzjvgzbFDXOIB1nPQxd0y7ftt6LIrD6a98Tt4wZbbfqg5uVT8PL3+GkQRGT5tf75kZlFnYZMmPLCUgz3CScGgAKXyLMhv0lEIKn+7a6uswbUkoIwPnGYw3T/kX6da6Nx/P71Ffh08z7c89ZK3PZS8tribjlFxR4CvKEyUdcUvvvUJ/j1y8tQefuMnK/IJ0LBB/0h0IWC/mD95OlFSfvN1Mwk7Uq9B4eiSMj0KdhvtGNvfQNrd1XjlheXuu4XjTP0SVNxQdg181XxG4e/QTlDdU3hN5cMstYTUAO5M6pIDXybvjxsqyGjBtzfXTrYtn0kRNZDt3TLPtd462Upsjh14RfWosFUmKo6vttM6pt/m4tRf5qFR99fi90+aye7hbmmCjXefziKrftqMLBLG/xoVB/c7jj3cAPLhROAWy88Lqm9q1a/SM1QvWzatdEYCsMhm/36+YWJUh2zXFYJBLy1368HXG9A5aIQgPnr9wRaD3navI22Z0zRsVURvjLkaAD2CYIb+7U1wdXzUFKY/Nvu0vrjFQX0wDur8fe5G5La3QRpJESYNm9jUl6Im/nur1eeBCCYhcEZGvzUnPWoOlib9Dxnm5wJBSIqJqJ5RLSEiJYR0W/M9nIiepuIVpn/22v7TCGi1US0kojOz1XfgqIPtHosvj6I6w5ExbVT52OfafMOE1k3tZPCSAgF4ZCnul4cCeFQXdSmScxbtweH62MOTSGE2voYYnHGOX9+z3fpTyCRZ6ALhYKw/XgAUOO4KdWNvKe6DtdPW2QJjd0H60AEfGtYD9v2xQVhy4RwwzNLMM+lREUq05k+EETCCa1K/TbK5/CfxcmZyR+btfhnr9ltKzECAEO6J6KYilxm5Wt3VXtqakBiMrBy+37cPHYAzh7QyfZ5gfn7+C3c89yCzVi8KdmkBhi/x6m9y5Pa9YnKm8u2Y/XOA0k27Tc/245ek1/Dxj2HkgTeX99NZDjfO8M9Osgrw35VmvkUISJ845G5aYXcxuOMWSt3YvSfZuEZs6xI62JjkkIeU28lFPfpQsF8XosiYduM/frRx1rVfIHga1Io3Na6eGnxVkx5YSkm/F+ipA0zY/LzyZO6003fXxB/weuf+j/HuSKXmkItgNHMPATAUABjiWg4gMkAZjJzXwAzzfcgooEAxgMYBGAsgIeIKLP1EbOELgj0WUMQH8DMFYa2EI0zOpQV4sx+iVLUapwrjIRMR3PiePqsb+u+Ggy87b+45qnEEp8q2kYfLEtMTcFIkjqInz2byE51Ul5WiFU7DuBwXcw2UOozJlX50hmW6pWAtn1/DTq2KkqKqCguCKdMsEuVD6GXUSjQnPKqb8o56KzWqVNbH0uarU4Y0RMn9zTmI707luHhb51k+/y+mats9XWcIasqO9lLoF81vCfu/sYJ+NapPVw/B4zFby7962zXz+65fIhV2lvH6Rv5z6KtSTbtZ01z4Udr96A0hXnNjWzlBqhZdTrC5IVFW3D1k59gbVU1fvH8UkTjcVcNYdzQxERLDbBKU1i544BV+I4ZOKpNMYZ0a4sHrjgRN53X3xY2HAlRoEgvRdDS4Ot3H7JpDn++3CjqqGSQn3Nacb+jYF9jkTOhwAbqbigw/xjAOABTzfapAC41X48D8DQz1zLzOgCrAQzLVf+C4BWrrx4aP+flenPQiJk5BTef39/6TA2URZFwkvnIbQbxwaoqm2AA7CFtKq9AyRM/08fZ/Tvh/S+qcMbd72DNrkSMvz4LUyYZp5/COQDqwlE9WP+4dhgGHNXaPE7I9YHWZ8+pBKyuwYRDoSShEAoZuQGfbPAulFcXi2OP6ZhWfo1IOISnJw7Hgl+eg+7lpbjg+C548runeB5D1TZSIcAKr+J+JYVhXF7Z3VaOPB2OblfiYday/wYPvrsaMx2mMzWzrjpYa5VAf0o7t1SJjtlYQwSAayhoKtQqcVZfYmxl/6o76bQ+HTBa08xqTC15p0thyJiZZDi4a1tLY9dn++EQ4a0bzsTxXdsG6l/QtUHWO0qanNKr3Po+Jw3J0s8FOe0NEYWJaDGAnQDeZuaPAXRm5m0AYP5Xv25XALpHdLPZ1mR4OVDVQOZXl0WF/qmZjj7jVQOwMh/pdkovtdIZxeCmKcQsYeXZLfQ/qhXqYnFrTYHE9ybOVZmPLv/bXEyds95qdwob3ZSh4r3P6FthLU9ZUhB2jcc+VYteUufOzK5rB+gaTIGW06FrCj+Ztghrd1Wjsmf7pP2NbeOW1qd8PAVmUlwHTUAV+YSUti8rRO+OZfjBWX2Sju2Hm3M0KG4rqwXRUttoWd+qzIguXJVgGaqZ0PpUlFmvD9REk0qZ+HHz2P6+n2dajkVphuo6qMMc26mVTXOrqY9jx/4aRONsCW9F3Eya1M/fOZdrV1qInh1K4YWuIc/SZv+v/Ph0z32c/oDEOdivxa+/MhCzJ4+23qvlfzNZ0S1b5FQoMHOMmYcC6AZgGBEN9tnc7c5JmtIQ0UQimk9E83ftanhyjB9es4ft+2pw34xVvqaPh2etwQerdqGmPo5QiGyDuDJ5FEVCKC0M41BdFBt3H0JNfSxweQl9sldUEEJNfdwaMP00hc4uEUCAfdUxPSLjf19eZmlEzjo0+szsHlM9BhIL9BRG3DUFHVV475gpr+MmxywRsA9mkTBhx/5aPD1vo5V8NH/Dl5Zf59sjerqfWzRmaTlqFu2WU+C2RnNM0wpLi8JJszpdU7jt4oHJ/U/TZq3jNpg+8/0R+PZw9/NUtC9NhFmq0Fb9WDX1MezYX2PzZ5zZrwLr77oI5WWFePqTTRh5V6IsSioGHe0/yw4iFNw2KQyHLE0bSEymmO2abU19DLvNSc65x9k1OSO7mW2/rT4pUZMwvyhCXRC/ptn5j+/mfd7OyZ3SdpzXolv7UnRsVYS5U0Zj9IBOltCY+Pf5nsfO9dILjaK3MPNeALNg+Ap2EFEXADD/qynwZgDdtd26AUiKV2TmR5m5kpkrKyrSWzIyXe674kS8/OORSe2T/r0Qf5nxhaeTUKHWUvhwVZXDDJLQFLbvq8Hs1btx5h/fxZWPfZSyGJlCn0WWmOYjy6zluL914VbhYT+tcdEUFHtNW62zHIdeA6dEEyQqiznO7rHqKpQ1RMYDlyqHQBEJhVB1sBaTX1iKlTsSUUuH62M4qk2xVWfJSXVt1CozofrmZkJxEwpKwKp6U84SD7//6vHW62tcMm713+mGc/q59s+LLm1LcPfXT8AzE4dbbaf0Kk+K8tJZs+ugTair/uoaV200josfsMfiK/+FPmgpgTdnTRX+9p69BLfO6T6Jk4ChlaWizGWlv+q6GOpiCZ+Cfivprw/XxyytfbBjIrd9fw1qojHb+f/f1adYJk6lcapn5zyHeRBIr4rx0s378L2p85Mi99T3OxPklCbZpW0JigtCWLnjAG6avsTXr5PN/A/X4+fqwERUQUTtzNclAM4BsALAywAmmJtNAPCS+fplAOOJqIiIjgHQF8A8NCGtiiK2WjsKZTIIqhbvPlhrNx+Z/4vCIVvp5oUb9/pGJby+NDFL0U0Eyqeg+yZ0J7L+3f06t3advZdrs0ung1Pla/g9HPrMTYW0xjXV343CSAi10ThunL7YcxunpqDQH7o91XXo3KYoyRms0DUapcW4XWdfTYGNGWvP8lLrAb95bH9ccHwXz77r/S8vK8SJPZLvJZ2u7Urwf1dX4okJlVbb5ad0x6m9OyRte81I9+S4W19cakuuVIOtHmWzfV9N0sJMSgPSB3A1+77ysY9x5xsr4OTs/hX4ypCjUz4HQTK9/TK4ncdnhwGhRjOdOhPc/vjflWC2t5cVRawKAGoyoyYmbnWd0llj/H+eXoQZy3ckmX/UNXCaU/XrrQb75xdu9r2muV6lLZeaQhcA7xLRpwA+geFTeBXAXQDOJaJVAM4134OZlwGYDuBzAG8CmMTM6a143sj86yP3OjdOQiF7NqYKp3SutQwkLx7fWntYfvSvhdZrPRFM+RR0s9Evnk9EIOkmqY6titCtfSLWffwp3fHIVSfh8sqEkuZ8MC7/21zMWVNlqu0+J2pSpmUxu81qlNApKQijLmZfDtOJPsNTph8gea2AVsURn5BF7btNLchVKLgMXsrsFI2xteBP2xLjdysNkGBWoNmS/QQkYAzMowd0xpjjkmerTpwOb0VxQdim0anBVs+QfuXT5OutroueX+EM43Xy5HeH4YErTrTef/VEdxdgKhMi4F9kz9IUzPdOJc8wnRqN4RBhxo1nJR3DKSxuHtsf940favm/1My8xPGbMjO+rA5eR0tpV/oyqt87/RhLQ03qVyTZggAAnc21Uq4Y1t1m7gISEUy5IpfRR58y84nMfAIzD2bm35rtu5l5DDP3Nf/v0fa5g5n7MHN/Zn7D++hNg9MZ9drSYHHEj0+otA1u6gZ0i2F3Jse4zaA6tirEOccloi+KTZ+Cbhfdr2XwOk1SrbTBtSAcwtjBXWwzmA5lhTbH477D9dZi6EEGQjU7jcfZdUD41cUD8f0ze6NnhzLXpUQfuSoRHqoP1O1KkoWo8ztToQSSm5nOT1PQbdvK1VXq8p2PXHUy/jMpYXJM2JKB1h4lOfy+3wuvpKvSwrBNUygvM76zgzYB2bY3OUlMXWe9zv+6qmr0mvya9T5V2XK1lrYTL2G4ac8hKwDAz6avnhcl9ONsP/8JT86zJkSRMLlGfDkFflEkjHFDu2rHNL7DmUleF4tj7+Hga2yr4BR90nK+h1kTgKuvEUjkEv123GA88/0RuPvrJ7hulwuaVyxUM+febw5Ne5+j2hRjwFFtXOu2lJcVWjMVhQqXUyF3HVsnD4TfGdHLHkJaaGgKbusDA8m5APrg5PbAEhF+4rB/q4clSAi7Mm1V10WtGahOeVkhplx4nKdPpq02+OvXzZkdrs9+3RLQ3LhmZC+UlxXiHJfZuJumsHXvYazZdRDz1u+xBkw1U3Wb3Y4dfJQtqkf1Mc6waWhupBOa6GVe6NS62KYFKW1Udz67aQDqu/Xj/vQZu1mvujaKbw/viTbF7gI43QVlzrj7XSvU2s80qc4nccszzurXCSNMs1pdNG4VrAyHyDbpUXgJUcVN5/VH745lST6J+hjbfGeA4aj+4OazXY/jltzmpynp97d+/ZRQKHAxO4Uo9TLBDUGEQhoc1dY9cucKRyavjkoe0n/8pycOx/Wjj0U4RHjkqpNt26tVln56Tl+s/f2FrrNj54BQHAmjLhq3reymb+MsgaEPpl513Z3OQVWAzq+chuK4LoYTj0CuCVip0Puuq9fOmaoeMqxmevq+rV20rL6dW2Phr851jcJym6lf/MCH+IrplFWzUfVdQYrWKQFZUx/zLXkCJPtydJ79wQj85ZuJmbjX7/bUnPW2xChVx0ofVNxKlqvv9hvAoqa29N7Pz/YcFN1wm3goR//ctbuxYvt+36g79RlpAYrhEOHWixJlQFTph0jI/Z5LFRo8tHs7vPOzUUlBFvPX70lKWOtd0Qrdy91DWN3Gaj//gC0AxUcD+MqQhO/qQG0UfW55Hf9JY33xdBChEICz+xtRTvpsS+fMvt7RF2r2GQ4R2pUW4NdfGYjhvTvgpvOM2O6yoojrDLFtSYERyuryQDlvMjU4fVdLcNNt+UpFV036Mb0GAeeg45wtKU7o1hb3jbdrUK2LC/Dny4fgH9cO8539Xn1aL9d2e2VU3afgM6iaD6N6sP502ZC0Q/e8zDdK+Ow7bAwOapbmZj5yogRZbTTu6fNI9f2AEXX01RO7We+DBDmMPLaDLTJI2f+dOSpAwqSYyu8RCRHalxV6Dope+zjRtdqx935gJRe6oYTCAHOyoQo3ul0Dr+uia29+XDjYHjjwxY4DqHWadANMBnT8JvV6d/3WWCiKhJPCdu+f2fBFjNwIZojNc5787jAws+dD7XdD6w/64tvOc90mEiLUwhAgyqegBhM31ds5oyh1mRHrN9Dh+hh+OKoPfjF2QNIxvR4i5+DgJRTOG9gZ44YmOxi/dpIxgPk5K380qg+e0pLjAOAHZ/WxawqacHIz1/zlm0NwwzNLEsXHzF29TBx+eEUvKfaZv7MSEnotfi9UtdQg+Sdezkg3UplDAOBf3xtue9+nwkiM2l1di67tSmwJauq7vX5nRdA1I3TctBpnLo1fkrW6dif1aI95t4xBJ5WV7nLvepmqurUPJsScA/Pu6jpUmM+bqtLrOzlxwW8RLd2vuGnPIc/tgGThMqJPclRaNhBNISB+s7xWRRHPAeXTzftSHls9aBccn3BIqRvP7eF3DuRuMxfV3/pYHPUxtkVVBJllOs9HmY9+OKqP7ftSmTb9NAW3mfFVw3vYcgj0vjojQ4zjG23OQaUgHEo5M3eSasD70hww/QIFnKgZeCp/AgDc/lW/3E476RZyAxIRRsyGKXTqNcOsdcJVP9WyrV7nFiSSKLmvyfu8uyJ44qkefNFJM/u5/V7qu57/4Wm29lQC34vdB+ssn4bSDFsHmHAURhJri3vVx1p354W2CMRUZUGcmeM/P98/kzxTRCg0kO+dfgwuGXI0Zt50VpIZJShq4Ndr8KubW934Kv1d/0zhFqGkTBzKB6Dbv1OZCIBkG6yaQR7XpQ3+9u1EHH2qBaR0odCqKILObRKDjZtQKAwnFrwZ4lD53Wz4ahZqObQ5+XuzhbqWSta0SxGNAxi/1aPfPhnTvz8i6TM9yKBf51aBhIwik8FZvyZtiiM4q18FbjzXGFh6dyyz+gt4m0QzKVnhNnhP+vfCpLbu5e6C00vLcrsG6lka3LWNow+ZRey8vGQrVpoZ82pSogsFvTCfTvvSAtx8fn8c1aY4qS8K56Ql1bX90ahjbc9AuhpLUEQoNJBzB3YGEaF7eamrGSUIF5kJUG71V9R9MqBL4sZK1hTcwlaNAWadWfROD7ULMqDo25zWp4M1WysMk61QoDORKOk4Yd18di5m/yJR50WfvbVxKT/h7KUzXHD0gE44u38FJp7Z2yoxofrjJnC8Bh2dWT8bha95xNsrXvnx6fjtuEGB1tkFgPMGHYWj2yV/95l9Exn591yW3oQik0FOv36dzDj4K0/tgfV3XYR2pr9siFm64daLjrMqe9q+NwOhECSjGTACNtbfdVFSu1cWvp9PwWlizaTfpx/bEXXROF5eshWF4ZB1zfXE0fvGn+i6b1EkjMpe5fjoljGBB+8gAld/9jKtKZUKEQoNJJ3Yci9+fckgLL7tXGu5S52XzMSuJVr4pvMGdyuPrOLwp7xg1HTXbftBTA+6HXi35piMhEI4SSs8F2Cp2cS+4ZBt0LcJAEokKFnuAcc97/QpPPadSkTCIdxy4XFWYTu1r9vvcssFyYvWOOnVsQwn93IvrKcY3LUtvjOiV8pj+TGkW1sMN0MqSwrCvnV0gtClbXFKn4SuKXRq4z7QPvadSjxy1Uno0KrI8gvpBIk+A2CLTgpalsHNxHP/FSdi8gUDXLd3E4xeNYbSNSUC9glPUSRkrdfgJaR0nFFMQQgiuHIYiWohQiFNLjrBHp2QDaFARGhXWui6BKNSTxmMr5sPqfMhcxsMXv10Gypvn2HVeNGzTb95SiJ7uavLDBawCwU916AgEkJ5WSF+Mqav2a/s8OCVJ+K0Ph3MWZiZrOTYxulTcJspKYdzYSSEP37jhEBho04y2Sddnp44wgpxDjrQ6jhnwn07t7aVY3j1+uQKnnqoptfA1qFVEcZqETgPXnkipl4zDL8bNwgAXFdHc0OPTpq3fg96TX4tZV0vp0P6N5cMwiVDjvZcmtQ1CMNRPK8hVB2wl4qZeEZvAMAJHoUyx2jlvHt3TJ7gKa4+rZfrc+e2UltTIEIhTSaPtc9asmm7Vje/Lmi+Y1b+NNZlNtqcN4/XIFZ1sBY10Rh6dyyzPaRDu7ezlgU8o597UUHdwa3bdJUp4Lsje+GCwUfhmpG9Up7XzWP7Y+o1/ktjnNG3Av++brj5UNv9KQo3R7MTdWWKImGcN+govPfz4PH06XxPQykpDFuLxg/JQEvo1KYYj377ZCtTeYFjVTtnEhZg/02DLixz8QlH46x+FdbKcoNTVET1w1kkLrl/uuYITPAIWVa4TQoyMRN58etLBlmviyIhXD+mL9bdeaFnQMJ5gxIJkX+87ATXbdRx9XLZiqAVknONCIU0cc5mCsPZG0DU7E/N7oHE7K4+zpbZx1lBUfcp3Dd+KE7RzB+rdx60Bh+di07ognV3XhhIU9CrrCoHdLvSQjx81cmWLdqPH406Fmd5CB8AOMlRKO7E7u1w3RnH4C+ODPKSwrCrw1bHaT7S1fig8zC3GegZfTti1s9GBTxCMMIhwn8mjUwpML04b9BRVgRKTYABhYisgbRNmk7Kbu1LseR/z/PMLXHDuWiRmsx4ZeMqofXSpJE235MXrpF5DVwbW2dEnw6Wf0llzPtpIMqP1619SUZOYGc9LudzAbj7HbONCIU0cd6I2TAfKY7r0ho/HNUHj2tVMpUmEtPqCDkfKj366NRjOtjMBF/sOOgZT+93g+tCQS/VnM2ZGAB8+uvzMG2iPZ4+FCLcetFAV4E17JjkdYvdUPbp4jRDcYHE4DVcWyP5lF7l6NWxzGuXwHxw89l4SauNNLR7u0CC1Qt1PYKWPfiZmTTp5vhOhUqo9OP2SwdjvGmePLt/JwzrlbiG1uJUHnH76lka0r1doP7lSlO446uJc1ChukEsAkooHPZYnCsV919xos2v+H9XJ68E+Or1p+Pknu0xdpB3PaWGIslraeIM1UwlFNIxDUTCISvBTGFpCrG49RA4BwDnOsbOPrVKUYjNDduazdrA6lVeIVPSnbEGRV0Dvb9u9Y7cUJdXhdDu2F9rW+y9IXQvL00rGzgVfiWn3fj+mb1xydCjPTXEhnKVYwGg3dWJkhrKPOIs+qhIN5fAz6fQEL51ak/gVOO1uveD+Jk6miY55+pvQbnw+C648PguVhFCt8lC6+KCpByMbCOaQpo4b1yn5uB8/6/r7LPgdFFOXlW6GfCfFRZEQkmzmgU+6xd74WU3zbZQyBXOa3BKr/aBBwwVV35ZZXe8+7NRhu/EZQGd5kC6TvFQiHImENzQS4EoDSHqkcyVbmSN2+/ZkJXu3Cg276Mgwrd1cQSvXn86HrjypJTb+hF0vehcIZpCmjgHRadTcv6t56I+Hkfl7TNQXlaYVukCN3TzkSUUfKIUCsOhJE1huMsiLanwmrVlmgTUWHRpW4xt+2ps/f/4ljFpaSTd2pfa4uUfdhQtbE4Eqb/UlHxlSBcs3WJk9R+oiWLbvsOea58HjWxy457LhuDBd1fbzFsrfjcWA371ZsbHBIwKxID/ZOi+8UPx3ILNKCuMuDr40+WlSSOzFtWXCc37jmqG6LOTLm2Lk2bUbc1KmPNuHZOVKJaEo9nbfKRTEE7WFPRqkkHxGvwzLReQTV758emei7I8+4MRWLRxr21w8FqXuiWQqwSmbHHdGb3x+9eNVdt++sxibNjtXd/HWUY+Hb5+cjd8/WR7XoVXKGs6FJvPn1+uxbihXTNOXHUjaFJkrhChkAFv3XAmzvvL+74PpMoYbShq1s8MT0ezTjjk5lNI/2f2ctg1B03BL9GrW/vSwMXPWgr/M/pYDO/dAYWREK6ftsjVQdlU6MEMToFw4fFH4fWl2wHANZM5G4w8toMt+TJdlGDJdoBFc0aEQgao2YNXKe1sEg4RLq/shkuHdrWSnLxqqSiWbd1ve59JIo/XPkeKTyGfuPG8RGG0uVPGNGFP3DmrX0XS8qlAduLy7/za8b42eGel2HRRPsJU6zG0JEQoZED38hLceuFxuHiI/6Lt2eLubyRq0MyZPNo1XO8f1w7Dh6urAACLNrqvaJYNsu3IE1o+Pzmnr6tQcMufSRe/Ba6ygaoO3DFAmfSWgjzhGUBEuO7M3ujStvGiOBRe8dtn9K3AFLO2zxNankO2aQ7mI+HIwuuOueHcfh6fNB9KTEf+cV38tfOWhGgKLZDRWg2WbCPmIyFdDtQkr1sMGOuXjxt6NHp2aHhSYK64clgPdCwrxNjBuUsWa26IUGiBZKMYmBdBVvwSBJ2Te7pXnSUiz9LTzYVwiHDB8Y1jJm4uiFBooaz43VjMXl3lOUtLl+d/OAKvL92eU4EjtEzSzboWmhb5tVooxQVhjAlY1sGLJ797ilVN8+Se5Ti5Z7C6Q4IgHLmIgVjw5Oz+nTCoAaWSBUHxp8uSV3ETmiciFARByDnfOLkbfplBZr3Q+IhQEAShUThTW1OjtfgZmi0iFARBaBT0sjBzpqReREdoFAGt6QAAB4NJREFUGkQoCILQKKj6QaWF4YxWJhMaBxEKgiA0Cs29oqtgIEJBEIRGwW2lNKH5Ib+SIAiNgmgKRwYiFARBaBSUT8Fn4UChGZAzoUBE3YnoXSJaTkTLiOgnZns5Eb1NRKvM/+21faYQ0WoiWklE5+eqb4IgND5hqZt1RJBLTSEK4CZmPg7AcACTiGgggMkAZjJzXwAzzfcwPxsPYBCAsQAeIqKGr6cnCEKzICx1s44IciYUmHkbMy80Xx8AsBxAVwDjAEw1N5sK4FLz9TgATzNzLTOvA7AawLBc9U8QhMZFrXNc6rG+ttA8aJS0QiLqBeBEAB8D6MzM2wBDcBCRKv7fFcBH2m6bzTbnsSYCmAgAPXrkdtUlQRCyR0lhGJMvGIBzBzasUKOQW3LuaCaiVgCeB/BTZt7vt6lLW5JLipkfZeZKZq6sqKhw2UUQhObKD87qgz4VrZq6G4IPORUKRFQAQyD8i5lfMJt3EFEX8/MuAHaa7ZsBdNd27wZgay77JwiCINjJZfQRAXgCwHJm/rP20csAJpivJwB4SWsfT0RFRHQMgL4A5uWqf4IgCEIyufQpjATwbQBLiWix2XYLgLsATCeiawFsBHAZADDzMiKaDuBzGJFLk5g5lsP+CYIgCA5yJhSY+UO4+wkAYIzHPncAuCNXfRIEQRD8kYxmQRAEwUKEgiAIgmAhQkEQBEGwEKEgCIIgWBAfwSULiWgXgA0NOERHAFVZ6s6RQL6dLyDnnC/IOadHT2Z2zf49ooVCQyGi+cxc2dT9aCzy7XwBOed8Qc45e4j5SBAEQbAQoSAIgiBY5LtQeLSpO9DI5Nv5AnLO+YKcc5bIa5+CIAiCYCffNQVBEARBIy+FAhGNNdeBXk1Ek5u6P9kiX9fFJqIwES0iolfN9y39fNsR0XNEtML8rUfkwTnfYN7TnxHRNCIqbmnnTET/R0Q7iegzrS3tcySik4loqfnZ/WbF6uAwc179AQgDWAOgN4BCAEsADGzqfmXp3LoAOMl83RrAFwAGArgbwGSzfTKAP5ivB5rnXwTgGPO6hJv6PDI47xsB/BvAq+b7ln6+UwF8z3xdCKBdSz5nGCswrgNQYr6fDuDqlnbOAM4EcBKAz7S2tM8RxpIDI2AUJH0DwAXp9CMfNYVhAFYz81pmrgPwNIz1oY94OA/XxSaibgAuAvC41tySz7cNjMHjCQBg5jpm3osWfM4mEQAlRBQBUApjAa4Wdc7M/D6APY7mtM7RXLisDTPPZUNC/F3bJxD5KBS6AtikvXddC/pIx29dbAD6uthH+rW4F8DNAOJaW0s+394AdgF40jSZPU5EZWjB58zMWwD8Ccb6K9sA7GPmt9CCz1kj3XPsar52tgcmH4VCoLWgj2SyvS52c4WILgawk5kXBN3Fpe2IOV+TCAwTw8PMfCKAahhmBS+O+HM27ejjYJhJjgZQRkRX+e3i0nZEnXMAvM6xweeej0KhRa8FnWfrYo8EcAkRrYdhBhxNRP9Eyz1fwDiHzcz8sfn+ORhCoiWf8zkA1jHzLmauB/ACgNPQss9Zke45bjZfO9sDk49C4RMAfYnoGCIqBDAexvrQRzz5ti42M09h5m7M3AvG7/gOM1+FFnq+AMDM2wFsIqL+ZtMYGEvYtthzhmE2Gk5EpeY9PgaGv6wln7MirXM0TUwHiGi4ea2+o+0TjKb2uDeRl/9CGJE5awDc2tT9yeJ5nQ5DVfwUwGLz70IAHQDMBLDK/F+u7XOreR1WIs0oheb0B2AUEtFHLfp8AQwFMN/8nf8DoH0enPNvAKwA8BmAf8CIumlR5wxgGgyfST2MGf+1mZwjgErzOq0B8CDMJOWgf5LRLAiCIFjko/lIEARB8ECEgiAIgmAhQkEQBEGwEKEgCIIgWIhQEARBECxEKAhCAIjoVrNK56dEtJiITiWinxJRaVP3TRCyiYSkCkIKiGgEgD8DGMXMtUTUEUZ10jkAKpm5qkk7KAhZRDQFQUhNFwBVzFwLAKYQ+AaMOjzvEtG7AEBE5xHRXCJaSETPmjWoQETriegPRDTP/DvWbL/MXB9gCRG93zSnJgh2RFMQhBSYg/uHMEo2zwDwDDO/Z9ZcqmTmKlN7eAFGZmk1Ef0CQBEz/9bc7jFmvoOIvgPgcma+mIiWAhjLzFuIqB0bJbAFoUkRTUEQUsDMBwGcDGAijLLVzxDR1Y7NhsNY+GQ2ES2GUaemp/b5NO3/CPP1bABPEdF1MBZ/EoQmJ9LUHRCEIwFmjgGYBWCWOcOf4NiEALzNzFd4HcL5mpl/QESnwlgkaDERDWXm3dntuSCkh2gKgpACIupPRH21pqEANgA4AGPZUwD4CMBIzV9QSkT9tH2+qf2fa27Th5k/ZubbAFTBXgpZEJoE0RQEITWtADxARO0ARGEsfTgRwBUA3iCibcx8tmlSmkZEReZ+v4RRjRcAiojoYxgTMaVN/NEUNgSjAuaSRjkbQfBBHM2CkGN0h3RT90UQUiHmI0EQBMFCNAVBEATBQjQFQRAEwUKEgiAIgmAhQkEQBEGwEKEgCIIgWIhQEARBECxEKAiCIAgW/w9pD9B1aJ1cZwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "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 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": 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": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEMCAYAAADDMN02AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAP60lEQVR4nO3de5BkZX3G8e8jK+EiirojRcDJYEKIiuJlSIxoVIgVYFWi5TVeiKEySZkYvMSIRRJzrVoqJoV3aoNkNaHQiBgvGBUvgCkB3SWEi+uF6IIbiYuSUgNWkcVf/uiGHYbZnd7Z7XN65/1+qqamz9un+/3tW7vPnj7n9PumqpAkteN+fRcgSeqWwS9JjTH4JakxBr8kNcbgl6TGGPyS1JhVfRcwitWrV9fMzEzfZUjSXmXjxo3fq6qphe17RfDPzMywYcOGvsuQpL1KkpsWa/dUjyQ1xuCXpMYY/JLUGINfkhpj8EtSY8YW/EnOS7I1yfXz2v4myVeTXJvkw0kOHlf/kqTFjfOIfz1w4oK2S4Cjq+qxwNeBN42xf0nSIsYW/FV1OXDbgrZPV9W24eaVwOHj6l+StLg+v8D1W8AHdvRkkjlgDmB6erqrmqSRzZxx8T2PN69d02Ml0q7p5eJukjOBbcD5O9qnqtZV1WxVzU5N3ecbx5KkZer8iD/JqcCzgBPKdR8lqXOdBn+SE4E3Ak+rqju67FuSNDDO2zkvAK4AjkqyJclpwDuAg4BLklyT5Jxx9S9JWtzYjvir6iWLNL9nXP1JkkbjN3clqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTFjC/4k5yXZmuT6eW0PSXJJkm8Mfz94XP1LkhY3ziP+9cCJC9rOAD5bVUcCnx1uS5I6NLbgr6rLgdsWNJ8CvHf4+L3Ar4+rf0nS4lZ13N8hVXULQFXdkuRhO9oxyRwwBzA9Pd1RedK9zZxx8T2PN69d02Ml0p4zsRd3q2pdVc1W1ezU1FTf5UjSitF18H83yaEAw99bO+5fkprXdfB/FDh1+PhU4CMd9y9JzRvn7ZwXAFcARyXZkuQ0YC3wzCTfAJ453JYkdWhsF3er6iU7eOqEcfUpSVraxF7clSSNh8EvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY3pJfiTvDbJDUmuT3JBkv36qEOSWtR58Cc5DPgDYLaqjgb2AV7cdR2S1Kq+TvWsAvZPsgo4APhOT3VIUnM6D/6q+i/gLcDNwC3AD6rq0wv3SzKXZEOSDbfeemvXZUrSitXHqZ4HA6cARwA/DRyY5GUL96uqdVU1W1WzU1NTXZcpSStWH6d6fhX4VlXdWlX/B1wEPLmHOiSpSX0E/83Ak5IckCTACcCmHuqQpCb1cY7/KuBC4GrgumEN67quQ5JataqPTqvqzcCb++hbklrnN3clqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGjBT8SY4bpU2SNPlGPeJ/+4htkqQJt9PZOZP8MoNFUqaSvG7eUw9ksEi6JGkvs9S0zPsCDxjud9C89h8Czx9XUZKk8dlp8FfVZcBlSdZX1U0d1SRJGqNRF2L5qSTrgJn5r6mq48dRlCRpfEYN/g8C5wDnAneNrxxJ0riNGvzbqurdY61EktSJUW/n/FiSVyU5NMlD7v4Za2WSpLEY9Yj/1OHvN8xrK+ARe7YcSdK4jRT8VXXEuAuRJHVjpOBP8orF2qvqfXu2HEnSuI16qufYeY/3A04ArgYMfknay4x6qufV87eTPAj4x7FUJEkaq+VOy3wHcORyO01ycJILk3w1yabhnECSpA6Meo7/Ywzu4oHB5GyPBP55N/p9K/DJqnp+kn2BA3bjvSRJu2DUc/xvmfd4G3BTVW1ZTodJHgj8CvCbAFV1J3Dnct5LkrTrRj3Hf1mSQ9h+kfcbu9HnI4BbgX9IcgywETi9qm6fv1OSOWAOYHp6eje6k+5r5oyL73m8ee2aHiuRujfqClwvBL4EvAB4IXBVkuVOy7wKeALw7qp6PHA7cMbCnapqXVXNVtXs1NTUMruSJC006qmeM4Fjq2orQJIp4DPAhcvocwuwpaquGm5fyCLBL0kaj1Hv6rnf3aE/9P1deO29VNV/A99OctSw6QTgK8t5L0nSrhv1iP+TST4FXDDcfhHwid3o99XA+cM7er4JvHI33kuStAuWWnP354BDquoNSZ4HPAUIcAVw/nI7raprgNnlvl6StHxLna45G/gRQFVdVFWvq6rXMjjaP3vcxUmS9rylgn+mqq5d2FhVGxgswyhJ2sssFfz77eS5/fdkIZKkbiwV/F9O8tsLG5OcxuCLV5KkvcxSd/W8BvhwkpeyPehngX2B546zMEnSeOw0+Kvqu8CTkzwDOHrYfHFVfW7slUmSxmLUuXo+D3x+zLVIkjqw3Pn4JUl7KYNfkhpj8EtSYwx+SWqMwS9JjRl1dk6pGfNX5xqlXdrbeMQvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUmN6CP8k+Sf49ycf7qkGSWtTnEf/pwKYe+5ekJvUS/EkOB9YA5/bRvyS1rK+FWM4G/gg4aEc7JJkD5gCmp6c7KktanvmLtGxeu6bHSqSldX7En+RZwNaq2riz/apqXVXNVtXs1NRUR9VJ0srXx6me44DnJNkMvB84Psk/9VCHJDWp8+CvqjdV1eFVNQO8GPhcVb2s6zokqVXexy9Jjenr4i4AVXUpcGmfNUhSazzil6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGtPrXD3SJJi/iMo43s+FWTRpPOKXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqTOfBn+ThST6fZFOSG5Kc3nUNktSyPqZl3ga8vqquTnIQsDHJJVX1lR5qkaTmdH7EX1W3VNXVw8c/AjYBh3VdhyS1qteFWJLMAI8HrlrkuTlgDmB6errTurTy7OnFVqS9WW8Xd5M8APgQ8Jqq+uHC56tqXVXNVtXs1NRU9wVK0grVS/AnuT+D0D+/qi7qowZJalUfd/UEeA+wqar+ruv+Jal1fRzxHwe8HDg+yTXDn5N7qEOSmtT5xd2q+jcgXfcrSRrwm7uS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JakyvK3BJu2tvWFlrfo2b167psRJpwCN+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4JekxvQS/ElOTPK1JDcmOaOPGiSpVZ0Hf5J9gHcCJwGPAl6S5FFd1yFJrerjiP8XgRur6ptVdSfwfuCUHuqQpCb1sRDLYcC3521vAX5p4U5J5oC54eb/JvlaB7XtzGrgez3XMCkci+12aSxy1hgr6Z9/L7ablLH4mcUa+wj+LNJW92moWgesG385o0myoapm+65jEjgW2zkW2zkW2036WPRxqmcL8PB524cD3+mhDklqUh/B/2XgyCRHJNkXeDHw0R7qkKQmdX6qp6q2Jfl94FPAPsB5VXVD13Usw8ScdpoAjsV2jsV2jsV2Ez0WqbrP6XVJ0grmN3clqTEGvyQ1xuBfYJTpJJI8Pck1SW5IclnXNXZlqbFI8qAkH0vyH8OxeGUfdY5bkvOSbE1y/Q6eT5K3Dcfp2iRP6LrGrowwFi8djsG1Sb6Y5Jiua+zKUmMxb79jk9yV5Pld1bYUg3+eUaaTSHIw8C7gOVX1aOAFnRfagRGn1vg94CtVdQzwdOBvh3dqrTTrgRN38vxJwJHDnzng3R3U1Jf17HwsvgU8raoeC/wlE36RczetZ+djcfe/o7MY3MwyMQz+extlOonfAC6qqpsBqmprxzV2ZZSxKOCgJAEeANwGbOu2zPGrqssZ/Nl25BTgfTVwJXBwkkO7qa5bS41FVX2xqv5nuHklg+/prEgj/L0AeDXwIWCicsLgv7fFppM4bME+Pw88OMmlSTYmeUVn1XVrlLF4B/BIBl/Auw44vap+0k15E2WUsWrRacC/9l1EX5IcBjwXOKfvWhbqY8qGSTbKdBKrgCcCJwD7A1ckubKqvj7u4jo2ylj8GnANcDzws8AlSb5QVT8cd3ETZqRpSFqS5BkMgv8pfdfSo7OBN1bVXYMPxZPD4L+3UaaT2AJ8r6puB25PcjlwDLDSgn+UsXglsLYGXwa5Mcm3gF8AvtRNiRPDaUjmSfJY4FzgpKr6ft/19GgWeP8w9FcDJyfZVlX/0m9ZnupZaJTpJD4CPDXJqiQHMJhZdFPHdXZhlLG4mcEnH5IcAhwFfLPTKifDR4FXDO/ueRLwg6q6pe+i+pBkGrgIePkK/BS8S6rqiKqaqaoZ4ELgVZMQ+uAR/73saDqJJL87fP6cqtqU5JPAtcBPgHOraqe3c+2NRhkLBndtrE9yHYPTHW+sqkmYinaPSnIBg7uWVifZArwZuD/cMw6fAE4GbgTuYPBJaEUaYSz+FHgo8K7hke62SZ6lcneMMBYTyykbJKkxnuqRpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4pQ4keXuSq5Mc23ctksEvjVmSA4GHAb8DPKvnciSDX7pbkj9L8oe78fqZJD9Ocs389uG8TocClwJvG+67/3AxnzuTrN6duqVdZfBLe9Z/VtXj5jckeShwAPAj4C6AqvrxcL9mJ3NTfwx+NS3JmcPlJT/DYJK5cfhj4C3ADQxWM5N6ZfCrWUmeyGDW0ccDzwP2+IXXJDPAk4EPMJjF9dF7ug9pVxn8atlTgQ9X1R3DxWPumXY6yQuSvDXJO5L89bDtI/Oe/+BwPdWl/BXwF8M1Cwx+TQSnZVbr7jM9bZLjgNmqOn24fU6SpwHz59i/X1XdtbM3TvI4Bp8knpLkncB+DJaolHrlEb9adjnw3OEdNgcBzx62nwa8fcG+TwAeNfxP4L2MdlH2LODZ8xbjOAaP+DUBPOJXs6rq6iQfYLBu8E3AF4ZP3Z/hJ4EkRzC4B/9W4PVV9eUka4Cpnb13kuOBA6vqs/P6+26SA5M8pKpu2/N/Imk0LsQiLZDkMcCZwFYG/wn8CfD3wIuq6s4kfw5cWFXXLXjdDPDxqjp6F/razOC00opbuUyTy+CX9pAkDwe+CHx/4b38i+y7P3AFg08Oj/ETgLpk8EtSY7y4K0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWrM/wOQAeGw/vRN5QAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "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 controlling the cutoff distance" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "neighbors = final_struct.get_neighbors(cutoff_radius=8)" ] }, { "cell_type": "code", "execution_count": 16, "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": 17, "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": 18, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEcCAYAAAAC+llsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAb90lEQVR4nO3dfZQldX3n8fdHhKCAIk6Dw5NDEsIGTEB3gg+YBCR6mAElusYwcY1x3YwazWo0biayMSZnY3DXRFc0skSJmihCVjEoo0KMkRhFHVjkIUBAMoRxRmYA5UFYzeB3/6gavTT3dt/p6vvQ4/t1Tp+u+6tfVX27pqY//au6typVhSRJXTxs0gVIkpY+w0SS1JlhIknqzDCRJHVmmEiSOjNMJEmdGSbSEJLcm+RHJ13HfJKsSFJJHr7A5d+Q5D2LXZd2fYaJJirJryW5Osl9Sb6R5N1J9h1iuYOTfDDJHUm+neTLSU4ZVZ1VtXdV3Tyq9U9CkuOTbOptq6o3V9V/nlRNWroME01MktcBbwFeDzwaeArweOCSJHvMsdx+wOeB7wJHAcuAtwEfSvL8Udfdp57dRrTeh4wuFjrikEbNMNFEJHkU8AfAb1bVp6rq36pqI/ACmkD5j3Ms/lvAvcBLq+obVXV/VZ0L/BHwJ0nSZ3s7Tv+sTbI5yZY2zHbMPzbJF5N8q533zt5Aa5f98Xb6fe0Ian2SbwMn9Nnefkn+ot3WN5N8rGferye5KcmdSS5McuCs7bwyyY3AjTtGD0l+J8k3gL9I8rAk65J8rR2Znd8GbL/9/JIk1yW5J8nNSV7Wtu8FfBI4sD2Fd2+SA5O8Kclf9Sz/nCTXtvvl75P8ZM+8jUl+O8lVSe5Kcl6SPef4d9MuzDDRpDwN2BP4aG9jVd1L80vumXMs+0zgI1X1vVnt5wOHAj8xx7InAIcDzwLWJfmFtv0BmpBaBjwVOBH4jTnW8ys04bUPzShptr8EHkkzctqfZuREkmcAf0wTmsuBW4APz1r2F4EnA0e2rx8H7EcTsmuB/9L2+XngQOCbwLsG1LkVOAV4FPAS4G1JnlRV3wZWAZvbU3h7V9Xm3gWT/ARwLvAaYAZYD3x81qjxBcBJwGHATwO/NqAO7eIME03KMuD2qtreZ96Wdv5cy24ZsNyO+YP8QVV9u6quBv4CWANQVZdX1WVVtb0dIf1vml/Wg/xNVf1jVX2vqv5f74wky2l+Ub+8qr7Zjro+185+IXBOVV1RVd8Bfhd4apIVPav446q6s6rub19/D/j9qvpO2/Yy4PSq2tSu403A8/udAquqi6rqa9X4HHAx8LNz/Fy9fhm4qKouqap/A94KPILmD4Ed3lFVm6vqTuDjwDFDrlu7GMNEk3I7sGzANYDl7XySnNVzGuYNPcsuH7DcjvmD3NozfQvNX/Yk+Ykkn2jfBHA38GbmDqVb55h3CHBnVX2zz7wD2+0C3x+J3QEcNMe6t80KrMcDF7Snnr4FXEczsjpg9saSrEpyWXtK7VvA6nl+rrlq/V5bW2+t3+iZvg/Ye8h1axdjmGhSvgh8B3heb2N7Ln8V8BmAqnp5z2mYN7fd/hb4D0lmH78voPll989zbPeQnulDgR2ndt4NXA8cXlWPAt4APOTaS4+5brd9K7DfgHelbaYJA+D7P+9jga/Pse7Zr28FVlXVvj1fe1ZV7zpI8iPAR2hGFAdU1b40p6p2/Fzz3TJ8dq2h2X9fH7iEfmgZJpqIqrqL5gL8mUlOSrJ7e6rnr4FNNNccBnkbzTWA9yZ5XJI9k6wBTgdeX3M/V+H3kjwyyVE01xDOa9v3Ae4G7k3y74BXdPjZttBc9/mzJI9pf7afa2d/CHhJkmPaX/ZvBr7Unlob1lnAHyV5PECSmSSn9um3B/AjwDZge5JVNNeKdrgNeGySRw/YzvnAyUlOTLI78DqaPwC+sBO16oeEYaKJqar/QTMCeCvNL/Iv0fzVfWJ7LWDQcncAT6e5gP9PNKeJXgu8qKrOG7Rc63PATTQjn7dW1cVt+2/TXFS/B/hzfhAyC/Ui4N9oRjtbaS5iU1WfAX6PZsSwBfgx4LSdXPf/Ai4ELk5yD3AZzQX7B6mqe2gu1p9Pc5H+V9rldsy/nuYC+83tKbMDZy1/A8276s6kOXX4bODZVfXdnaxXPwTiw7H0w6Ad9fwLsPuAi/6SOnBkIknqzDCRJHXmaS5JUmeOTCRJnRkmkqTOdqk7kC5btqxWrFgx6TIkacm4/PLLb6+qma7r2aXCZMWKFWzYsGHSZUjSkpHklvl7zc/TXJKkzgwTSVJnhokkqTPDRJLUmWEiSepsZGGS5JAkn22fP31tkle37fsluSTJje33xwxY/qQkN7TPyl43qjolSd2NcmSyHXhdVf0k8BTglUmOBNYBn6mqw2luA/6QoEiyG80zrVfRPAd7TbusJGkKjSxMqmpLVV3RTt9D82jRg4BTgfe33d4P/GKfxY8Fbqqqm9tnJ3y4XU6SNIXGcs2kfZbEE2kefnRA+yS6HU+k27/PIgfx4Odgb+LBz53uXffaJBuSbNi2bdtili1JGtLIwyTJ3jRPlXtNVd097GJ92vre3riqzq6qlVW1cmam8x0BJEkLMNIwaZ8b/RHgg1X10bb5tiTL2/nLaR5pOtsm4JCe1wcDm0dZqyRp4Ub5bq4A7wWuq6o/7Zl1IfDidvrFwN/0WfwrwOFJDkuyB80zsi/s00+SNAVGOTI5DngR8IwkV7Zfq4EzgGcmuRF4ZvuaJAcmWQ/QPqP7VcCnaS7cn19V146wVklSByO7a3BVfZ7+1z4ATuzTfzOwuuf1emD9aKqTJC0mPwEvSerMMJEkdWaYSJI6M0wkSZ0ZJpKkzgwTSVJnhokkqTPDRJLUmWEiSerMMJEkdWaYSJI6M0wkSZ0ZJpKkzgwTSVJnhokkqbORPc9EkqbNinUXDZy38YyTx1jJrmdkYZLkHOAUYGtVPaFtOw84ou2yL/Ctqjqmz7IbgXuAB4DtVbVyVHVKkrob5cjkfcA7gQ/saKiqX94xneRPgLvmWP6Eqrp9ZNVJkhbNKB/be2mSFf3mJQnwAuAZo9q+JGl8JnUB/meB26rqxgHzC7g4yeVJ1o6xLknSAkzqAvwa4Nw55h9XVZuT7A9ckuT6qrq0X8c2bNYCHHrooYtfqSRpXmMfmSR5OPA84LxBfapqc/t9K3ABcOwcfc+uqpVVtXJmZmaxy5UkDWESp7l+Abi+qjb1m5lkryT77JgGngVcM8b6JEk7aWRhkuRc4IvAEUk2JXlpO+s0Zp3iSnJgkvXtywOAzyf5KvBl4KKq+tSo6pQkdTfKd3OtGdD+a33aNgOr2+mbgaNHVZckafF5OxVJUmeGiSSpM8NEktSZYSJJ6swwkSR1ZphIkjozTCRJnRkmkqTODBNJUmeGiSSpM8NEktSZYSJJ6swwkSR1ZphIkjozTCRJnRkmkqTODBNJUmejfGzvOUm2Jrmmp+1NSb6e5Mr2a/WAZU9KckOSm5KsG1WNkqTFMcqRyfuAk/q0v62qjmm/1s+emWQ34F3AKuBIYE2SI0dYpySpo5GFSVVdCty5gEWPBW6qqpur6rvAh4FTF7U4SdKimsQ1k1cluao9DfaYPvMPAm7teb2pbZMkTalxh8m7gR8DjgG2AH/Sp0/6tNWgFSZZm2RDkg3btm1bnColSTtlrGFSVbdV1QNV9T3gz2lOac22CTik5/XBwOY51nl2Va2sqpUzMzOLW7AkaShjDZMky3tePhe4pk+3rwCHJzksyR7AacCF46hPkrQwDx/VipOcCxwPLEuyCfh94Pgkx9CcttoIvKzteyDwnqpaXVXbk7wK+DSwG3BOVV07qjolSd2NLEyqak2f5vcO6LsZWN3zej3wkLcNS5Kmk5+AlyR1ZphIkjozTCRJnRkmkqTODBNJUmcjezeXNIwV6y4aOG/jGSePsRJJXTgykSR1ZphIkjozTCRJnRkmkqTODBNJUmeGiSSpM8NEktSZYSJJ6swwkSR1ZphIkjozTCRJnY0sTJKck2Rrkmt62v5nkuuTXJXkgiT7Dlh2Y5Krk1yZZMOoapQkLY5RjkzeB5w0q+0S4AlV9dPAPwO/O8fyJ1TVMVW1ckT1SZIWycjCpKouBe6c1XZxVW1vX14GHDyq7UuSxmeS10z+E/DJAfMKuDjJ5UnWjrEmSdICTOR5JklOB7YDHxzQ5biq2pxkf+CSJNe3I51+61oLrAU49NBDR1KvJGluYx+ZJHkxcArwwqqqfn2qanP7fStwAXDsoPVV1dlVtbKqVs7MzIyiZEnSPMYaJklOAn4HeE5V3Tegz15J9tkxDTwLuKZfX0nSdBjlW4PPBb4IHJFkU5KXAu8E9qE5dXVlkrPavgcmWd8uegDw+SRfBb4MXFRVnxpVnZKk7kZ2zaSq1vRpfu+AvpuB1e30zcDRo6pLkrT4JnIBXhqlFesuGjhv4xknj7ESqb9d8Rj1diqSpM4ME0lSZ0OFSZLjhmmTJP1wGnZkcuaQbZKkH0JzXoBP8lTgacBMktf2zHoUsNsoC5MkLR3zvZtrD2Dvtt8+Pe13A88fVVGSpKVlzjCpqs8Bn0vyvqq6ZUw1SZKWmGE/Z/IjSc4GVvQuU1XPGEVRkqSlZdgw+WvgLOA9wAOjK0eStBQNGybbq+rdI61EkrRkDfvW4I8n+Y0ky5Pst+NrpJVJkpaMYUcmL26/v76nrYAfXdxyJElL0VBhUlWHjboQSdLSNVSYJPnVfu1V9YHFLUeStBQNe5rrZ3qm9wROBK4ADBNJ0tCnuX6z93WSRwN/OZKKJElLzkJvQX8fcPhcHZKck2Rrkmt62vZLckmSG9vvjxmw7ElJbkhyU5J1C6xRkjQmw96C/uNJLmy/LgJuAP5mnsXeB5w0q20d8JmqOhz4TPt69rZ2A94FrAKOBNYkOXKYOiVJkzHsNZO39kxvB26pqk1zLVBVlyZZMav5VOD4dvr9wN8DvzOrz7HATe2z4Eny4Xa5fxqyVknSmA17zeRzSQ7gBxfib1zg9g6oqi3tOrck2b9Pn4OAW3tebwKePGiFSdYCawEOPfTQBZalabQrPicbdt2fa6HcH7uGYU9zvQD4MvBLwAuALyUZ1S3o06etBnWuqrOramVVrZyZmRlRSZKkuQx7mut04GeqaitAkhngb4H/s5Pbuy3J8nZUshzY2qfPJuCQntcHA5t3cjuSpDEa9t1cD9sRJK07dmLZXhfyg1uzvJj+F/G/Ahye5LAkewCntctJkqbUsCOTTyX5NHBu+/qXgfVzLZDkXJqL7cuSbAJ+HzgDOD/JS4F/pTltRpIDgfdU1eqq2p7kVcCnaR4NfE5VXbtzP5YkaZzmewb8j9NcNH99kucBT6e5pvFF4INzLVtVawbMOrFP383A6p7X65knrCRJ02O+U1VvB+4BqKqPVtVrq+q3aH7Rv33UxUmSlob5wmRFVV01u7GqNtA8wleSpHnDZM855j1iMQuRJC1d84XJV5L8+uzG9gL65aMpSZK01Mz3bq7XABckeSE/CI+VwB7Ac0dZmCRp6ZgzTKrqNuBpSU4AntA2X1RVfzfyyiRJS8aw9+b6LPDZEdciSVqiFvo8E0mSvs8wkSR1ZphIkjozTCRJnRkmkqTODBNJUmeGiSSpM8NEktTZsA/HkjpZse6iSZcwEnP9XBvPOHnJbmshxl3foO1Nw75YKhbz/6UjE0lSZ2MPkyRHJLmy5+vuJK+Z1ef4JHf19HnjuOuUJA1v7Ke5quoG4BiAJLsBXwcu6NP1H6rqlHHWJklamEmf5joR+FpV3TLhOiRJHUw6TE4Dzh0w76lJvprkk0mOGrSCJGuTbEiyYdu2baOpUpI0p4mFSZI9gOcAf91n9hXA46vqaOBM4GOD1lNVZ1fVyqpaOTMzM5piJUlzmuTIZBVwRfsArgepqrur6t52ej2we5Jl4y5QkjScSYbJGgac4kryuCRpp4+lqfOOMdYmSdoJE/nQYpJHAs8EXtbT9nKAqjoLeD7wiiTbgfuB06qqJlGrJGl+EwmTqroPeOystrN6pt8JvHPcdUmSFmbS7+aSJO0CDBNJUmeGiSSpM8NEktSZYSJJ6swwkSR1ZphIkjozTCRJnRkmkqTODBNJUmcTuZ2K1NWKdRct6nIbzzi5Szk7ta2lbtp/rnHWtxS2NYpjux9HJpKkzgwTSVJnhokkqTPDRJLUmWEiSerMMJEkdTaRMEmyMcnVSa5MsqHP/CR5R5KbklyV5EmTqFOSNJxJfs7khKq6fcC8VcDh7deTgXe33yVJU2haT3OdCnygGpcB+yZZPumiJEn9TSpMCrg4yeVJ1vaZfxBwa8/rTW3bQyRZm2RDkg3btm0bQamSpPlMKkyOq6on0ZzOemWSn5s1P32WqX4rqqqzq2plVa2cmZlZ7DolSUOYSJhU1eb2+1bgAuDYWV02AYf0vD4Y2Dye6iRJO2vsYZJkryT77JgGngVcM6vbhcCvtu/qegpwV1VtGXOpkqQhTeLdXAcAFyTZsf0PVdWnkrwcoKrOAtYDq4GbgPuAl0ygTknSkMYeJlV1M3B0n/azeqYLeOU465IkLdy0vjVYkrSEGCaSpM4ME0lSZ4aJJKkzw0SS1Nkkb/SoXcyKdRdNuoQFm6v2jWecPBV1LPZy4/y5FmopH1MLtVR/ZkcmkqTODBNJUmeGiSSpM8NEktSZYSJJ6swwkSR1ZphIkjozTCRJnRkmkqTODBNJUmeGiSSps0k8A/6QJJ9Ncl2Sa5O8uk+f45PcleTK9uuN465TkjS8SdzocTvwuqq6Isk+wOVJLqmqf5rV7x+q6pQJ1CdJ2kljH5lU1ZaquqKdvge4Djho3HVIkhbPRK+ZJFkBPBH4Up/ZT03y1SSfTHLUHOtYm2RDkg3btm0bUaWSpLlMLEyS7A18BHhNVd09a/YVwOOr6mjgTOBjg9ZTVWdX1cqqWjkzMzO6giVJA00kTJLsThMkH6yqj86eX1V3V9W97fR6YPcky8ZcpiRpSJN4N1eA9wLXVdWfDujzuLYfSY6lqfOO8VUpSdoZk3g313HAi4Crk1zZtr0BOBSgqs4Cng+8Isl24H7gtKqqCdQqSRpCdqXf0StXrqwNGzZMuowlY9Czpud6NvhSfT61do7HwIMN2h+7wr645S2nXF5VK7uux0/AS5I6M0wkSZ0ZJpKkzgwTSVJnhokkqTPDRJLUmWEiSerMMJEkdWaYSJI6M0wkSZ0ZJpKkzgwTSVJnhokkqTPDRJLUmWEiSerMMJEkdWaYSJI6m0iYJDkpyQ1Jbkqyrs/8JHlHO/+qJE+aRJ2SpOGMPUyS7Aa8C1gFHAmsSXLkrG6rgMPbr7XAu8dapCRpp0xiZHIscFNV3VxV3wU+DJw6q8+pwAeqcRmwb5Ll4y5UkjSch09gmwcBt/a83gQ8eYg+BwFbZq8syVqa0QvAd5Jcs3iljsQy4PZJFzGXvAVYAnW2rHNxLQNub4+BaTbW/dlhfyyFf/cjFmMlkwiT9GmrBfRpGqvOBs4GSLKhqlZ2K2+0lkKNYJ2LzToXl3UuniQbFmM9kzjNtQk4pOf1wcDmBfSRJE2JSYTJV4DDkxyWZA/gNODCWX0uBH61fVfXU4C7quohp7gkSdNh7Ke5qmp7klcBnwZ2A86pqmuTvLydfxawHlgN3ATcB7xkyNWfPYKSF9tSqBGsc7FZ5+KyzsWzKDWmqu+lCEmShuYn4CVJnRkmkqTOpj5MkhyS5LNJrktybZJX9+kz8PYr8926Zcx1vrCt76okX0hydM+8jUmuTnLlYr1Vr0Odxye5q63lyiRv7Jk3Tfvz9T01XpPkgST7tfPGtT/3TPLlJF9t6/yDPn0menwOWeM0HJvD1DkNx+YwdU782OypZbck/zfJJ/rMW7xjs6qm+gtYDjypnd4H+GfgyFl9VgOfpPl8ylOAL7XtuwFfA34U2AP46uxlx1zn04DHtNOrdtTZvt4ILJuS/Xk88Ik+y07V/pzV/9nA301gfwbYu53eHfgS8JRpOj6HrHEajs1h6pyGY3PeOqfh2OzZ3muBDw3Yb4t2bE79yKSqtlTVFe30PcB1NJ+G7zXo9ivD3LplbHVW1Req6pvty8toPj8zVkPuz0Gman/OsgY4dxS1zKU95u5tX+7efs1+V8tEj89hapySY3OYfTnIOI/Nna1zIscmQJKDgZOB9wzosmjH5tSHSa8kK4An0vwl0GvQ7VcGtY/UHHX2einNXwQ7FHBxksvT3CJm5Oap86ntMP6TSY5q26ZyfyZ5JHAS8JGe5rHtz/Y0wpXAVuCSqpq643OIGntN7Ngcss6JH5vD7s9JH5vA24H/CnxvwPxFOzYncTuVBUmyN80/yGuq6u7Zs/ssUnO0j8w8de7ocwLNf9in9zQfV1Wbk+wPXJLk+qq6dEJ1XgE8vqruTbIa+BjNHZyncn/SnEb4x6q6s6dtbPuzqh4AjkmyL3BBkidUVe894iZ+fA5RIzD5Y3OIOqfi2Bx2fzLBYzPJKcDWqro8yfGDuvVpW9CxuSRGJkl2p/mF8sGq+mifLoNuvzLW27IMUSdJfppmyHlqVd2xo72qNrfftwIX0AwzJ1JnVd29YxhfVeuB3ZMsYwr3Z+s0Zp1GGOf+7Nnmt4C/p/lLtNdUHJ8wZ41TcWzOV+e0HJvz1dljksfmccBzkmykOU31jCR/NavP4h2bc11QmYYvmoT8APD2OfqczIMvIn25bX84cDNwGD+4iHTUBOs8lOZT/U+b1b4XsE/P9BeAkyZY5+P4wQdajwX+tV1uqvZn2+/RwJ3AXhPanzPAvu30I4B/AE6ZpuNzyBqn4dgcps5pODbnrXMajs1ZtRxP/wvwi3ZsLoXTXMcBLwKubs9RAryB5uCn5rj9Sg24dcsE63wj8Fjgz5IAbK/mjqIH0AyVoflH/FBVfWqCdT4feEWS7cD9wGnVHGHTtj8BngtcXFXf7ll2nPtzOfD+NA99exhwflV9IkPcHmiMx+cwNU7DsTlMndNwbA5TJ0z+2OxrVMemt1ORJHW2JK6ZSJKmm2EiSerMMJEkdWaYSJI6M0wkSZ0ZJpKkzgwTaYokOTPJFUl+ZtK1SDvDMJGmRJK9gP2BlwGnTLgcaacYJtJOSvKmJL/dYfkVSe7v+WQ/AO0npZfT3OvpHW3fR7QPUfpuew8qaSoZJtJkfK2qjultSPJY4JHAPcADAFV1f9tv5DctlLowTKQhJDm9fYTp3wJHjGgz/w14K3AtcOSItiGNhGEizSPJv6e5lfgTgecBi35xvH0A2NOA82ieKnnUXP2labMU7hosTdrPAhdU1X0ASS7cMSPJL9E8SGo34K6qOr1f2xDb+O/AH1ZVJTFMtOQYJtJwHnJ77STHASur6tXt67OS/HyftiOq6oZBK05yDM2I5+lJ3gXsCVw9ih9CGhVPc0nzuxR4bvvOqn1oHsUKzeNtz5zV90192r47z/rfAjy7qlZU1QrgaByZaIlxZCLNo6quSHIecCVwC82T9QB2px2xJDmM5jMim2a3VdW/DFp3kmfQPInvMz3buy3JXkn2qwc/O1yaWj4cS1qgJD8FnA5spQmW36P5nMiD2qrq9lnLraB5hOoTdmJbG2lOn90+X19pEgwTacySHELz7O87Zn/WpE/fRwBfpHnu+E85UtG0MkwkSZ15AV6S1JlhIknqzDCRJHVmmEiSOjNMJEmdGSaSpM4ME0lSZ4aJJKmz/w9kewThHB8TDQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "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(r\"d$_{OO}$ [$\\AA$]\")\n", "plt.ylabel(\"Count\");" ] }, { "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.7.6" }, "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 }