{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Geometry Optimization\n", "This tutorial will demonstrate the use of PyBigDFT for performing geometry optimization. We will start with a simple example of a molecular system. Then we will move to an advanced example involving a slab of NaCl.\n", "\n", "## Molecule Example\n", "Let's begin with a simple example of an Aflatoxin B1 molecule." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "cano = \"CN1C=NC2=C1C(=O)N(C(=O)N2C)C\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can use openbabel to create a system from its cannonical smiles representation." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from openbabel.openbabel import OBMol, OBConversion, OBBuilder\n", "conv = OBConversion()\n", "conv.SetInFormat(\"can\")\n", "\n", "mol = OBMol()\n", "conv.ReadString(mol, cano)\n", "mol.AddHydrogens()\n", "\n", "builder = OBBuilder()\n", "builder.Build(mol)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then we convert from babel to PyBigDFT. An openbabel structure after calling build is usually not a good starting point. Normally you would do some combination of optimization and conformer search. For our case, we will do geometry optimization followed by a short molecular dynamics run. This will give us a decent geometry, but one that is not the actual minimum." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/wddawson/Documents/CEA/binaries/bds/install/lib/python3.7/site-packages/BigDFT/IO.py:371: UserWarning: Unsupported bond type had to be set to 1 (i.e. aromatic)\n", " UserWarning)\n", "/Users/wddawson/Documents/CEA/binaries/bds/install/lib/python3.7/site-packages/BigDFT/IO.py:371: UserWarning: Unsupported bond type had to be set to 1 (i.e. aromatic)\n", " UserWarning)\n" ] } ], "source": [ "from BigDFT.Interop.BabelInterop import convert_babel_to_system, molecular_dynamics, optimize_system\n", "sys_start = convert_babel_to_system(mol)\n", "sys_opt = optimize_system(sys_start)\n", "sys = molecular_dynamics(sys_opt, 10000, 300)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "application/3dmoljs_load.v0": "
\n

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n jupyter labextension install jupyterlab_3dmol

\n
\n", "text/html": [ "
\n", "

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n", " jupyter labextension install jupyterlab_3dmol

\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from BigDFT.Visualization import InlineVisualizer\n", "viz = InlineVisualizer(400, 300)\n", "viz.display_system(sys_start, colordict={x: \"black\" for x in sys_start}, show=False)\n", "viz.display_system(sys, colordict={x: \"blue\" for x in sys}, show=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we are ready to optimize the geometry of this system suing BigDFT. First, set up a basic calculator and input file." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "from BigDFT.Calculators import SystemCalculator\n", "code = SystemCalculator(verbose=False, skip=True)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "from BigDFT.Inputfiles import Inputfile\n", "inp = Inputfile()\n", "inp.set_xc(\"PBE\")\n", "inp.set_hgrid(\"0.37\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "BigDFT has a number of built-in geometry optimization methods which we can probe from the documentation." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Help on function optimize_geometry in module BigDFT.InputActions:\n", "\n", "optimize_geometry(inp, method='FIRE', nsteps=50, betax=4.0, frac_fluct=1.0, forcemax=0)\n", " Optimize the geometry of the system\n", " \n", " Args:\n", " nsteps (int): maximum number of atomic steps.\n", " method (str): Geometry optimizer. Available keys:\n", " * SDCG: A combination of Steepest Descent and Conjugate Gradient\n", " * VSSD: Variable Stepsize Steepest Descent method\n", " * LBFGS: Limited-memory BFGS\n", " * BFGS: Broyden-Fletcher-Goldfarb-Shanno\n", " * PBFGS: Same as BFGS with an initial Hessian obtained from a force\n", " field\n", " * DIIS: Direct inversion of iterative subspace\n", " * FIRE: Fast Inertial Relaxation Engine as described by Bitzek et\n", " al.\n", " * SBFGS: SQNM minimizer, keyword deprecated, will be replaced by\n", " SQNM in future release\n", " * SQNM: Stabilized quasi-Newton minimzer\n", " betax (float): the step size for the optimization method.\n", " This stepsize is system dependent and it has therefore to be\n", " determined for each system.\n", " frac_fluct (float): Fraction of force fluctuations. Stop if\n", " fmax < forces_fluct*frac_fluct.\n", " forcemax (float): Max forces criterion when stop.\n", "\n" ] } ], "source": [ "from BigDFT.InputActions import optimize_geometry\n", "help(optimize_geometry)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's activate the SQNM method. The important parameter to adjust is `betax` which is the step length, and the appropriate value can be sensitive to the method and system being computed. BigDFT automatically computes a measure of the force fluctations that come from the discretization and SCF error. This is used as a convergence measure. If you want to converge tighter, you should consider decreasing hgrid and the wavefunction convergence (`set_wavefunction_convergence`)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "inp.optimize_geometry(method=\"SQNM\", betax=1.0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And run." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Found 27 different runs\n" ] } ], "source": [ "log = code.run(input=inp, posinp=sys.get_posinp(), name=\"caf\", run_dir=\"work\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will also do a calculation using implicit solvent." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "from copy import deepcopy\n", "inp2 = deepcopy(inp)\n", "inp2.set_implicit_solvent()" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Found 55 different runs\n" ] } ], "source": [ "log_imp = code.run(input=inp2, posinp=sys.get_posinp(), name=\"caf-imp\", run_dir=\"work\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first thing we want to do is extract the energies at each step of the geometry optimization and plot them." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "energy = []\n", "forces = []\n", "for run in log:\n", " energy.append(run.energy)\n", " forces.append(run.forcemax)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdUAAAEJCAYAAADCRb2jAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de5xdVX338c+XDCFAuElAI9EAhSJ3CEnARx50AuFWBCMwQi0BW56QaLT2RqEaEsBWSb2UiCai0oLSYsTSRrkJZCxgBSdIAoSAxsglAiEBJEjkEvJ7/tjrJDsnZ2b2mTlnzsyZ7/v12q9z9t5r7bN2zmR+s9ZeF0UEZmZm1ntbNboAZmZmzcJB1czMrEYcVM3MzGrEQdXMzKxGHFTNzMxqpKXRBWgGW221VWy77baNLoaZ2YCybt26iIimqtw5qNbAtttuy6uvvtroYpiZDSiS/tDoMtRaU/2FYGZm1kgOqmZmZjXioGpmZlYjfqZqVmNvvvkmK1eu5LXXXmt0UfqVYcOGMWrUKLbeeutGF6XhBtvPyGD67uW5f3tv++23D3dUspLf/OY37LDDDuy6665IanRx+oWI4IUXXuCVV15hr732anRxGm4w/Yx09d1LWhcR2zeoaHXh5t8GmD17Nu3t7dnOrFkAtLe3M3v27MYVymrmtddeGxS/LKshiV133XXQ1My6M5h+Rgbbd++g2gDjxo2jra0tC6yXXkp7ezttbW2MGzeu0UWzGhkMvyyr5X+TzQ2mf4/BdK8Oqg3Q2trK/PnzaWtrA6CtrY358+fT2tra4JKZmVlvOKg2wqxZtE6YwOo1awBYvWYNrRMmbGwKNjOrtzlz5rD//vvz0Y9+tNFFaSoOqo0waxbtCxey24gRAOw2YgTtCxc6qA5Cmz1fT/rr8/X169c3ughWw98RX//617nlllu4/vrru03r7744B9UGKD1DnT9/PsDGpuDyX67W/DZ7vg41fb7+3e9+l/Hjx3PYYYdxwQUX8NZbbzF8+HA+85nPcOihh3LUUUexatUqAFavXs3pp5/OuHHjGDduHD/96U8BmDVrFlOmTOH4449n8uTJrF69mokTJzJmzBguuOACRo8ezZo1a5gxYwZXXnnlxs/+zGc+w5w5c3p9D1bm0ktrcpmpU6eyYsUKTj31VL70pS/xoQ99iEMOOYSjjjqKhx56CNjyu1+1ahWTJk3i0EMP5dBDD+V///d/gco/Z4NaRHjr5bbddttFNa644opYuHBhtjNzZkRELFy4MK644oqqrmP906OPPlpV+oULF8aIESNixowZMWLEiE0/G70swymnnBJvvPFGRERMmzYtrr322gBiwYIFERHxd3/3d3H55ZdHRMTZZ58d99xzT0REPPnkk/Ge97wnIiJmzpwZY8aMiTTxeXziE5+If/qnf4qIiFtvvTWAWL16dfzmN7+Jww8/PCIi3nrrrdh7771jzZo1Fctlvfh3gJqVYfTo0bF69eqYPn16zJo1KyIi7rrrrjj00EMjYsvvvq2tLb7yla9ERMT69evjd7/7Xac/Z5VUumfg1egHv8NruXnyhwa48MILN+2k5pzW1lZ3VBqkWltbmTZtGpdffjkzZsyoyc/BXXfdxQMPPLCxxvuHP/yB3XffnaFDh3LKKacAcMQRR3DHHXcAcOedd/Loo49uzL927VpeeeUVAE499VRKqzDde++93HTTTQCceOKJ7LLLLgDsueee7Lrrrjz44IOsWrWKww8/nF133bXX92FkvyPyNdRST9qZM2vSHHzvvffygx/8AIAJEybwwgsv8PLLLwObf/cLFy7kuuuuA2DIkCHstNNOfOc736n4c1ZLkk4ErgSGAN+KiC+UnVc6fzKwDjgvIn4haT/ge7mkewOXRMS/1LSAZRxUzRqsvb2duXPnMmPGDObOnVuTP7AignPPPZfPf/7zmx3/4he/uHF4w5AhQzY+K9uwYQM/+9nPqLSE4fbbbxqbn1UuKjv//PP5t3/7N5577jn+/M//vFflt5xZszYFTwm6+A56otJ3WvoZyX/3neWt9HNWK5KGAF8DJgIrgQ5JCyLi0Vyyk4B903YkMBc4MiIeBw7LXee3wE11KWiOn6maNVD++fpll11Ws+frxx57LDfeeCPPP/88AC+++CJPPvlkp+mPP/54rrrqqo37ixcvrpju6KOP3tgX4Mc//jEvvfTSxnOTJk3itttuo6OjgxNOOKFX5be+c8wxx2zsrPSTn/yEESNGsOOOO26R7thjj2Xu3LkAvPXWW6xdu7bqn7MeGA8sj4gVEfEGcANwWlma04DrUovyfcDOkkaWFx/4dUTUtHCVOKiaNVBHR8dmY5RLY5g7Ojp6dd0DDjiAz33ucxx//PEccsghTJw4kWeffbbT9HPmzGHRokUccsghHHDAAcybN69iupkzZ/LjH/+YMWPGcOuttzJy5Eh22GEHAIYOHUprayttbW0MGTKkV+W3TsycWfNLzpo1a+N3f9FFF3HttddWTHfllVfS3t7OwQcfzBFHHMHSpUur/jmroEXSotw2pez8HsDTuf2V6Vi1ac4C/qOagvWU5/6tAc/9a3nLli1j//33b3Qx6uL1119nyJAhtLS08LOf/Yxp06ZtrNVu2LCBMWPG8P3vf5999923Yv5m/repxmD8d6h0z93N/SvpTOCEiDg/7Z8DjI+IT+bS3Ax8PiLuTft3ARdGxANpfyjwDHBgRKyq8W1twc9Uzaywp556ira2NjZs2MDQoUP55je/CcCjjz7KKaecwqRJkzoNqGY9sBJ4V25/FFmArCbNScAv+iKggoOqmVVh33335cEHH9zi+AEHHMCKFSsaUCJrch3AvpL2IutodBbwp2VpFgDTJd1A1lHp5YjIt0GfTR81/YKDqlldRMSgmkS8CD9q2txg+hnp6XcfEeslTQduJxtSc01ELJU0NZ2fB9xCNpxmOdmQmo+V8kvajqzn8AW9uoEqOKia1diwYcN44YUXBs3SXkVEZGtqDhs2rNFF6RcG089Ib7/7iLiFLHDmj83LvQ/gE53kXQf06YBpd1SqAXdUsrw333yTlStXDpr1I4saNmwYo0aNYuutt250URpusP2MdPbdN+Mi5Q6qNeCgamZWvWYMqh6namZmViMND6qSzpS0VNIGSWNzxydKekDSw+l1Qu7cbZKWpHzz0hRU5dcdL2lx2pZImpSObyfpZkmPpfxfyOXZRtL3JC2XdL+kPet792Zm1kwaHlSBR4APA3eXHV8DfDAiDgbOBb6TO9cWEYcCBwG7AWd2ct2xEXEYcCLwDUmljllfjIj3AIcD75N0Ujr+F8BLEbEP8BXgil7fnZmZDRoND6oRsSxNfFx+/MGIKA3gXQoMk7RNOrc2HW8BhgJbPBiOiHURUVpZd1gpTTrent6/AfyCbLAwZHNIlubouhE4Vs3eNc/MzGqm4UG1oNOBByPi9dIBSbcDzwOvkAXALUg6UtJS4GFgai7Ils7vDHwQuCsd2jiHZEr7Mp10x5Y0pTRfZWmlDzMzG9z6JKhKulPSIxW28tUGKuU9kKwZdrPBuxFxAjAS2AaYUCErEXF/RBwIjAMulrRxoFRqCv4PYE5ElKaCqVQrrdg9OiKujoixETG2pcXDfc3MrI8mf4iI43qST9IosvXvJkfErytc9zVJC8iabe/o4vOXSXqV7BnsonT4auBXZQvWluaQXJmC7k7Aiz0pu5mZDT79tvk3Nc3eDFwcET/NHR9eWisvBb6Tgccq5N+r1DFJ0mhgP+CJtP85soD56bJsC8g6RQGcASwMD+Q1M7OCGj75Qxrq8lWyXry/AxZHxAmSPgtcDPwql/x4sibaH5E1+w4BFgJ/leaIPJWsx+8laYmgi4A3gQ3AZRHxX6n2+zRZIC49o70qIr6Vmoe/Q9Yr+EXgrFzTcKc8+YOZWfWacfKHhgfVZuCgamZWvWYMqv22+dfMzGygcVA1MzOrEQdVMzOzGnFQNTMzqxEHVTMzsxpxUDUzM6sRB1UzM+u3JJ0o6fG0JOdFFc5L0px0/iFJY3LndpZ0Y1rqc5mk99a7vA6qZmbWL6W1sr8GnAQcAJwt6YCyZCcB+6ZtCjA3d+5K4La01OehwLJ6l9lB1czM+qvxwPKIWJGW6ryBbK73vNOA6yJzH7CzpJGSdgSOAb4N2VKfEfG7ehfYQdXMzBqlpbSEZtqmlJ3fuBxnsjIdK5Jmb2A18K+SHpT0LUl1n73JQdXMzBplfWkJzbRdXXa+yHKcnaVpAcYAcyPicOBVsvng68pB1czM+qvScpwlo4BnCqZZCayMiPvT8RvJgmxdOaiamVl/1QHsm5byHAqcRbZEZ94CYHLqBXwU8HJEPBsRzwFPS9ovpTsWeLTeBe6TRcrNzMyqlZb0nA7cTrbU5zURsVTS1HR+HnAL2bray4F1wMdyl/gkcH0KyCvKztWFl36rAS/9ZmZWPS/9ZmZmZp1yUDUzM6sRB1UzM7MacVA1MzOrEQdVMzOzGml4UJV0pqSlkjZIGps7PlHSA5IeTq8Tcuduk7Qk5ZuXJl0uv+54SYvTtkTSpHR8O0k3p1ULlkr6Qi7PeZJW5/KdX+/7NzOz5tEfxqk+AnwY+EbZ8TXAByPiGUkHkY1TKs352BYRayWJbJaMM8kmWi6/7tg0zmkksETSD9O5L0ZEexq7dJekkyLi1nTuexExvba3aGZmg0HDg2pELAPI4uNmxx/M7S4FhknaJiJej4i16XgLMJQt54IkItbldoeV0qTj7en9G5J+QTatlZmZWa80vPm3oNOBByPi9dIBSbcDzwOvkNVWtyDpSElLgYeBqRGxvuz8zsAHgbvyn5UWur1RUn4+yfJrTymtrLB+/frOkpmZ2SDSJ0FV0p2SHqmwla+LVynvgcAVwAX54xFxAjAS2AaYUCErEXF/RBwIjAMuljQsd90W4D+AORGxIh3+IbBnRBwC3Alc21m5IuLq0soKLS0Nr/CbmVk/0CfRICKO60k+SaOAm4DJEfHrCtd9TdICskVq7+ji85dJehU4CFiUDl8N/Coi/iWX7oVctm+SBXMzM7NCugyqqTZ3KvAnwKHAzsDvgCXArcB/lTep1kpqmr0ZuDgifpo7PhzYISKeTeU7GbinQv69gKdTR6XRwH7AE+nc54CdgPPL8oyMiGfT7qnAsprfmJmZNa1OJ9SXdAHwGbLA8j/p9RVgB2B/4P3p9Z/SSgE9K0A21OWrwG5kAXtxRJwg6bPAxcCvcsmPJ1uQ9kdkzb5DgIXAX6XgeSpZj99LJJ1DtiDtm8AG4LKI+K9U+30aeAwoPaO9KiK+JenzZMF0PfAiMC0iHuvuHjyhvplZ9ZpxQv2uguqXgH9Oa9J1lmYk8DcR8bd1Kt+A4KBqZla9QRVUrTgHVTOz6jVjUC3cUUnS/sAZwDsi4hOS3gMMjYiH6lY6MzOzAaTQkBpJZ5I9V90DOCcdHg58uU7lMjMzG3AKNf9KWgacHRGLJb0UEbtI2hp4JiJ2q3sp+zk3/5qZVa8Zm3+LTv6wO9kwGtg0JWBQYXpAMzOzwapoUH2ATc2+JWcBP69tcczMzDaRdKKkxyUtl3RRhfOSNCedf0jSmNy5J9JKZ4slLSrPWw9FOyp9CvixpL8Atk/z7v4x2bhRMzOzmkvLen4NmAisBDokLYiIR3PJTgL2TduRwNz0WtIaEWv6qMjFgmpEPJZ6+55CNvHC08CPIuL39SycmZkNauOB5aX52SXdQDYtbT6ongZcF1kHofsk7Vw2O16fKjyhfloy7WfAPRFxgwOqmZn1Uktpta+0TSk7vwdZJa5kJZvW1S6SJshaWR+ocO26KFRTlfRushVdDiMr5HBJZwAnRsT5XWY2MzOrbH1EjO3ivCocK+8g21Wa90XEM5J2B+6Q9FhE3N2TghZVtKb6DbLJ7Xcgm0sXslVhJtajUGZmZmS1zvy61qOAZ4qmiYjS6/NkK56Nr1tJk6JBdTzwhYjYQPoLICJeJlvpxczMrB46gH0l7SVpKNmokwVlaRYAk1Mv4KOAl9MqZttL2gFA0vZkHWsfqXeBi/b+XQXsA/yydEDSAcBT9SiUmZlZWn1sOnA72apk10TEUklT0/l5wC1kS4AuB9YBH0vZ3w7cJAmyWPfvEXFbvctcdEalPydbRu3zwJXABcA/kNVer69rCQcAz6hkZla9/jKjUpoh8Ci2XDf8voh4s6u8W1yr6Co1kj4ETAFGk9VQvxER/1XNhzUrB1Uzs+o1OqhKGkFWYTyXbA3tx9h83fBdgGvJKpCFxrp22/ybBt/eBZzgIGpmZk3kHuDbwGER8dvyk5LeCXwUuBs4oMgFizb/Pgm8JyL+UFVxBwnXVM3MqtcPaqpDI+KNWqWD4r1/LwXmShotaYikrUpbwfxmZmb9StFAWTQdFK+pbihdO384+6wYUvTDmpVrqmZm1Wt0TbUISV+PiI8XTV90SM1ePSyPmZnZQFZpxqZOFW2+PTMinizfgNOrL9/mJJ0paamkDZLG5o5PTPM1PpxeJ+TO3SZpSco3L3WmKr/u+LTcz+KUdlJ3+SVtI+l7aQmh+yXt2dv7MzOzgSsiplWTvmjz79qI2LHC8Rcj4m3VfGCFa+wPbCCbCvFvI2JROn44sCrN23gQcHtE7JHO7RgRa5WN6r0R+H5E3FB23e2AN9Lg4ZFkY47emfYr5pf0ceCQiJgq6SxgUkR8pLt7cPOvmVn1+kPzr6S9Kxx+E/htmkWwKl02/+Zqh0MktbJ5NXhvsvE8vRIRy9JnlR9/MLe7FBgmaZuIeD0i1qbjLcBQtpxgubSqTsmwfJou8p8GzErvbwSukqQoOpjXzMwGmuVkMaC8mXe9pO8DH0/T8hbS3TPVb6fXYcA1ueMBPAd8sugH9dLpwIMR8XrpQFoofTxwK1kA3IKkI8nKPRo4JyLWd5N/4xJCqUb7MrArsMWg37SM0BSAoUOH9vL2zMysESJii8egklrIKo7/CPwz6Xd9EV0G1YjYK33A9RHx0eqKulkB7wTeUeHUZyLiv7vJeyBwBdlkyPmynSBpGHA9MIFs1RzK0twPHJiamK+VdGtEvNZF/iLLDJWufTVwNWTNv13dg5mZDRypAvZLSRcAD1WTt+iMSh8uNb32sIDH9SSfpFFky/VMjohfV7jua5IWkDXbbhFUc+mWSXoVOAhY1EX+0hJCK9NfKjuRTV1lZmaDz1pgu2oydNv7NyLeIludZtceFqpHJO1MtobrxRHx09zx4anjUamKfjLZfI3l+fdK55E0GtgPeKKb/AvI5oAEOANY6OepZmaD1kfI+vQUVnSc6vXAjyRdSVaby3f6WVjNB5ZLQ12+CuwG3CxpcUScAEwnW25uhqQZKfnxZE20CyRtQ7YU0EJgXrrWqcDYiLgEOBq4SNKbZL2LPx4RayS9vbP8ZM+QvyNpOVkN9aze3JuZmfVvkr7Dlo/5tgb2JKuMnVzV9QoOqflNJ6ciIip1Rx5UPKTGzKx6/WRIzcwKh9eTrcZ2W0Ssrup6bt3sPQdVM7Pq9YegWmueEN/MzKwTkt5XTfpCz1Ql7Ug2KcL7gRHkhp5ExLur+UAzM7MB5DayRcsLKVpT/TowBrgMeBvZpA9PAV+ptnRmZmZFSTpR0uNpTvaLKpyXpDnp/EOSxpSdHyLpQUk/6snnR0ThgArFg+rxwOlpooa30utHgHOqLJ+ZmVkhaZ6ErwEnAQcAZ0s6oCzZScC+aZsCzC07/5fAsh589laS/kTS/GryFQ2qWwGluQ9/n8aQPks25MXMzKwexgPLI2JFWij8BrLJevJOA66LzH3Azrm5CEYBfwJ8q+gHSjpU0peB3wL/DjxfTYGLjlNdQvY89S7gHrK/HH5PNimEmZlZPWycjz1ZCRxZIM0eZBW/fwEupJtnomn+go+STf5zAHA3MBw4OCKeqKbARWuq/w8oXfhTwB+AnYHJ1XyYmZlZToukRbmtfOL6IvOxV0wj6RTg+Yh4oKsCpGetTwN/ClwLvDsijiWrOK7rKm8lhWqqEbEi9341cH61H2RmZlZmfUSM7eJ8aT72klHAMwXTnAGcKulkspXWdpT03Yj4s7L8HyCb4/dW4JaIeLbqu8gpup5qp3o7TaGZmVknOoB9Je1F9ozzLLIaZd4CYLqkG8iahl9OgfHitCHpA8DfVgioALuTBeBzgX+QtIRsat6t6WSVsq4UXU+15F1s3nYdZGvOmZmZ1VRa13o6cDvZXO3XRMRSSVPT+XnALWTz8y4na679WJWfsQ64DrhO0rvJHmtOIRs++h1JcyLilqLXq2qaQkkvRcQu1RR4MPA0hWZm1evP0xRKei9wHnBGRBRepa3aoPpiRLyt+uI1NwdVM7Pq9eegWlLtWuKe+9fMzAYlSV+W9I5uku2Sxq0WUnScqpmZWbN5HPi5pGXA/6T9V8jGtf4xWc/g/YDPFb1gl82/ku5h895P7wV+lk8TEccU/bBm5eZfM7Pq9YfmX0lbk83KdBJwMNkcDC8BD5F1gvphRKwvfL1uguq53V0gIq4t+mHNykHVzKx6/SGo1poXKa8BB1Uzs+o1Y1DttKOSpFOLXKBoOjMzs2bXVe/fsyQ9IuliSf9H0q6ShqbX90q6SNIjQFtfFdbMzKw/6+6Z6sHABWQPcPdiU6elX5M9wP1mRCytdyH7Ozf/mplVb1A1/wJExMMRMT0i/ohsGZx3ATtExB9HxKdrEVAlnSlpqaQNksbmjk+U9ICkh9PrhNy52yQtSfnmpYVsy687XtLitC2RNKm7/JLOk7Q6l88LB5iZDRKS9pc0Q9LX0v57JB1SzTUKT/4QEesi4pk0T2ItPQJ8mGz9urw1wAcj4mCyiY6/kzvXFhGHAgcBuwFndnLdsRFxGHAi8A1JLQXyfy8iDktb4YVt62n27Nm0t7dvdqy9vZ3Zs2c3qERmZs1F0plkY1X3AM5Jh4cDhSd+gH4wo1JELIuIxyscfzAiSkv8LAWGSdomnVubjrcAQ6mwkkD6I6A0tmhYPk2R/P3JuHHjaGtr2xhY29vbaWtrY9y4cQ0umZlZ07gMOD4ipgJvpWNLgEOruUjDg2pBpwMP5udflHQ78DzZ7Bc3Vsok6UhJS4GHgan5Abxd5D9d0kOSbpT0LjohaUppYd316wuPC+6R1tZW5s+fT1tbG5dccgltbW3Mnz+f1tbWun6umdkgsjtZEIVNFa2gykpXnwRVSXemnsTl22kF8h4IXEHWYWqjiDgBGAlsA1Rc9zUi7o+IA4FxwMWShnWT/4fAnhFxCHAn2SrwFUXE1RExNiLGtrTUf7bH1tZWpk2bxuWXX860adMcUM3MausBNjX7lpwF/LyaixQKqpI+JWlENRfOi4jjIuKgCtt/d/O5o4CbgMkR8esK132NbIHaLoNzRCwDXiV7htpp/oh4IVcb/iZwRLE7rL/29nbmzp3LwmOOYe7cuVs8YzUzs175FPA5Sf8DbJ9aMy8H/qqaixStqR4HPCHpR5I+Unq2WU+SdgZuBi6OiJ/mjg+XNDK9byFbnPaxCvn3KnVMkjSabFLkJ7rKXzqenAosq8e9Vav0DHX+/Pm03n33xqZgB1Yzs9qIiMeA9wBfAz4L/CtwcET8qprrFAqqEXEqMBq4Ffg08Jykb0nq9WT6kiZJWkk2Wf/N6a8DgOnAPsCM3BCX3YHtgQWSHiJr/34emJeudaqky1L+o4ElkhaT1XY/HhFrusoPfCoNs1lC9lfLeb29v1ro6OjY7Blq6RlrR0dHg0tmZtYcJO0BbBMR8yPinyPiBmBrSe+s6jo9mfs3jdv5Dllz6tNkTaVXRsTvq75YE6j75A+zZsGll255fObM7JyZ2QDUnyZ/kNQB/HlEPJw7djDwrYg4svB1qgmqko4F/ozsGeQiso48TwF/Cbw9Iv5v4Ys1kT6dUUkCL4JgZk2gnwXVlyNip6LHO1O0o9IXUxPtHLLnjwdHxPERcX1E3AOcDRxe9EPNzMyKkHSipMclLZd0UYXzkjQnnX9I0ph0fJikn+dmz6vQ3LeZ1ZL2Kbv2PsAL1ZS36FiQYcCkiKj4EC8i3sxPMWh1NHNmo0tgZtYn0hSyXwMmAiuBDkkLIuLRXLKTgH3TdiQwN72+DkyIiN8rW4j8Xkm3RsR9nXzcNcAPJH0GWAH8EVnv36pm1isUVCNieoE0W/TAtTrwM1QzGzzGA8sjYgWApBvIHj/mg+ppwHWRPcu8T9LOkkZGxLNAqZ/P1mnr6tnZF4A3gS+SzXP/FPBtqpymsFBQlXRPJ4V5neyvh/+MiB9W88FmZjbotUhalNu/OiKuzu3vQdYZtmQlWS2UbtLsATybaroPkI0k+VpE3F+pECndNcCUiPjnHt1JUnSc6k+APckmG/5ueh1N1llpFXCNpAt7UxAzMxt01pdmpkvb1WXnVSFPeQWv0zQR8VZaVGUUMF7SQRXSEhFvAccDG6or/paKPlM9HjghzUwEgKTrgWsj4khJ/wncAHjZFDMzq5WVZE2xJaOAZ6pNExG/k/QTshXLHunks74CXCppVkS80dMCF62pvofswW3ek2SzFBERPyebjNjMzKxWOoB90wx5Q8nm4l1QlmYBMDn1Aj4KeDkinpW0W5qZD0nbks0M2FXfn08CfweslfS0pKdKWzUFLlpTvRv4V0mXkP1VMAqYBdybCnww8Gw1H2xmZtaViFgvaTpwOzAEuCYilkqams7PA24hm252ObAO+FjKPhK4Nj0v3QqYHxE/6uLj/qwWZS40+YOktwFfJ1tMfAiwHvhP4JMRsUbSfsAOEbGoi8s0rT6d/MHMrEn0p8kfaqXb5t8U5T9NNg/uMOCdwLYRcXaaS5eIeHywBlQzMxv4JG0t6VJJKyS9ll4vTc3OhXUbVFOvqE8Ab0TEhohYFRG97iFlZmbWj8wme+46FTg0vU4gW8+7sKLNv18mG4D79erL2fzc/GtmVr3+1PybpuI9NCJeyB0bASyJiD2KXqdoR6XxwCfTWNSnyY0TioheL/9mZmbWYJXGu3Z1vKKiQfWbaTMzM2tG3wd+mCbef4psgqPPAvOruUiP1lO1zbn518ysev2s+XcoWRD9U7IOub8lm9TocxHxeuHrFHymKuB8siXeRkTEITd1dQkAABZuSURBVJKOAd4REVVF8WbkoGpmVr3+EFQlvSMinqvV9YrOqHQZ8BfA1cC707GVwN/XqiBmZmYN8Mv8Tpp2t8eK1lSfBg5PEz28FBG7pNrrixGxS28K0AxcUzUzq14/qam+EhE75PZfjIi39fR6RWuqQ9i0Ll0pCg/PHTMzMxuIatqxqGhQvQX4sqRtYOMz1suBXq+hKulMSUslbZA0Nnd8oqQHJD2cXifkzt0maUnKNy/N+lR+3fGSFqdtiaRJFdIskPRIbn8bSd+TtFzS/ZL27O39mZlZv9YiqVXShBRnNtvPx54iijb/7ghcR7ZsztbAa8CPgckR8Ur197DZtfcnW8PuG8DflqY7lHQ4sCoinklr4N1eGoAraceIWJuC+43A9yPihrLrbkc2C9R6SSOBJcA7I2J9Ov9h4AzgkIg4KB37eNqfKuksYFJEfKS7e3Dzr5lZ9fpJ8+8TdF1bjYjYu+j1Co1TjYi1wIck7U42dufpWvWWKq3RmsXHzY4/mNtdCgyTtE1EvJ7KA1n5h1LhHyQi1uV2h+XTSBoO/DUwhc3HIJ1GtvoOZMH6KkkKjzsyM2tKEbFnLa9XtPk37wVgO0l7SyocvXvpdODB/FghSbcDzwOvkAXALUg6UtJS4GFgaqmWStZ0/SWyZYLy9iCbMYqU9mVg106uPUXSIkmL1q9fXymJmZkNMoWCqqQTJf0WeI5szbrS9quC+e+U9EiF7bQCeQ8km9D4gvzxiDiBbL28bcgmPd5CRNwfEQcC44CLJQ2TdBiwT0TcVOnjKl2mk2tfHRFjI2JsS0vRianMzKyZFY0GXyOr3V0bEX+o9kMi4rhq8wBIGgXcRPbs9tcVrvuapAVkzbZ3dPH5yyS9ChxEFmCPSO3oLcDukn4SER8gG3v7LmClpBZgJ+DFnpTdzMwGn6LNv7sA3+hJQO0pSTsDNwMXR8RPc8eHp45HpMB3MvBYhfx7pfNIGg3sBzwREXMj4p2pHf1o4JcpoAIsAM5N788AFvp5qpmZFVU0qH4b+Fg9CiBpUlpy573AzelZKcB0YB9gRm5ozO7A9sACSQ+R9eh9HpiXrnWqpMtS/qOBJZIWk9V2P15aVL0L3wZ2lbScrCPTRbW7UzMza3ZFh9TcQ7b825Nkz1U38tJvHlJjZtYTRYbUSDoRuJJsEqJvRcQXys4rnT+ZrPPpeRHxC0nvIhsK+g6yYZtXR8SVdbiNzRR9pvqttNkAMHv2bMaNG0dra+vGY+3t7XR0dHDhhRc2sGRmZsWliX2+Bkwk6/PSIWlBRDyaS3YSsG/ajgTmptf1wN+kALsD8ICkO8ry1lzRcarX1rMQVlvjxo2jra2N+fPn09raSnt7+8Z9M7MBZDywPCJWAEi6gaxjaj4wngZcl/q/3CdpZ0kjI+JZ4FmAiHhF0jKyYZN1DapdPlOVNKds/y/K9n9Qj0JZ77S2tjJ//nza2tq45JJLNguwZmYDyMa5A5KV6VhVadKUs4cD99e8hGW666h0Xtn+P5ftT6xdUayWWltbmTZtGpdffjnTpk1zQDWz/qilNIlO2qaUnS8yd0CXadIMej8APp2bja9uumv+LS9spcJbP9Te3s7cuXOZMWMGc+fOpbW11YHVzPqb9RExtovzpbkDSkYBzxRNI2lrsoB6fUT0ap3UorqrqZb/ReAxmwNA/hnqZZddtrEpuL29vdFFMzOrRgewb5p3YChwFtl8AnkLgMnKHAW8HBHPpl7B3waWRcSX+6rA3dVUWyS1sqmGWr6/xZJr1ngdHR2bPUMtPWPt6OhwbdXMBoy0yth04HayeHNNRCyVNDWdn0e2NOnJZFPnrmPTnArvA84BHk7zFQD8Q0TcUs8ydzlOtcCSOETEXjUu04DjcapmZtXrD0u/1VqXNdVaL4ljZmbWzHqy9JuZmZlV4KBqZmZWIw6qzW7WrEaXwMxs0Cg0ob51rV93VJLA37GZ9UPN2FHJNVUzM7MacVBtRrNmZTVUpeHEpfduCjYzqys3/9aAm3/NzKrn5l9rOrNnz95i+sL29nZmz57doBKZmQ1cDqrNbubMLk+X1l4tBdbSvMHjxo3ri9KZmTUVN//WQL9u/i2gFEinTZvG3LlzvfaqmfUJN/9aU/Laq2ZmtdHwoCrpTElLJW2QNDZ3fKKkByQ9nF4n5M7dJmlJyjdP0har5UgaL2lx2pZImlQhzQJJj+T2z5O0Opfv/Hrcc39Tvvaql4gzM+uZ7pZ+6wuPAB8GvlF2fA3wwYh4RtJBZEv/7JHOtUXE2rRe3o3AmcANFa47Ni0dNBJYIumHEbEeQNKHgd9XKM/3ImJ6Te5sAMivvVpayDy/b2ZmxTW8phoRyyLi8QrHH4yI0grvS4FhkrZJ59am4y3AUCosTxcR60oBFBiWTyNpOPDXwOdqdiMDVFdrr5qZWXX6Q021iNOBByPi9dIBSbcD44FbyWqrW5B0JHANMBo4JxdkLwe+RLag7RafJekY4JfAX0XE051cewowBWDo0KE9uad+4cILL9ziWKnGamZm1emTmqqkOyU9UmE7rUDeA4ErgAvyxyPiBGAksA0woUJWIuL+iDgQGAdcLGmYpMOAfSLipgpZfgjsGRGHAHcC13ZWroi4OiLGRsTYlpaB8reJmZnVU58E1Yg4LiIOqrD9d1f5JI0CbgImR8SvK1z3NWAB0GVwjohlwKvAQcB7gSMkPQHcC/yxpJ+kdC/kasPfBI6o5j4HPE9jaGbWKw1/ptoZSTsDNwMXR8RPc8eHp45HSGoBTgYeq5B/r3QeSaOB/YAnImJuRLwzIvYEjgZ+GREfSOlG5i5xKrCsHvfWb116aaNLYGY2oDU8qEqaJGklWQ3y5vSsFGA6sA8wIzfEZXdge2CBpIeAJcDzwLx0rVMlXZbyH03W43cxWW334xGxppvifCoN01kCfAo4r3Z3amZm1ZJ0oqTHJS2XdFGF85I0J51/SNKY3LlrJD2fHzpZ9/J6RqXeG9AzKs2aVbmGOnNmxebg2bNnM27cuM06MrW3t9PR0VGx05OZWWe6m1EpzUHwS2AisBLoAM6OiEdzaU4GPknWankkcGVEHJnOHUM2dPK6iDiobjeS0/CaqjXYrFnZKjalP65K7zt5vuq5gs2sD40HlkfEioh4g2w+gvI+NKeRBc2IiPuAnUuP8iLibuDFviywg6pVpTSOta2tjUsuucQTRZhZb7RIWpTbppSd3wPID2tcyaZJgKpJ02ccVG2Tbla0KfFcwdZj7mFum1tfGpqYtqvLzqtCnvJnlkXS9BkHVduk4C88zxVsPeYe5ladlcC7cvujgGd6kKbPOKhaVfJzBV+21VYbm4IdWM2sDjqAfdMQyaHAWWRzE+QtACanXsBHAS9HxLN9XdASB1WrymZzBV96abdzBc+ePXuLgNve3s7s2bP7orjWH8yaBVK2wab3bgq2bqSpZaeTLaiyDJgfEUslTZU0NSW7BVgBLCebtOfjpfyS/gP4GbCfpJWS/qLeZfaQmhoY0ENqekPa1Gu4E+Wr4JTv2yBT4GfGBg8vUm5WZa3DvYXNbDBxTbUGXFPt3iWXXMJWl1/OhhkzuOyyy7rPYM1p1iw3+9pGzVhTdVCtAQfVrpWafFevWcNuI0a4pmpmQHMGVTf/Ws8VGNeaf4YKdNtb2B2bzGwgc021BgZtTbWAn06cyPvuvHPL48cdx/vuuGOL45t1ZPqf/6H9/e/3c1izJtWMNVUH1RpwUC3IzcVWDT9/bXrNGFTd/Gv9TmkaRKDLaRDdVNzkPPuSDUAOqtZ3iswtnIbsXHb55QDZaydDdrxijpn1Nw6q1ncKNOW1v//97DZiBO0LF2b7Cxdm++9//xZp82Ngizx7dc22DmrdPOvZl2yAc1C1fmWzaRCh22kQS03FrXff3e2KOa7Z1kGtm2irXN93i7xmjRYR3nq5bbfddmF1MHNmt0kWLlwYI0aMiIAYMWJELFy4sFD6hccc0236K664YtP5VJaFCxfGFVdcUfQOmh/0n2vXsyxWF8Cr0Q9+h9dya3gBmmFzUG2MFZMnR2yqy2zcVkye3GW+GTNmRED22oWNAXjhwgjYfL/MoArAM2dW/Hcv8kdQ1Z9TDQfVAcdB1VvFzUG1MTYLZOkXaneBrKc12+7SVxOAt1DrYNSXGh3I+irAW104qHqruDmo9gMFfrlXXbOt8hd2NQG72j8I+m1NuNFBNa/aslQTeB2k68JBtR4FgDOBpcAGYGzu+ETgAeDh9Dohd+42YEnKNw8YUuG644HFaVsCTMqd+wnweO787un4NsD3yNblux/Ys8g9OKj2AwV+6dW1ZtvDAFy0ZlvPpujN0uc+r1DA7k/Bpp7PYPvTHw9NxEG1PkF1f2C/FOjyQfVw4J3p/UHAb3PndkyvAn4AnFXhutsBLen9SOD53P5mn5XL83FgXnp/FvC9IvfgoDoAFfgl2dPAV6hpuQfNlvVqit7s/MyZ3aavNgj3KmhXo57PYOtZCx7EHFTrG1wrBrrYFDxfALYpO7418EPgI91cey9gVYGgejvw3vS+BVhDmsqxq81BdQCqtmbbTY2vp89Ui3aaqmdTdLXpy++tqqDdTfp6B+x7jzuu4r/jvccd16u0WxhMNdte/AHhoFrPgnQdVM8A7iw7djvwEvDvlZp/U5ojUxPx7ys0/z6cmn5nsGkO5EeAUbl0vwZGdHLtKcAiYNHQoUMr/sDY4NGTZ551C3w97LxTOMDnyjJjxoyqyt5d+noG7C3Od/PHT7V/KNX7OXk16et57S3SF3yUUomDas8D5p0pYJVvp+XSdFZ7PDAFtz+qcG5Yav6d2M3n7w/8HBiW9vdIrzsAPwYmp/2lFYLqrt3dn2uqVq1qm1z7W1N0RBaEKRiEq0lfr4Bdnr6aWnmRtNV2hKvnc/V6Xrsn6TtTJKgCJ5L1gVkOXFThvIA56fxDwJiieeuxNTyq525+i6AKjAJ+Cbyvi3znAlcVuH57J0H7vFJ+N/9aX+lVM2eNm6J78guy3oGvXgE7n35mwfTVpK1ns3u16etalhoNZeouqAJDUuVmb2AoWafTA8rSnAzcmoLrUcD9RfPWY2t4MM39w2wWVIGd0z/C6WXphgMjY1Pg+x4wvcL19mLTM9TRwDPAiJRnRDq+NXAjMDXtf4LNOyrNL1J2B1XrT+rdlNdXTbT1rqkWSV/ttSPq95y8qvT1vHZP7rUTBYLqe4Hbc/sXAxeXpfkGcHZu/3Gyzqnd5q3H1h+C6SRgJfA6WWei29PxzwKvsmnYy2Jgd+DtQEeq5i8FvpoLnqcCl6X356Tzi4FfAB9Kx7cnG6JTyn8l6ZksWXPy98maCn4O7F3kHhxUrWlU24ErqVVnoj59ptpN+mqvvVmaAtNg5tMPuJpqD9JXUiCongF8K7d/DmUtk8CPgKNz+3cBY4vkrcfW8KDaDJuDqllt1Lv3bzXpq712rwL8QH6mWqBPQGdSZWpRbpsSmwfMMysExq+Wpbm5QlA9okjeemwND0jNsDmomlk9n5NXm75Pe/8mXaXvTDM2/5aGklgvbL/99vHqq682uhhmZgOKpHURsX0X51vIOqseC/yW7NHfn0bE0lyaPwGmk3VYOhKYExHji+Sth5Z6XtzMzKynImK9pOlkIzOGANdExFJJU9P5ecAtZAF1ObAO+FhXeetdZtdUa8A1VTOz6nVXUx2Itmp0AczMzJqFg6qZmVmNOKiamZnViJ+p1oCkDcAfepi9BVhfw+L0Z4PlXgfLfcLgudfBcp/Qt/e6bUQ0VeXOQbXBJC2KiLGNLkdfGCz3OljuEwbPvQ6W+4TBda/10FR/IZiZmTWSg6qZmVmNOKg23tWNLkAfGiz3OljuEwbPvQ6W+4TBda8152eqZmZmNeKaqpmZWY04qJqZmdWIg2oDSTpR0uOSlku6qNHlqRdJT0h6WNJiSYsaXZ5aknSNpOclPZI79jZJd0j6VXrdpZFlrJVO7nWWpN+m73axpJMbWcZakPQuSe2SlklaKukv0/Gm+l67uM+m+077kp+pNoikIWTLEk0EVpItS3R2RDza0ILVgaQngLERsabRZak1SccAvweui4iD0rHZwIsR8YX0x9IuEfH3jSxnLXRyr7OA30fEFxtZtlqSNBIYGRG/kLQD8ADwIeA8muh77eI+22iy77QvuabaOOOB5RGxIiLeAG4ATmtwmaxKEXE38GLZ4dOAa9P7a8l+UQ14ndxr04mIZyPiF+n9K8AyYA+a7Hvt4j6tFxxUG2cP4Onc/kqa9wc6gB9LekDSlEYXpg+8PSKehewXF7B7g8tTb9MlPZSahwd0k2g5SXsChwP308Tfa9l9QhN/p/XmoNo4qnCsWdvi3xcRY4CTgE+kZkRrDnOBPwIOA54FvtTY4tSOpOHAD4BPR8TaRpenXircZ9N+p33BQbVxVgLvyu2PAp5pUFnqKiKeSa/PAzeRNX03s1XpeVXpudXzDS5P3UTEqoh4KyI2AN+kSb5bSVuTBZrrI+I/0+Gm+14r3Wezfqd9xUG1cTqAfSXtJWkocBawoMFlqjlJ26dOEEjaHjgeeKTrXAPeAuDc9P5c4L8bWJa6KgWZZBJN8N1KEvBtYFlEfDl3qqm+187usxm/077k3r8NlLqq/wswBLgmIv6xwUWqOUl7k9VOIVtS6t+b6T4l/QfwAWAEsAqYCfwXMB94N/AUcGZEDPgOPp3c6wfImgkDeAK4oPTccaCSdDRwD/AwsCEd/gey541N8712cZ9n02TfaV9yUDUzM6sRN/+amZnViIOqmZlZjTiompmZ1YiDqpmZWY04qJqZmdWIg6rZICbp92nYk5nVgIOqWQOlZfGOk3SepHvr/Fk/kXR+/lhEDI+IFfX8XLPBxEHVrAlIaml0GczMQdWsP9gfmAe8NzXH/g5A0jaSvijpKUmrJM2TtG069wFJKyX9vaTngH+VtIukH0laLeml9H5USv+PwP8FrkqfcVU6HpL2Se93knRdyv+kpM9K2iqdO0/Svak8L0n6jaST+vxfyqyfc1A1a7xlwFTgZ6k5dud0/Argj8mmjNuHbGnAS3L53gG8DRgNTCH7//yvaf/dwB+AqwAi4jNkU9JNT58xvUI5vgrsBOwNvB+YDHwsd/5I4HGyaQpnA99O88eaWeKgatYPpWD1/4C/iogX0yLS/0S28ELJBmBmRLweEX+IiBci4gcRsS6l/0ey4Fjk84YAHwEujohXIuIJsiW/zsklezIivhkRb5Et0j0SeHsvb9Wsqfg5jFn/tBuwHfBArjIossUXSlZHxGsbT0rbAV8BTgRKC0vvIGlICoRdGQEMBZ7MHXuSrHZc8lzpTUSsS+UaXvSGzAYD11TN+ofylS3WkDXfHhgRO6dtp4gY3kWevwH2A46MiB2B0mLw6iR9+ee9SdZ0XPJu4LdV3IPZoOegatY/rAJGpbV1yS0Q/RVJuwNI2kPSCV1cYweyQPw7SW8jW5qt/DMqjklNNdn5wD9K2kHSaOCvge/24p7MBh0HVbP+YSGwFHhO0pp07O+B5cB9ktYCd5LVRDvzL8C2ZLXO+4Dbys5fCZyReu/OqZD/k8CrwArgXuDfgWt6djtmg5PXUzUzM6sR11TNzMxqxEHVzMysRhxUzczMasRB1czMrEYcVM3MzGrEQdXMzKxGHFTNzMxqxEHVzMysRv4/yGz2bDweY/MAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from matplotlib import pyplot as plt\n", "fig, axs = plt.subplots(1,1)\n", "axs2 = axs.twinx()\n", "axs.plot(energy, 'kx', label=\"energy\")\n", "axs2.plot(forces, 'r+', label=\"force\")\n", "axs.set_xlabel(\"Iteration\", fontsize=12)\n", "axs.set_ylabel(\"Energy (Hartree)\", fontsize=12)\n", "axs2.set_ylabel(\"Force (A.U.)\", fontsize=12)\n", "axs.legend(loc=\"upper center\")\n", "axs2.legend(loc=\"upper right\")\n", "axs.ticklabel_format(useOffset=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have to take some care with plotting the values from the calculation using implicit solvent. This is because sometimes BigDFT restarts a calculation if it senses that the extrapolated geometry guess is not working good. In that case, one of our logfile instances won't have an energy attribute. A try catch can handle this smoothly." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "energy = []\n", "forces = []\n", "for run in log_imp:\n", " try:\n", " energy.append(run.energy)\n", " forces.append(run.forcemax)\n", " except AttributeError as e:\n", " pass" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdUAAAENCAYAAABZ1P+1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dfbxVZZ3//9c7EMlbCrRICvWbY97kLTf6m74aGKiMo5FJViNaOQhFzbemHK0QxO6kmtJqIDMbNRtDHRvyNu0cS03tQICKmIPkDWkKmqKSN8jn98e6Nqyz2fuctTn7nH3Y5/18PPZjn3Wta611Xfsc+OzrWte6LkUEZmZm1nVvaHQBzMzMmoWDqpmZWZ04qJqZmdWJg6qZmVmdOKiamZnViYOqmZlZnTQ8qEo6SdIySRskjcilj5O0SNJ96X1sbt9Nkpam4+ZJ6lfhvKMkLUmvpZImdna8pNMkrc4dd3p319/MzJqHGv2cqqR9gA3AD4HPR8TClH4w8FREPCFpf+DmiNgt7dspItZKEnA1cFVEXFl23u2AVyNivaShwFLgbWm74vGSTgNGRMT0nqm9mZk1k/6NLkBELAfI4lu79MW5zWXAQEnbRsQrEbE2pfcHBgCbfTOIiHW5zYH5PEWONzMzq1XDg2pBJwKLI+KVUoKkm4FRwI1krc3NSBoNXAIMB06JiPUFjj9R0hHAQ8BnI+Lxzgr3hje8Id74xjfWXCkzs75s3bp1ERENvw1ZTz3S/SvpVuCtFXZ9KSL+J+W5jVz3b+7Y/YAFwPiIeLhs30DgCmBeRNzSwfX3AS4FjoiIl6sdL2kw8GJEvCJpKjApIsZWOecUYArAgAEDDn3llVcqZTMzsyokrYuI7Rtdjnpq+D3VkkpBVdIwoAX4WETcWeW4U4GRnd0HldQKfKFC0K54fBq89GxE7NxZ2bfffvt46aWXOstmZmY5zRhUe22zW9Ig4Hrg7HxAlbRDGniEpP7ABODBCsfvkfYjaTiwN/BIR8eX0pPjgeXdUTczM2tODb+nmh51+R6wC3C9pCURcTQwHXgnMEPSjJR9PCBggaRtgX5kLdl56VzHk43ePQd4D3CWpNfIRhd/MiLWSHpLteOBz6RzrAeeBU7r3tqbmVkz6TXdv1szd/9a3muvvcaqVat4+eWXO8/chwwcOJBhw4axzTbbNLooDdfX/kaq/e6bsfvXQbUOHFQt709/+hM77rgjgwcP3uxRsb4qInjmmWd44YUX2GOPPRpdnIbrS38jHf3umzGo9tp7qs1szpw5tLa2ZhuzZgHQ2trKnDlzGlcoq5uXX365T/xnWQtJDB48uM+0zDrTl/5G+trv3kG1AUaOHMmkSZOywHruubS2tjJp0iRGjhzZ6KJZnfSF/yxr5c+kvb70efSlujqoNsCYMWOYP38+kyZNAmDSpEnMnz+fMWPGNLhkZmbWFQ6qjTBrFmPGjmX1mjUArF6zhjFjx27sCjYz624XXngh++yzDx/96EcbXZSm4qDaCLNm0drSwi5DhgCwy5AhtLa0OKj2Qe3urye99f76+vXrO89k3auO/0f8x3/8BzfccANXXHFFp3n9uy/OQbUBSvdQ58+fD7CxK7j8P1drfu3ur0Nd76//9Kc/ZdSoURx00EGcccYZvP766+ywww586Utf4sADD+Swww7jqaeeAmD16tWceOKJjBw5kpEjR3Lnndl8K7NmzWLKlCmMHz+eyZMns3r1asaNG8chhxzCGWecwfDhw1mzZg0zZszgggsu2HjtL33pS1x44YVdroOVOffcupxm6tSprFy5kuOPP55vf/vbvP/97+eAAw7gsMMO49577wU2/90/9dRTTJw4kQMPPJADDzyQ3/3ud0Dlv7M+LSL86uJru+22i1qcf/750dLSkm3MnBkRES0tLXH++efXdB7rnR544IGa8re0tMSQIUNixowZMWTIkE1/G10sw3HHHRevvvpqRERMmzYtLr300gBiwYIFERHxhS98Ic4777yIiPjwhz8ct99+e0REPProo/Gud70rIiJmzpwZhxxySKSJz+NTn/pUfO1rX4uIiBtvvDGAWL16dfzpT3+Kgw8+OCIiXn/99dhzzz1jzZo1FctlXfgcoG5lGD58eKxevTqmT58es2bNioiIX//613HggQdGxOa/+0mTJsV3vvOdiIhYv359PPfcc1X/ziqpVGfgpegF/4fX89XwGZX6ojPPPHPTRurOGTNmjAcq9VFjxoxh2rRpnHfeecyYMaMufwe//vWvWbRo0cYW79/+9jd23XVXBgwYwHHHHQfAoYceyi23ZOtQ3HrrrTzwwAMbj1+7di0vvPACAMcffzylVZjuuOMOrr32WgCOOeYY3vSmNwGw++67M3jwYBYvXsxTTz3FwQcfzODBg7tcDyP7PyLfQi2NpJ05sy7dwXfccQfXXHMNAGPHjuWZZ57h+eefB9r/7ltaWrjssssA6NevHzvvvDOXX355xb+zepJ0DHAB2Qx4F0fEN8r2K+2fAKwDTouIP0jaG/h5LuuewDkR8d26FrCMg6pZg7W2tjJ37lxmzJjB3Llz6/IFKyI49dRT+frXv94u/Vvf+tbGxxv69eu38V7Zhg0buOuuu6i0hOH22296Nj9rXFR2+umn85//+Z/85S9/4eMf/3iXym85s2ZtCp4SdPA72BKVfqelv5H8777asZX+zuolLWzyA2AcsApok7QgIh7IZTsW2Cu9RgNzgdER8UfgoNx5/gxc2y0FzfE9VbMGyt9fnz17dt3urx911FFcffXVPP300wA8++yzPProo1Xzjx8/nu9///sbt5csWVIx33ve856NYwF+9atf8de//nXjvokTJ3LTTTfR1tbG0Ucf3aXyW8854ogjNg5Wuu222xgyZAg77bTTZvmOOuoo5s6dC8Drr7/O2rVra/472wKjgBURsTIiXgWuBE4oy3MCcFnqUb4bGFS2OArAUcDDEVHXwlXioGrWQG1tbe2eUS49w9zW1tal8+6777585StfYfz48RxwwAGMGzeOJ598smr+Cy+8kIULF3LAAQew7777Mm/evIr5Zs6cya9+9SsOOeQQbrzxRoYOHcqOO+4IwIABAxgzZgyTJk2iX79+XSq/VTFzZt1POWvWrI2/+7POOotLL720Yr4LLriA1tZW3v3ud3PooYeybNmymv/OtsBuwOO57VUprdY8JwP/Vc+CVeO5f+vAc/9a3vLly9lnn30aXYxu8corr9CvXz/69+/PXXfdxbRp0za2ajds2MAhhxzCVVddxV577VXx+Gb+bGrRFz+HSnWW9CpwXy7pooi4KLf/JODoiDg9bZ8CjIqIT+fyXA98PSLuSNu/Bs6MiEVpewDwBLBfRDzVLZXL8T1VMyvsscceY9KkSWzYsIEBAwbwox/9CIAHHniA4447jokTJ1YNqGYVrI+IER3sXwW8Pbc9jCxA1pLnWOAPPRFQwUHVzGqw1157sXjx4s3S9913X1auXNmAElmTawP2krQH2UCjk4GPlOVZAEyXdCXZQKXnIyLfB/1heqjrFxxUzbpFRPSpScSL8K2m9vrS38iW/u4jYr2k6cDNZI/UXBIRyyRNTfvnATeQPU6zguyRmo+Vjpe0HdnI4TO6VIEaOKia1dnAgQN55pln+szSXkVEZGtqDhw4sNFF6RX60t9IV3/3EXEDWeDMp83L/RzAp6ocuw7o0QemPVCpDjxQyfJee+01Vq1a1WfWjyxq4MCBDBs2jG222abRRWm4vvY3Uu1334yLlDuo1oGDqplZ7ZoxqPo5VTMzszppeFCVdJKkZZI2SBqRSx8naZGk+9L72Ny+myQtTcfNS1NQlZ93lKQl6bVU0sQKeRZIuj+3va2kn0taIekeSbvXv8ZmZtasGh5UgfuBDwC/LUtfA/xjRLwbOBW4PLdvUkQcCOwP7AKcVOW8IyLiIOAY4IeSNg7MkvQB4MWyYz4B/DUi3gl8Bzh/i2tlZmZ9TsODakQsTxMfl6cvjojSA7zLgIGStk371qb0/sAAYLMbwxGxLiJKK+sOzOeRtAPwOeArZYedAJTm6LoaOErNPjTPzMzqpuFBtaATgcUR8UopQdLNwNPAC2QBcDOSRktaRjYN1tRckD0P+DbZM015G+eQTHmfp4eHY5uZ2darR4KqpFsl3V/hVb7aQKVj9yPrhm338G5EHA0MBbYFxlY4lIi4JyL2A0YCZ0saKOkg4J0RUWkJoEqt0orDoyVNkbRQ0sLS8llmZta39cjkDxHxvi05TtIwsvXvJkfEwxXO+7KkBWTdtrd0cP3lkl4iuwc7EjhU0iNk9d9V0m0R8V42zSG5Kt1/3Rl4tso5LwIuguyRmi2pn5mZNZde2/0raRBwPXB2RNyZS9+htFZeCnwTgAcrHL9HaWCSpOHA3sAjETE3It4WEbsD7wEeSgEVsjkkT00/fxBoCT/Ia2ZmBTU8qEqaKGkVcDhwfbpXCjAdeCcwI/dozK7A9sACSfcCS8nuq85L5zpe0ux0/HuApZKWkLV2PxkRazopzo+BwZJWkA1kOqt+NTUzs2bnGZXqwDMqmZnVzjMqmZmZWVUOqmZmZnXioGpmZlYnDqpmZmZ14qBqZmZWJw6qZmZmdeKgamZmVicOqmZm1mtJOkbSH9M615tNyKPMhWn/vZIOye0bJOlqSQ9KWi7p8O4ur4OqmZn1SpL6AT8AjgX2BT4sad+ybMcCe6XXFGBubt8FwE0R8S7gQGB5d5fZQdXMzHqrUcCKiFgZEa8CV5ItoJJ3AnBZZO4GBkkaKmkn4Aiy6WeJiFcj4rnuLrCDqpmZ9VYb17hOVqW0Inn2BFYDP5G0WNLFkrp9SkQHVTMza5T+pXWp02tK2f4ia1xXy9MfOASYGxEHAy/RA4uk9Mh6qmZmZhWsj4gRHewvrXFdMgx4omCeAFZFxD0p/Wp6IKi6pWpmZr1VG7BXWh97AHAy2brXeQuAyWkU8GHA8xHxZET8BXhc0t4p31HAA91dYLdUzcysV4qI9ZKmAzcD/YBLImKZpKlp/zzgBmACsAJYB3wsd4pPA1ekgLyybF+38HqqdeD1VM3Mauf1VM3MzKwqB1UzM7M6cVA1MzOrEwdVMzOzOml4UJV0kqRlkjZIGpFLHydpkaT70vvY3L6bJC1Nx81L80OWn3eUpCXptVTSxAp5Fki6P7d9mqTVueNO7446m5lZc+oNj9TcD3wA+GFZ+hrgHyPiCUn7kw2pLk1PNSki1koS2QO9J5HNCVl+3hFpSPZQYKmkX0bEegBJHwBerFCen0fE9LrUzMzM+pSGB9WIWA6Qxcd26Ytzm8uAgZK2jYhXImJtSu8PDGDzaauIiHW5zYH5PJJ2AD5HtqLB/DpUw8zMrPHdvwWdCCyOiFdKCZJuBp4GXiBrrW5G0mhJy4D7gKmlVipwHvBtsgeFN7tWWpPvaklvr7DfzMysoh4JqpJulXR/hVf5Ej6Vjt0POB84I58eEUcDQ4FtgbEVDiUi7omI/YCRwNmSBko6CHhnRFxb4ZBfArtHxAHArcClHZRrSmkS6PXr11fLZmZmfUivmVFJ0m3A5yNiYS5tGNACfCwi7qxy3KnAyM7ug0pqBb5AFmBnAK+SdR/vCvwuIt5blr8f8GxE7NxZ2T2jkplZ7ZpxRqUO76lK6g8cD/wD2arpg4DngKXAjcAvcl2qdSVpEHA9cHY+oKb7oTtGxJOpfBOA2yscvwfweBqoNBzYG3gkBe25Kc/uwHWlgCppaEQ8mU5xPD2wSryZmTWPqt2/ks4gm4D4DOBh4KvA1PT+MPDPwMrSxMZbStJESauAw4Hr071SgOnAO4EZuUdcdgW2BxZIupcsuD8NzEvnOl7S7HT8e8hG/C4BrgU+GRFrOinOZ9JjOkuBzwCndaVuZmbWt1Tt/pX0beCbafmcanmGAv8aEZ/vpvJtFdz9a2ZWu2bs/u0191S3Zg6qZma1a8agWvg5VUn7AB8E3hoRn5L0LmBARNzbbaUzMzPbihR6pEbSScBvyGY0OiUl7wD8ezeVy8zMbKtTqPtX0nLgwxGxRNJfI+JNkrYBnoiIXbq9lL2cu3/NzGrXjN2/RSd/2JVspC1smu4vqDA9oJmZWV9VNKguYlO3b8nJwO/rWxwzM7NNJB0j6Y+SVkg6q8J+Sbow7b9X0iG5fY+klc6WSFpYfmx3KDpQ6TPAryR9Atg+PUv6d8D4biuZmZn1aWlmux8A44BVQJukBRHxQC7bscBe6TWabHKf0bn9YwrMUVA3hYJqRDyYRvseB1wHPE42E1GlpdPMzMzqYRSwIiJWAki6EjgByAfVE4DLIhsgdLekQWWz4/WowhPqp6XU7gJuj4grHVDNzKyb7UbWiCtZxaZ1tYvkCbJe1kWSpnRbKXMKtVQlvQP4L+AgskLuIOmDwDERcXo3ls/MzJpX/7J7nRdFxEW5bZUfwOYDZDvK8/cR8USa4vYWSQ9GxG+7UN5OFb2n+kOyye3/L/BMSruFbE1SMzOzLbE+IkZ0sH8VkF/XehjwRNE8EVF6f1rStWTdyd0aVIt2/44CvhERG0jfACLieaDTZdHMzMy2UBuwl6Q9JA0ge+pkQVmeBcDkNAr4MOD5tIrZ9pJ2BJC0PdnA2vu7u8BFW6pPka0Y81ApQdK+wGPdUSgzM7O0dOd04GagH3BJRCwrrY4WEfOAG8iWAF0BrAM+lg5/C3CtJMhi3c8i4qZK10mTGR3G5kuc3h0Rr9VS5qIzKn0cOAv4OnAB2XJwXyRrvV5RywWbkWdUMjOrXaNnVJI0hCy2nQo8CzwIvADsCOwDvAm4lCzWFXosp+gjNZdIehaYQjbKajIwIyJ+UWslzMzMeonbgR8DB0XEn8t3Snob8FGy+7D7Fjlhpy3V9PDtr4GjI+KVWkvcF7ilamZWu17QUh0QEa/WKx8UGKgUEa8DexTJa2ZmtrUoGiiL5oPigfJcYK6k4ZL6SXpD6VX0QmZmZlsbSf9RS/6io38vTu/5SfVF9nhNv1ouaGZmthWpNLlEVUWD6h5bUBAzM7OtWkRMqyV/0e7bkyLi0fIXcGLtRWxP0kmSlknaIGlELn1cmq/xvvQ+NrfvJklL03Hz0mCq8vOOSsv9LEl5J+b23ZaWEirt3zWlbyvp52kJoXsk7d7V+pmZWe8lac8Kr7dv6e3Nos+pro2InSqkPxsRb96SC+fOsQ+wgWwqxM9HxMKUfjDwVJq3cX/g5ojYLe3bKSLWKnuq92rgqoi4suy82wGvpoeHh5I9yPu2tH1b/lq5Yz4JHBARUyWdDEyMiA91VgeP/jUzq12jR/+mMpRmCizv5l0PXAV8Ms0gWEiH3b+51mE/SWPKLron2UOyXRIRy9O1ytMX5zaXAQMlbRsRr0TE2pTeHxjA5hMsl1bVKRlYKU8FJwCz0s9XA9+XpCjyzcPMzLY6EbFZi1RSf7IY91Xgm2RzNBTS2T3VH6f3gcAl+XIAfwE+XfRCXXQisDj/nGxaKH0UcCNZANyMpNFk5R4OnBIR63O7fyLpdeAa4CspcG5cQii1aJ8HBgM9tsCtmZk1VooVD0k6A7i3lmM7DKoRsQeApCsi4qNbWkBJtwJvrbDrSxHxP50cux9wPtlkyPmyHS1pIHAFMJZs1RzK8twD7Je6mC+VdGNEvAx8NCL+nCZbvoZsVPNlFFtmqFSuKaRvLwMGDOioCmZmtnVaC2xXywGdjv5Ng4A+UOp63ZJSRcT7tuQ4ScOAa4HJEfFwhfO+LGkBWbftZkE1l2+5pJeA/YGFpemoIuIFST8ja/FexqYlhFal5v/OZPNBVjrnRcBFkN1T3ZL6mZlZr/YhstuPhXUaVCPidUkPkXWDlq9j120kDSJbw/XsiLgzl74DsGNa2qc/2eoEt1c4fg/g8dSNOxzYG3gkHTMoItaklQmOA25Nhy0gm1j5LuCDQIvvp5qZNS9Jl7N5j+Q2wO5kcWNCLecr+pzqFcB1ki4ga81tLEBEtNRywXLpUZfvAbsA10taEhFHA9PJlpubIWlGyj6erIt2gaRtySaeaAHmpXMdD4yIiHOA9wBnSXqNbHTxJ1Mg3R64OQXUfmQB9Ufp/D8GLpe0gqyFenJX6mZmZr3eigpp68mWlLspIlbXcrKij9T8qcquiIg9a7lgM/IjNWZmtesNj9TUW9Gl3zyjkpmZ9TmS/j5/C7IznhDfzMysuptqyVyopSppJ7JJEY4EhpB79CQi3lHLBc3MzLYWEbFjLfmLtlT/AzgEmA28mWzSh8eA79RUOjMzs61AWt70HyTNr+W4okF1PHBimqjh9fT+IdovBWdmZlZXko5JC6CskHRWhf2SdGHaf6+kQ8r295O0WNJ1Ba93oKR/B/4M/Ax4upbyFn2k5g1AaULhF9MzpE+SPfJiZmZWd2nyoR8A48ge52yTtCAiHshlOxbYK71GA3PTe8m/AMuBzRaFyV3nLcBHyeYp2Bf4LbAD8O6IeKSWMhdtqS4lu58K2UQLP0gFf6iWi5mZmdVgFLAiIlZGxKvAlWQz6OWdAFwWmbuBQWllstKsfP8AXFztAqkF+zjwEeBS4B0RcRTwIrCu2nHVFA2q/ww8kn7+DPA3YBAwudYLmpmZFbRxkZNkVUormue7wJlkEwBV816yOX5vBG6IiCe7UN7Cz6muzP28Gji9Kxc1MzMD+kvKr2t9UZpXvaTIIicV80g6Dng6IhZJem8HZdiVbFraU4EvSlpKNovgNhWu1ami66lW1dVpCs3MrM9aHxEjOthfWuSkZBibz0FfLc8HgeMlTSBbvnQnST+NiH/KH5zW3r4MuEzSO8h6YKeQPelyuaQLI+KGohXqcJrCCtMTvp32zWxPU4inKTQz2xKdTVOYFkB5CDiKbDRuG/CRiFiWy/MPZHPFTyAboHRhRIwqO897gc9HxHE1lO1w4DTggxExuOhxhdZTzV3kr56y0MzMekJaZWw6cDPZAiiXRMQySVPT/nlkE99PIJsYfx3wsTpd+y7gLkmfqeW4QhPqb8wsPRsRb661cM3OLVUzs9o1ekL99DzqnIj4Swd53gqcGRGfK3LOos+pmpmZNZs/Ar+XtBz4Tdp+AdgR+DuykcF7A18pekK3VOvALVUzs9o1uqWayrAN2bOuxwLvJntc9K/AvWRdy7+MiPWFz9fJQKXbaT+k+HDgrnyeiDii6MWalYOqmVntekNQrbfOun/LZ6H4cXcVxMzMbGtXU/evVeaWqplZ7ZqxpVp1mkJJxxc5QdF8ZmZmza6juX9PlnS/pLMl/X+SBksakN4Pl3SWpPuBST1VWDMzs96salCNiI8AHyabmPhyYDXZRPpPk83k/1bgQ+VTPtVK0kmSlknaIGlELn2cpEWS7kvvY3P7bpK0NB03Ly0PVH7eUZKWpNdSSRNz+25L6/OV9u+a0k+TtDqX7jmOzcz6CEn7SJoh6Qdp+12SDqjpHEXvqUrajmyo8XNprsS6kLQP2QoCPySbRmphSj8YeCoinpC0P3BzROyW9u0UEWslCbgauCoirqxQ3lfTjBxDyZave1vavi1/rdwxpwEjImJ6LXXwPVUzs9r1pnuqkk4iW9b0v8mmQtwpNfS+ERHvK3qewpM/pEBat2CaO+9ygCw+tktfnNtcBgyUtG1EvBIRa1N6f2AAFVYSKAv8AyvlMTMzS2YD4yNiiaQPpbSlwIG1nKToeqqNdiKwOCJeKSVIupmsK/oFstbqZiSNlrQMuA+YWvYA709SF+8MtY/oJ0q6V9LVkt6OmZn1BbuSBVHY1AgLamyQ9UhQlXRrGvRU/ipfwb3SsfsB5wNn5NMj4mhgKLAtUHGJuoi4JyL2A0YCZ0samHZ9NCLeDfzf9Dolpf8S2D0iDgBuJbt3XK1cUyQtlLRw/frCk22YmVnvtIhNsaDkZOD3tZykR4JqRLwvIvav8Pqfjo6TNAy4FpgcEQ9XOO/LwAKyKaY6uv5y4CVg/7T95/T+AvAzYFTafibXGv4RcGgH57woIkZExIj+/eszhfKcOXNobW1tl9ba2sqcOXPqcn4zM6vqM8BXJP0G2D71hp4HfLaWkxQKqpI+I2lI7WXccpIGAdcDZ0fEnbn0HdLAo9JaexOAByscv0faj6ThZJMiPyKpf6kuac7H44D70/bQ3CmOB5Z3R92qGTlyJJMmTcoC66xZtLa2MmnSJEaOHNmTxTAz63Mi4kHgXWSDlb4M/AR4d0T8by3nKTT6V9ICsi7W28ger/lF/v5mV6RHXb4H7AI8ByyJiKMlfRk4G8hXaDwg4Dqybt9+QAvw2TSq93iy0bvnSDoFOAt4jWx08eyI+IWk7YHfAtuk428FPhcRr0v6OlkwXQ88C0xLH3SH6jn6txRIV69Zwy5DhjB//nzGjBlTl3ObmfUmvWz0727Auoj4ay7tTcAbI+KJwuep4ZGawWT9y/9EFs2vAS6LiN/WUvBmVO9Has455xxmn3ce58yYwezZs+t2XjOz3qSXBdU24OMRcV8u7d3AxRExuvB5tmTu3/Qw7OVk9ygfJ7v/eEFEvFjzyZpA3YLqrFlw7rmbp8+cme0zM2sivSyoPh8ROxdNr6amgUqSjpL0E7Ju4KeAyWSjpQ4GbqzlXLa51iOPZJchQ2htacm2W1qy7SOPbHDJzMya3mpJ78wnpO1najlJ0YFK35K0CriQbFDQuyNifERcERG3k01neHAtF7bNtbW1tbuHOmbMGObPn09bW1uDS2Zm1vQuAa6RdJykfSX9I9kcCOVLoHao6ECl7wOXRkTV/90lvavIoJ5m1C3TFM6a5S5fM2tqRbp/JR0DXEA2sPTiiPhG2X6l/RPIZv07LSL+kOYl+C3ZoNb+wNURMbOD67wB+FfgE8DbgcfI1hD/94jYULhOXk+16zz3r5lZ7ToLqmmxlIeAccAqoA34cEQ8kMszAfg0WVAdTTa+Z3QKtttHxIvp8ck7gH+JiLurXOcSYEpXn2wpNGuBpNupPFXTK2QV/e+I+GVXCmJmZlZmFLAiIlYCSLqSbLKfB3J5TiB7EiWAuyUNkjQ0Ip4ESoNnt0mviq3I9EjleLLHL7uk6ECl24Ddgd8AP03vw4GFZAOWLpF0ZlcLY2ZmlrMb2RMmJatSWqE8kvpJWkI2T02U+UAAABorSURBVPwtEXFPB9f6DnCupAFdKXDR+fXGA0eXVpQBkHQF2X3W0ZL+G7gS8Hx6ZmZWVH9J+SU4L4qIi3LbKj+AzVubVfNExOvAQWmGvmsl7R8R91cpy6fJ1gn/nKTV+etExDs6qcdGRYPqu4CVZWmPkk39R0T8vrTQt5mZWUHrI2JEB/tXkQ0aKhkGlM9u1GmeiHguraN9DGla2gr+qUiBO1M0qP6WbKm0c8gqMAyYRXbjtzTrxJP1KJCZmVnSBuwlaQ/gz2Sz+n2kLM8CYHq63zoaeD4inpS0C/BaCqhvBN5HtuJZRRHxm3oUuOg91VNT3gfIVntZRja8+bS0/1WyZ1XNzMzqIq2BPR24mWyBk/kRsUzSVElTU7YbyHpSV5DN7vfJlD4UaJV0L1lwviUirqt2LUnbSDpX0kpJL6f3mu+xdvpITRpqPBP4Glnw3AVYXctzO83Oj9SYmdWul01T+B2y0cbnkt3eHA7MABZGROHl34pO/vAMsIsDaWUOqmZmtetlQXUVcGBEPJNLGwIsjYjyEcdVFe3+vRSY2mkuMzOzrVOlUcQdpVdUdKDSKODT6VnUx2k/1PiIWi5oZmbWC10F/FLSuWRTFA4nW6x8fi0nKdr9e2q1fRFxaS0XbEbu/jUzq10v6/4dQBZEPwK8jWy08ZXAV2qZutBz/9aBg6qZWe16Q1CV9NaI+Eu9zld06TdJ+mdJLWl4MpKOkDSpXgUxMzNrgIfyG2mGwC1WdKDSbLLlcC4CStM1rQL+rSsXNzMza7DygUjv7crJigbV04DjIuJKNg1S+hOwZ1cubmZm1mB1vQdaNKj2Y9MSOqUC7JBL22KSTpK0TNIGSSNy6eMkLZJ0X3ofm9t3k6Sl6bh5aYKK8vOOkrQkvZZKmpjbN0DSRZIekvSgpBNT+raSfi5phaR7JO3e1fqZmVmv1l/SGEljU5xpt52PPUUUHf17MdlsSp8lm+N3MNkyOQMi4pMdHVvg3PuQrWH3Q+DzEbEwpR8MPBURT0jaH7i59ACupJ0iYm1ahPZq4KrUis6fdzvg1YhYL2kosBR4W9o+F+gXEV9Oq72/OSLWSPokcEBETJV0MjAxIj7UWR08UMnMrHa9ZKDSI3TcWo2IKNwrW/Q51c8BlwHPky30+iLwK2By0QtVU1pOLouP7dIX5zaXAQMlbRsRr0TE2pTeHxhAhQ8kItblNgeW5fk42co7pFmi1qT0E8gWCoAsWH9fksJDpM3MmlJE7F7P8xXq/o2ItRHxfrJBSocB/yciJkbEC/UsTAdOBBbnnxWSdDPZwrMvkAXAzUgaLWkZcB8wNbVSB6Xd50n6g6SrJL0lpW1c7DZN5Pw8WavczMysU0XvqeY9A2wnaU9JhZrEkm6VdH+F1wkFjt2PbLmeM/LpEXE02SoE2wIV+7wj4p6I2A8YCZwtaSBZ63YYcGdEHALcBXyrdLlKp6lSrimSFkpauH79+s6qYWZmfUCh7l9JxwA/JgtieUE2iKlDEfG+2osGkoYB1wKTI+LhCud9WdICsm7bWzq4/nJJLwH7A4uAdem8kE1N9Yn0c2mx21WS+gM7A89WOedFZI8Ysf3227t72MzMCrdUfwCcB2wfEW/IvToNqFsqddNeD5wdEXfm0ndIA49IgW8C8GCF4/dI+5E0HNgbeCTdH/0lm55FOopsnVjIFrstTcn4QaDF91PNzKyooqN/nwUGd0eASY+6fI9sndbngCURcbSkLwNnA/+byz6erIv2OrJu335AC/DZdL/0eGBERJwj6RTgLOA1stHFsyPiF+maw4HLgUHAauBjEfFY6h6+HDiYrIV6ckSs7KwOHv1rZla73jD6t96KBtVvAssj4pLuL9LWx0HVzKx2fTmo3k62/NujQLuJh730m4OqmdmWaMagWvQ51YvTy8zMrMekgbIXkN3uuzgivlG2X2n/BLJBqKdFxB8kvZ1sfoW3kt0CvCgiLuju8hYKql4z1czMelqagvYHwDiypzPaJC2IiAdy2Y4F9kqv0cDc9L4e+NcUYHcEFkm6pezYuutw9K+kC8u2P1G2fU13FMrMzIzstuOKiFgZEa+SLRpePr/BCcBlkbkbGCRpaEQ8GRF/AEgTFS0nm+CnW3X2SM1pZdvfLNseV7+iWDVz5syhtbW1XVpraytz5sxpUInMzHrExlnuklVsHhg7zZMWRzkYuKfuJSzTWVAtn2Go0oxD1s1GjhzJpEmTNgbW1tZWJk2axMiRIxtcMjOzLulfmpkuvaaU7S8yy12HeSTtAFwD/L/cvPHdprN7quWF90QIDTBmzBjmz5/PpEmTmDZtGnPnzmX+/PmMGTOm0UUzM+uK9RExooP9pVnuSoYBTxTNI2kbsoB6RUT8d9eL27nOgmp/SWPY9E2gfLvbZlSy9saMGcO0adM477zzmDFjhgOqmfUFbcBekvYA/gycDHykLM8CYLqkK8kGKD0fEU+mUcE/Jptj4d97qsAdPqdaYJ05ImKPOpdpq9MTz6mWunzdUjWzZlHkOVVJE4DvkjXiLomIr0qaChAR81Lw/D5wDNkjNR+LiIWS3gPcTrZK2YZ0ui9GxA3dVJ2svJ7atuu6O6iWAmopkJZvm5ltjZpx8octWfrNelhbW1u7AFq6x9rW1tbgkpmZWZ5bqnXgaQrNzGrnlqo13qxZjS6BmZlV4ZZqHfRoS1UC/87MrAm4pWpmZmZVOahuDWbNylqoSo8Hl352V7CZWa/i7t86cPevmVnt3P1rZmZmVTmobm1mzmx0CczMrAoH1a1N2X1ULwtnZtZ7OKhu5bwsnJlZ79HwoCrpJEnLJG2QNCKXPk7SIkn3pfexuX03SVqajpsnabPVciSNkrQkvZZKmpjbN0DSRZIekvSgpBNT+mmSVueOO727699V+WXhzjnnHM8JbGbWQJ0t/dYT7gc+APywLH0N8I8R8YSk/YGb2bSa+6SIWJtWJ7gaOAm4ssJ5R0TEeklDgaWSfhkR64EvAU9HxN9JegPw5txxP4+I6XWtYTfzsnBmZr1Dw1uqEbE8Iv5YIX1xRJQWo10GDJS0bdpXWr29PzCACsvTRcS6FEABBpbl+Tjw9ZRvQ0SsqUtlGqS1tZW5c+cyY8YM5s6du9k9VjMz6xkND6oFnQgsjohXSgmSbgaeBl4ga61uRtJoScvI1tObmlqtg9Lu8yT9QdJVkt6Sv5akeyVdLentm59147mnSFooaeH69eurZet2+WXgZs+evbEr2IHVzKzn9UhQlXSrpPsrvE4ocOx+wPnAGfn0iDgaGApsC4ytcCgRcU9E7AeMBM6WNJCsdTsMuDMiDgHuAr6VDvklsHtEHADcClxarVwRcVFEjIiIEf37N64X3cvCmZn1Hr1mRiVJtwGfj4iFubRhQAvZSu53VjnuVGBkZ/dBJbUCXwAWAS8CO0bEhtQavSkF33z+fsCzEbFzZ2X30m9mZrXzjEo9KHXTXg+cnQ+oknZIA4+Q1B+YADxY4fg90n4kDQf2Bh6J7FvEL4H3pqxHAQ+kfENzpzgeWF7napmZWRNreEs1PeryPWAX4DlgSUQcLenLwNnA/+ayjwcEXEfW7duPrCX72XS/9HiyEb/nSDoFOAt4DdgAzI6IX6RrDgcuBwYBq8lawo9J+jpZMF0PPAtMi4jNAnY5t1TNzGrXjC3VhgfVZtBrguqsWV65xsy2GkWCqqRjgAvIGlEXR8Q3yvYr7Z8ArANOi4g/pH2XAMeRPUK5fzdUYfPyOqh2Xa8Jql7Bxsy2Ip0F1TS25SFgHLAKaAM+HBEP5PJMAD5NFlRHAxdExOi07wiyMTSX9VRQ7bX3VM3MrM8bBayIiJUR8SrZJD/lT42cQBY0IyLuBgaVxsdExG/JbuX1GAfVrV2VBczvHDfOE+2b2dZuN+Dx3PYqNs2sV0ueHuOgurWbNSvr8i11+6afX/3iFz3Rvpn1dv1Lk+ik15Sy/apwTPk9riJ5ekxvmPvXukF+ov1p06Yxd+5cT7RvZr3N+ogY0cH+VUB+ZrthwBNbkKfHuKXaTMoWMM9PtD9t2jQHVDPb2rQBe6V5BwYAJwMLyvIsACYrcxjwfEQ82dMFLXFQbSZlj9N4on0z25qlRVGmk61SthyYHxHLJE2VNDVluwFYCawAfgR8snS8pP8im4p2b0mrJH2iu8vsR2rqoNc8UpOTn2h/zG9+Q+uRR3qtVTPrVZpx8gffU21S7SbaHzuWMRHMnz+fb37zmwDtAmtrayttbW2ceeaZjSqumVlTcEu1DnpjS7Wd3KQQ7VqwY8Zstm1m1lOasaXqe6rNqsrzq2N+85tNa666S9jMrK7cUq2DramlWnLOOecw+7zzOGfGDGbPnt2ggplZX+aWqjWF0qhgwKOCzczqyEG1L8g9v/qnU09lzNixrF6zBoDVa9YwZuxYfnnooZ7W0MysixxU+4Lc86tX7bcfrS0t7aY1bG1p4boRIzZNazhrlqc1NDPbAr6nWge9/p5qNWX3WkuBdPWaNewyZIgHMJlZt/I9VWsuVaY1BDZOazhnzhx3C5uZFeSWah1stS3VvFmz4NxzN0v+0+TJjLrhBs/MZGZ155aqNa3WI49klyFDsvutQGtLC7sMGcIjp5228blWzj13Y0Bta2vb1IJN92w3a8GWzUVsZtbsHFQNKJvWkE1Lx7W1tVXsFh45cuSmgU3nnlt5YFOFli/gYGtmzSsiGvoCTgKWARuAEbn0ccAi4L70Pja37yZgaTpuHtCvwnlHAUvSaykwMaXvmEtfAqwBvpv2bQv8nGy1g3uA3YvUYbvttoumMnNm+583LYO+6TVzZrS0tMSQIUMiIIYMGRItLS1x/vnnR0tLS3YsRERES0tLHHvssRXTzz///J6rl5n1KsBL0eAYVO9X4wsA+wB7A7eVBdWDgbeln/cH/pzbt1N6F3ANcHKF824H9E8/DwWeLm2X5VsEHJF+/iQwL/18MvDzInVouqCaUwqcLS0tEbBxe+XkyRWD7bMHHFAx/Xfjx1c8z8ZAG9E+mBdRa/4y7b4A5OrrQG/WM5oxqDa8+zcilkfEHyukL46I0urty4CBkrZN+9am9P7AAGCz0VYRsS6ytfgABlbKI2kvYFfg9pR0AnBp+vlq4CipNHlu31StW7j0vOsuQ4YAbLwf+6alSyumH3744RsnmgA2TkAx4Gtf23RvNnUXdzS6uN1o5AL5O9KuC5tNjxQ9/PDDFUc8T5gwoaaR0LWOnG7kSOtq165W53p9Ft19/np+dr2tDh6Z30s1OqqXXpS1VMv2fRC4tSztZuCvwM+o0P2b8owmC8gvkrp/y/afA3wrt30/MCy3/TAwpMq5pwALgYUDBgyIPiHXMqzWgi21/GbMmBEB2XtOeXq180yZMmVTKzJdt6WlJaZMmVL1uu1anrlj2nU959JL+YcMGRItRxyx8TztypTr5v72t79dMb1QWbshf7V6dVTfap9RtWtXq3O19Go9D7V+pu3Sc8e3+yxy525X/mr5C34W1T6/ateoVtZuqUOB/OV5Swr1zJT1/hT6jKqdqwCasKXaUwHz1hSwyl8n5PJUDKrAfim4/Z8K+waSdf+O6+T6+wC/BwaWpT8AHJrbXlYhqA7urH7N3P1bTbV/bNUCVWl/+T3YavdsV06eXFO382b/seeOafefW8EvABXLWiW9oy8Y3Zm/Wr06qm+t1671s6h2Tz3/d1Hr+Yt84Wm3r1r+gp9Foc+vyt92kfQu1aFA/lq/OLUL2lV+Zx1+RtW+VBXgoNqdBakQVIFhwEPA33dw3KnA9wucv5X292wPBB4qy3MzcHj6uT/ZICZ1du6+GFSrKf/HVfQfYS2BrVr+jo6pNb2ja9SjrPXKvyX1qledK6V3FMTqcf56f9bd/XfRE3Wo5XdQ65e8Lf2MinJQ7c6CbD5QaRDZqN0Ty/LtAAxNP/cnG607vcL59mDTQKXhwBP5rlzgG8C5Zcd8ivYDleYXKbuD6ibVupg66i6q+I+zxlHHEVH9mCOPrJhebVBVtfRq59lq0mfOrPkzamRZ6xJ4Ovg7qvmz6KBM9QrO3f1lq+J5tuQzSi3gal8kinJQ7Y4CwERgFfAK8BRwc0r/MvAS7R9/2RV4C9AG3Ju6a7+XC57HA7PTz6ek/UuAPwDvL7vuSuBdZWkDgavIHqn5PbBnkTo4qG65Il1b+f8YOsqfP1+R/2SqdVVWu29bpFuwalm7IX+t9S3/zItce0u6mCOKt56qned348dv/h85tX8RWjl5ct1aqtVuPTy25541pddah3qld/hlq8LvrMPPqJNgW1SRoAocA/wx/b98VoX9Ai5M++8FDil6bHe8eiRwNvvLQXXLVWvZVgts1QZtFL7/UyEARESUgupmZcoF22rP2hYqazfkr+c91WrXrlbnjp47ruULTK3nr9dntyVfJOr1GXX330u130FXvrRt0b+pAjoLqkC/NLZlT7InPZYC+5blmQDcmILrYcA9RY/tjlfDA1IzvBxU66+jgVC1HtPRaNiNqn3DLpBeqKzdkL+eo3/r9Vm0+8+1wKjgQucv8IWnXfmr5e/i6N9q16hW1m6pQ4H8tX5xqjaAabMvQ9U+o2plKqBAUD281HuZts8Gzi7L80Pgw7ntP5LNTdDpsd3xanhAaoaXg6pZplsm1OhCkO829SpTN6R36UtbUtPvrAufd4Gg+kHg4tz2KZQNTAWuA96T2/41MKLIsd3x8io1ddAUq9SYmfUwSa+STUVbclFEXJTbfxJwdEScnrZPAUZFxKdzea4Hvh4Rd6TtXwNnknX7dnhsd+jfnSc3MzPrwPqIGNHB/lXA23Pbw8ie5CiSZ0CBY+uu4dMUmpmZVdEG7CVpD0kDyB51XFCWZwEwWZnDgOcj4smCx9adW6pmZtYrRcR6SdPJJubpB1wSEcskTU375wE3kI0AXgGsAz7W0bHdXWbfU60D31M1M6udpHURsX2jy1FP7v41MzOrE7dU60DSBuBvW3h4f2B9p7mai+vcN7jOfUNX6vzGiGiqxp2DaoNJWtjJ6Lem4zr3Da5z39AX69yRpvqGYGZm1kgOqmZmZnXioNp4F3Wepem4zn2D69w39MU6V+V7qmZmZnXilqqZmVmdOKg2kKRjJP1R0gpJZzW6PN1B0iWSnpZ0fy7tzZJukfS/6f1NjSxjPUl6u6RWScslLZP0Lym9mes8UNLvJS1NdT43pTdtnUsk9ZO0WNJ1abup6yzpEUn3SVoiaWFKa+o618pBtUEk9QN+ABwL7At8WNK+jS1Vt/hP4JiytLOAX0fEXmTLNDXTF4r1wL9GxD5kCyZ/Kv1em7nOrwBjI+JA4CDgmDQHazPXueRfgOW57b5Q5zERcVDuMZq+UOfCHFQbZxSwIiJWRsSrwJXACQ0uU91FxG+BZ8uSTwAuTT9fCry/RwvVjSLiyYj4Q/r5BbL/cHejuescEfFi2twmvYImrjOApGHAPwAX55Kbus5V9MU6V+Wg2ji7AY/ntleltL7gLWkVCdL7rg0uT7eQtDtwMHAPTV7n1A26BHgauCUimr7OwHfJ1u3ckEtr9joH8CtJiyRNSWnNXueaeJWaxlGFNA/FbhKSdgCuAf5fRKyVKv26m0dEvA4cJGkQcK2k/Rtdpu4k6Tjg6YhYJOm9jS5PD/r7iHhC0q7ALZIebHSBehu3VBunyOK7zeopSUMB0vvTDS5PXUnahiygXhER/52Sm7rOJRHxHHAb2X30Zq7z3wPHS3qE7NbNWEk/pbnrTEQ8kd6fBq4lu43V1HWulYNq4zRkAd1eYgFwavr5VOB/GliWulLWJP0xsDwi/j23q5nrvEtqoSLpjcD7gAdp4jpHxNkRMSwidif7t9sSEf9EE9dZ0vaSdiz9DIwH7qeJ67wlPPlDA0maQHZfprSA7lcbXKS6k/RfwHuBIcBTwEzgF8B84B3AY8BJEVE+mGmrJOk9wO3AfWy61/ZFsvuqzVrnA8gGqPQj+6I+PyJmSxpMk9Y5L3X/fj4ijmvmOkvak6x1Ctmtw59FxFebuc5bwkHVzMysTtz9a2ZmVicOqmZmZnXioGpmZlYnDqpmZmZ14qBqZmZWJw6qZn2YpBfToxJmVgcOqmYNlJbSep+k0yTd0c3Xuk3S6fm0iNghIlZ253XN+hIHVbMmIMnzeJv1Ag6qZo23DzAPODx1xz4HIGlbSd+S9JikpyTNS9MAIum9klZJ+jdJfwF+IulNkq6TtFrSX9PPw1L+rwL/F/h+usb3U3pIemf6eWdJl6XjH5X0ZUlvSPtOk3RHKs9fJf1J0rE9/kmZ9XIOqmaNtxyYCtyVumMHpfTzgb8jW/j7nWRLA56TO+6twJuB4cAUsn/PP0nb7wD+BnwfICK+RDZ94vR0jekVyvE9YGdgT+BIYDLwsdz+0cAfyaacnAP8WM2+/I5ZjRxUzXqhFKz+GfhsRDybFjz/Gtnk7SUbgJkR8UpE/C0inomIayJiXcr/VbLgWOR6/YAPAWdHxAsR8QjwbeCUXLZHI+JHaZm3S4GhwFu6WFWzpuL7MGa90y7AdsCiXGNQZJPWl6yOiJc37pS2A75Dtuzam1LyjpL6pUDYkSHAAODRXNqjZK3jkr+UfoiIdalcOxStkFlf4JaqWe9QvrLFGrLu2/0iYlB67RwRO3RwzL8CewOjI2In4IiUrir5y6/3GlnXcck7gD/XUAezPs9B1ax3eAoYltbWJSI2AD8CviNpVwBJu0k6uoNz7EgWiJ+T9GayZfbKr1HxmdTUkp0PfFXSjpKGA58DftqFOpn1OQ6qZr1DC7AM+IukNSnt34AVwN2S1gK3krVEq/ku8EayVufdwE1l+y8APphG715Y4fhPAy8BK4E7gJ8Bl2xZdcz6Jq+namZmViduqZqZmdWJg6qZmVmdOKiamZnViYOqmZlZnTiompmZ1YmDqpmZWZ04qJqZmdWJg6qZmVmdOKiamZnVyf8PZMyEwWXzG2kAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from matplotlib import pyplot as plt\n", "fig, axs = plt.subplots(1,1)\n", "axs2 = axs.twinx()\n", "axs.plot(energy, 'kx', label=\"energy\")\n", "axs2.plot(forces, 'r+', label=\"force\")\n", "axs.set_xlabel(\"Iteration\", fontsize=12)\n", "axs.set_ylabel(\"Energy (Hartree)\", fontsize=12)\n", "axs2.set_ylabel(\"Force (A.U.)\", fontsize=12)\n", "axs.legend(loc=\"upper center\")\n", "axs2.legend(loc=\"upper right\")\n", "axs.ticklabel_format(useOffset=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next we will extract the geometry." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "systems = []\n", "for step in log:\n", " systems.append(deepcopy(sys))\n", " systems[-1].update_positions_from_dict(step.astruct)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "systems_imp = []\n", "for step in log_imp:\n", " systems_imp.append(deepcopy(sys))\n", " systems_imp[-1].update_positions_from_dict(step.astruct)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can create a picture now with all three computed geometries overlapping." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "application/3dmoljs_load.v0": "
\n

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n jupyter labextension install jupyterlab_3dmol

\n
\n", "text/html": [ "
\n", "

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n", " jupyter labextension install jupyterlab_3dmol

\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "viz = InlineVisualizer(400, 300)\n", "viz.display_system(systems[0], colordict={x: \"black\" for x in systems[-1]}, show=False)\n", "viz.display_system(systems[-1], colordict={x: \"red\" for x in systems[-1]}, show=False)\n", "viz.display_system(systems_imp[-1], colordict={x: \"blue\" for x in systems[-1]}, show=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Slab Example\n", "In this section, we will study a more sophisticated example using a slab of NaCl as our test case. Here will show how to do a constrained optimization by fixing the positions or our slab, and allowing a molecule to move freely on top of it. First, will create a slab of NaCl built with the helper of the Atomic Simulation Environment." ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "from ase.build import bulk\n", "\n", "atoms = bulk('NaCl', 'rocksalt', a=5.64, orthorhombic=True)\n", "atoms *= [4, 4, 4]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we convert to the BigDFT system format. Note that we set the middle cell value to infinity so that we can explore this system as a surface (instead of the bulk)." ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "from BigDFT.Interop.ASEInterop import ase_to_bigdft\n", "from BigDFT.Systems import System\n", "from BigDFT.UnitCells import UnitCell\n", "sys = System()\n", "sys[\"SUR:1\"] = ase_to_bigdft(atoms)\n", "sys.cell = UnitCell([float(atoms.cell[0, 0]), float(\"inf\"), float(atoms.cell[2, 2])], units=\"angstroem\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We won't be optimizing the lattice constant, instead we want to optimize something that is sticking to the surface." ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "from BigDFT.IO import XYZReader\n", "from BigDFT.Fragments import Fragment\n", "with XYZReader(\"O2\") as ifile:\n", " sys[\"ABS:2\"] = Fragment(xyzfile=ifile)\n", "sys[\"ABS:2\"].translate([x - y for x, y in zip(sys[\"SUR:1\"].centroid, sys[\"ABS:2\"].centroid)])\n", "sys[\"ABS:2\"].translate([0, 0.4*sys.cell[2, 2], 0])" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "application/3dmoljs_load.v0": "
\n

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n jupyter labextension install jupyterlab_3dmol

\n
\n", "text/html": [ "
\n", "

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n", " jupyter labextension install jupyterlab_3dmol

\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "viz = InlineVisualizer(400, 300)\n", "viz.display_system(sys)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Also, we will try fixing the positions of our NaCl atoms and just allow the optimizer to operate on the oxygen molecule." ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "sys[\"SUR:1\"].frozen = \"fxyz\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Build our input file." ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "inp = Inputfile()\n", "inp.set_xc(\"LDA\")\n", "inp.set_hgrid(\"0.37\")\n", "inp.optimize_geometry(method=\"SQNM\", betax=1.0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And now we're ready to run." ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Found 51 different runs\n" ] } ], "source": [ "log = code.run(input=inp, posinp=sys.get_posinp(), name=\"geom-opt-slab\", run_dir=\"work\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's look at how the geometry optimization procedure went. We know from the above print statement \"Found 51 different runs\" that the optimization did not converge, and instead just ran the maximum number of states. For a real study, we would need to continue the run, or perhaps tweak our optimization parameters.\n", "\n", "Let's look at how the energy changes at each step. We can also monitor the distance between the O2 molecule and the surface." ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [], "source": [ "from BigDFT.Fragments import pairwise_distance\n", "\n", "systems = []\n", "energy = []\n", "distance = []\n", "for step in log:\n", " systems.append(deepcopy(sys))\n", " systems[-1].update_positions_from_dict(step.astruct)\n", " energy.append(step.energy)\n", " distance.append(pairwise_distance(systems[-1][\"SUR:1\"], systems[-1][\"ABS:2\"]))" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdUAAAEMCAYAAACSiCwQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de5zVVb3/8debGQlQFAQpT6joEQ9SCSpjlidjNG/pT/PSJKdS6yg6WnYyIz02gHDSGCutNNTMS3m8kJVppzJhBkPNGki8QXlBVDIVFA0Rjcvn98d37eE7e/be890ze8++fZ6Px3rM3ut7W9+57M+s9V0XmRnOOeec67sBpS6Ac845Vy08qDrnnHMF4kHVOeecKxAPqs4551yBeFB1zjnnCsSDqnPOOVcgHlSdc85VJEl1kh6W9Kss2ydLWirpCUn39UeZ6vvjIs4551wRfAlYDmyfvkHSMOAHwJFm9rykUf1RIA+qRTRgwAAbPHhwqYvhnHMV5a233jIzy9mSKmk0cDTwDeC8DLv8B/BzM3sewMxeKXhBM/CgWkSDBw9m/fr1pS6Gc85VFEkbEux2BTANGJpl+17ANpIWhn2+a2Y/LkwJs/Og6pxzrtzUS1oce3+tmV2beiPpGOAVM1siaXK2cwD7A4cCg4E/SHrIzJ4sVqFTF3XOOefKySYzm5Rj+0HAsZI+DgwCtpd0s5l9JrbPKmCNma0H1kv6PTABKGpQ9d6/zjnnKoqZXWhmo81sDHAy0JYWUAF+CXxEUr2kIcAHiTo1FZXXVJ3rJxs3bmTVqlW8/fbbpS5KWRk0aBCjR49mm222KXVResV/rr1X6J+9pLMAzOxqM1su6bfAo8AW4Doze7wgF8pVBl/6rXi23XZb845KLuXZZ59l6NChjBgxAkmlLk5ZMDNeffVV1q1bx+67717q4vSK/1x7J9fPXtJbZrZtiYrWJ978W2ZaW1tpb2+P3sycCUB7ezutra2lK5QriLfffts/eNNIYsSIERVdy/Ofa+9Uw88+k5IHVUmfDLNdbJE0KZY/UNINkh6T9Ei8h5ekT0l6NBzXmna+JknLwrZbMlxvaJhhI5XWSLoibZ+TJFlaeU6V9FRIpxb0mxDT0NBAU1NTFFgvvpj29naamppoaGgo1iVdP/IP3u6q4XtSDfdQCtX4fSuHZ6qPAycA16TlnwFgZh8IM2H8RlIDMBy4DNjfzFZLuknSoWa2QNJY4ELgIDNbm2kGDTNbB0xMvZe0BPh57P1Q4Fzgj7G8HYEZwCTAgCWS7jKztQW4/y4aGxuZN28eTU1NrAaampqYN28ejY2Nhb6Uc865Ait5TdXMlpvZXzNsGg8sCPu8ArxOFNT2AJ40s9Vhv/nAieH1GcBVqWDX0wwaIQiPAhbFsmcDrUC8TeII4F4zey2c+17gyMQ3mY+ZM2k85BBWr1kDwOo1a2g85JDOpmDnnEtXV1fHxIkTed/73seECRP4zne+w5YtWwBYvHgx5557btZjV65cyS23dGvUc71U8qCawyPAcaE79O5Eg3h3AZ4GxkkaI6ke+ETIh2gGjb0kPSDpIUk9Bb4pwO0WemtJ2hfYxczSJ2d+L/BC7P2qkFd4M2fS3tbGTiNHArDTyJG0t7V1DaoeYKtel2frQbk+W9+0aVOpi1CZCvh3PHjwYJYuXcoTTzzBvffey69//WsuvvhiACZNmsT3vve9rMd6UC2sfgmqkuZLejxDOi7HYdcTBa/FRNNRPUg0IHgt0AzcTlTDXAmk/qrrgbHAZKKAeV2YVDmbk4FbQxkHAJcDX8l0CxnyMnabljRV0mJJi3vzYZN6hjpv3jyAzqbgLh+w4Y+lGw+2VaPLs3Uo6LP1m2++mQMOOICJEydy5plnsnnzZrbbbjsuuugiJkyYwIEHHsjLL78MwOrVqznxxBNpaGigoaGBBx54AICZM2cydepUDj/8cE455RRWr17NYYcdxn777ceZZ57Jbrvtxpo1a2hpaeG73/1u57UvuuiinB/wNSXb33EfjRo1imuvvZYrr7wSM2PhwoUcc8wxANx3331MnDiRiRMnsu+++7Ju3TouuOACFi1axMSJE7n88stZuXIlH/nIR9hvv/3Yb7/9ePDBBwFYuHAhkydP5qSTTmLcuHF8+tOfJjV6pKOjgw9/+MNMmDCBAw44gHXr1rF582a++tWv0tDQwD777MM116Q/4atSZlYWCVgITMqx/UFgfIb8qUBreH01cFps2wKgIcv5JhA1I6fe7wCsIQrSK4maf18kanKeAlwT2/caYEpP9zRkyBDL15w5c6ytrS16M2OGmZm1tbXZnDlztu4EmQ/Olu/KwrJly/Lav62tzUaOHGktLS02cuTIrb8XfSzDMcccY//85z/NzKy5udluuukmA+yuu+4yM7OvfvWrNnv2bDMzmzJlii1atMjMzJ577jkbN26cmZnNmDHD9ttvPwsTn9s555xjl1xyiZmZ/eY3vzHAVq9ebc8++6ztu+++Zma2efNm22OPPWzNmjUZy1Wpel32Av69brvttt3yhg0bZi+99JK1t7fb0UcfbWZmxxxzjN1///1mZrZu3TrbuHFjl+1mZuvXr7cNGzaYmdmTTz5p+++/v5mZtbe32/bbb28vvPCCbd682Q488EBbtGiRvfPOO7b77rvbn/70JzMze+ONN2zjxo12zTXXdP4evf3227b//vvbihUrupUz0/cPWG9lEJd6k8qho1JGYQYMmdl6SYcR1VKXhW2jzOwVScOBs4GmcNidRAHwRkkjiZqDV2S5xBRCLRXAzN4ARsauvxA438wWS1oBXBKuB3A4UYeogps2bdrWN6Hm2djYSON990G8p1zq9YwZXkOtUo2NjTQ3NzN79mxaWloK0lltwYIFLFmypLPGu2HDBkaNGsXAgQM7azP7778/9957LwDz589n2bJlncf/4x//YN26dQAce+yxpFZhuv/++/nFL34BwJFHHsnw4dGfypgxYxgxYgQPP/wwL7/8Mvvuuy8jRozo831UrJkzu9ZQi/h3bBnmIDjooIM477zz+PSnP80JJ5zA6NGju+2zceNGvvCFL7B06VLq6up48smts/odcMABncdMnDiRlStXssMOO7Dzzjt3/k5tv320Ctvvfvc7Hn30Ue644w4A3njjDZ566qmKHY+cVMmDqqTjge8DOwH/J2mpmR1B1IHoHklbgL8Bn40d9l1JE8LrWbZ1guR7gMMlLQM2A181s1fDdZaa2cTYOZqAjycpo5m9Jmk20BG75mt532xfzJy59Y9OgtQfzMyZHmyrVHt7O3PnzqWlpYW5c+dG/1z1MbCaGaeeeiqXXnppl/xvfetbncMb6urqOp+TbtmyhT/84Q9kWsJw2223js3P9AGecvrpp3PjjTfy0ksv8fnPf75P5a942f6OC2zFihXU1dUxatQoli/fOjPfBRdcwNFHH82vf/1rDjzwQObPn9/t2Msvv5x3v/vdPPLII2zZsoVBgwZ1bnvXu97V+Tr1e2JmGYfGmBnf//73OeKIIwp8d+Wt5B2VzOwXFs3h+C4ze3cIqJjZSjP7NzPb28w+ZmbPxY6ZYmbjQ7otlm9mdl7I/0Datolp193DzP6So1yTzWxx7P31ZrZnSDcU6v77bObM6A8z9ceZeu0BtaLFn63PmjUr87P1Xjj00EO54447eOWVqGP8a6+9xnPPPZd1/8MPP5wrr7yy8/3SpUsz7vfv//7vnf0Afve737F27dbRZscffzy//e1v6ejoqLkP2FJYvXo1Z511Fl/4whe6BbtnnnmGD3zgA3zta19j0qRJ/OUvf2Ho0KGdrQ8Q1Sh33nlnBgwYwE9+8hM2b96c83rjxo3jxRdfpKMjqnOsW7eOTZs2ccQRRzB37lw2btwIwJNPPlkTS2GWvKbqemHGjFKXwBVZR0dHl/HJqfHLHR0dfaqtjh8/nv/5n//h8MMPZ8uWLWyzzTZcddVVWff/3ve+xznnnMM+++zDpk2bOPjgg7n66qu77TdjxgymTJnC7bffzkc/+lF23nlnhg6NlrkcOHAgjY2NDBs2jLq6ul6XveoU8O94w4YNTJw4kY0bN1JfX89nP/tZzjuv+7rdV1xxBe3t7dTV1TF+/HiOOuooBgwYQH19PRMmTOC0007j7LPP5sQTT+SnP/0pjY2NXVokMhk4cCC33347X/ziF9mwYQODBw9m/vz5nH766axcuZL99tsPM2OnnXbizjvvLNg9lyuf+7eICjn3b2trKw0NDV0+UNvb2+no6Nj6HDbetOTKzvLly9l7771LXYyieOedd6irq6O+vp4//OEPNDc3d9Zqt2zZwn777cdPf/pTxo4dm/H4Sv7eVHLZy0Gm75/P/euKLtEQCw+orkSef/55GhoamDBhAueeey4//OEPAVi2bBl77rknhx56aNaA6lw18ebfChGfvrC5uZm5c+d2Ngemtqd0q8E6V2Rjx47l4Ycf7pY/fvx4VqzI1gHfuerjNdUKEh9i0dzcTGNjY1EnCXCF549buquG70k13EMpVOP3zYNqBUkfYtHe3t6lBjt9+nSfgL+MDRo0iFdffbUqP0h6yyxaUzM+bKPS+M+1d6rhZ5+JN/9WiPgQi9R4xfj7Qk8S4Apv9OjRrFq1itWrV/e8cw0ZNGhQxkkIKoX/XHuv0n/2mXjv3yLqr96/qSbg+LNWD6zOuUpVyb1/PagWUSGDajbpNdj09845V2kqOaj6M9UKl2uSAOecc/3La6pF1B81VeecqzZeU3XOOeecB1XnnHOuUDyoOueccwXiQbVKtba2dlsmrL29ndbW1hKVyDnnqp8H1Srl0xc651z/86BapbJNX9jR0eE1WOecKxIPqlXMJ+B3zrn+5eNUi6jU41RTATN9+sJs+c45Vw58nKorO/HpCmfNmtXZFJxa2Sa9Busdm5xzru88qFapXNMXZlpCzpuFnXOu77z5t4hK3fybSa4J+AFvFnbOlZw3/7qKkasGm6lZ2DnnXHJeUy2icqyp5uIdmJxz5cBrqq7i5erY5JxzLhkPqg7wdVmdc64QvPm3iCqt+dc558pB0uZfSXXAYuBvZnZMln0agIeAT5nZHYUtaXdeU3XOOVepvgQsz7YxBN05wD39VSAPqs455yqOpNHA0cB1OXb7IvAz4JV+KRQeVJ1zzpWfekmLY2lqhn2uAKYBWzKdQNJ7geOBq4tYzm7q+/NizjnnXAKbzGxSto2SjgFeMbMlkiZn2e0K4GtmtllSMcqYkddUXU4+J7BzrgwdBBwraSVwG3CIpJvT9pkE3Bb2OQn4gaRPFLtgJQ+qkj4p6QlJWyRNiuUPlHSDpMckPRL/b0TSpyQ9Go5rTTtfk6RlYdstGa43VNLSWFoj6Yq0fU6SZGnl2Rw75q6CfhPKmM8J7JwrN2Z2oZmNNrMxwMlAm5l9Jm2f3c1sTNjnDuBsM7uz2GUrh+bfx4ETgGvS8s8AMLMPSBoF/CZ0jR4OXAbsb2arJd0k6VAzWyBpLHAhcJCZrQ3HdWFm64CJqfeSlgA/j70fCpwL/DHt0A1mNpEaE1/s3Gdacs6VM0lnAZhZvz5HjSt5UDWz5QAZ2rzHAwvCPq9Iep2oOm/Ak2a2Ouw3Hzgx7HsGcJWZrU0dl+vaIQiPAhbFsmcDrcD5vb+r6hKfE7ilpcUDqnOubJjZQmBheJ0xmJrZaf1VnpI3/+bwCHCcpHpJuwP7A7sATwPjJI2RVA98IuQD7AXsJekBSQ9JOrKHa0wBbrcwA4akfYFdzOxXGfYdFHqhPdQf7fLlJNNScc4557rrl5qqpPnAezJsusjMfpnlsOuBvYlmy3gOeJCoR9haSc3A7URdqR8E9gjH1ANjgcnAaGCRpPeb2etZrnEy8NlQxgHA5cBpWfbd1cxelLQH0CbpMTN7JsO9TgWmAgwcODDLqfrZzJlR6oX0peIaGxu7vHfOObdV2UxTKGkhcL6ZLc6y/UHgdDNblpY/FdjTzKZJuhp4yMxuDNsWABeYWbcJbCVNAH5qZnuF9zsAzwBvhl3eA7wGHJteJkk3Ar/qacqrspmmUIJe/pxbW1tpaGjoEkDb29vp6Ohg2rRphSqhc851quRVaso2qEoaQlS+9ZIOA1rM7OCwbVR4zjocaAeazOzJ0Nw7xcxOlTQSeBiYaGavZrjeN4F3zGxGT+UJ13nLzN4J5/0DcFx6gE9XDUHVOef6WyUH1ZI/U5V0vKRVwIeA/5OUmqNxFPBnScuBrxGaaYPvSloGPAB808yeDPn3AK+Gbe3AV1MBVdLStEs3AbcmLObewGJJj4TzfrOngFpyM2dGwTTVASz1upfNwM4553pWNjXVauQ1Veecy5/XVJ1zzjnnQbUmzMj42Ng551yBefNvEZVN869zzlWQSm7+zTlONUyucCzRmnUTgGHA60QTM/wGuNPMNhW7kM4551wlyFpTlXQmcBHRqur3ha/rgKFEvWE/Gr5eUsp5FsuZ11Sdcy5/1VpT3Qs4wMxeyrDtF8AlknYGvlKUkrni68NMS84557rzZ6pFVPY1VZ9pyTlXhiq5ppq496+kvSW1SLoqvB8naZ/iFc2VM19n1TlXDSSNlHSepAVhfe2N4esCSedL2imf8yUKqpI+SfRc9b1sndloO+A7eZXelV6BZlqKr7M6ffp0n2TfOVdxJF1KNJ3tvwE/Ag4j6it0WHg/lmhmv28mPmeS5t8wVeAUM1sqaa2ZDZe0DfCimeUVxWtJNTf/pkyfPr1zndVZs2YVqGDOuVrWX82/kr4IXGtm7+TYZxDRYi5XJjln0ubfUUTDaCBaJDz11R/I1rBM66y2trZ2W281le+cc+XEzL4fFkqpk/R5Se/KsM/bSQMqJA+qS+g6oT1Ea5H+KemFXBnqw0xL8XVWZ82a1dkUXF9fn/FZ6zPPPOPB1jlXlsxsM/CdXDXWpJI2/44Dfgc8CxwILCQacnO4mT3V10JUq7Jv/u2DXL1/U52YmpubmTt3LvPmzQPo8tw1HpRTx3hPYucclKb3r6SfAPPM7O4+ncjMEiVgCNFyaV8lqqVul/TYWk1DhgyxWtXS0mKAtbS0dOa1tbXZyJEjraWlxUaOHGltbW1d8rO9d87VFmC99fPnNfBT4B2iSuNPgB+nUj7nSTykxszeIlqce5GZ3WZmb/YpmruqlelZK0Q9hpubm5k9ezbNzc2dNdNsPYk7Ojq8ydg5118eBy4hWjP7aeCZWEouYQTflWhB8PXAmyHvJOC6/v5vopJSxdZUZ8zo9aG5ap3Zaqop6bVbr8E6V5soQU21UClpUP0N8N9EHZvWhrwdgOdKfQPlnCo2qEKvD50zZ063oNfW1mZTp05N1OSbrWk4np/tGnPmzOl1uZ1z5aNUQZVovGoT8Pl4yuscCS/0KjAgvH4tlv96KW68UlItBtVscgXCnmqkXoN1rraUIqiGiuPbwB9DE3AqteV1noQXWwbsZbGgCowHHu3vG6+kVFFBdcaM6NchPfWhKTipJAE3SQ3WOVcdShRUXwH26fN5El7s88CTwOeAfwBTgMeAT/f3jVdSqqigGleEmmpv5FuDdc5VhxIF1eeAgX09T+JVaiR9ApgK7AY8D1xjZncmOrhGVew41QJMX1gI+Y6F9fGuzlWHfpymMD4C5jPAQcBM4OX4fma2JfFJE0TvOqJxO+/q7/8cKj1VbE21H5p8+yJbDfbb3/52xvypU6d6xybnKgj9VFMFtgCbQ9qS9j6Vtzmfc/Y4TtWi6Zt2J49l4lyFK/OFyzs6OrqsiJMa57pp06aM411PPvnkvKZO/PjHP+7jY52rDbsDe4S0e9r7PWLvk0sYzT8P3EjU9FtHFGAHEHoEe6qymmqFSzqbk9d4nStPlHCcaohtO/c2viW9SEGqxbWWPKj2v1y9gnszdWKSIOzB1rnCKkVQBbYnmpbwnyG+vQPcBOyQ13kSXmy3bKm/b7ySkgfV/tXb2Zyy9SLua43Xh/k41zslCqo3Aj8jWizmXeHrT4Gb8jpPwoudnyX/vP6+8UpKHlT7V29mc+rNONikwdZnfnKud0oUVF8ChqTlbQe8nNd5El7sH1nyX8vnYrWWPKiWh3yDbbZnqvnWeL252LneKVFQXZne+gqMAZ7P6zw9XOSQkNYDjbH3hwCn43P/elCtYNmC7VFHHVXUGq83FzuXW4mC6teJJjk6CzgqfP0r8PW8ztPDRZ4NaXPs9bPACuBB4Nj+vvFKSh5Uq0tvFwvIp4OUNxk7lzyohtEoDwO/yrDt08CjIT0ITOjhXAojXeYTTc07H/hPiCZJSpqSRvD/zeeknjyo1pLezF1sljvYepOxq2V5BNXzgFuyBNUPA8PD66OAPyY5Z19TkkLXARvwGZU8qLq89LY3sjcZu1qXJKgCo4EF4XFkt6Catu9w4G85tgt4T6pWChwBfA+Y2lM5up0r0U7wCPAv+Z484bk/CTxBNC5oUix/IHAD0cT9jwCTY9s+Far0TwCtaedrClX3J4BbMlxvKLA0ltYAV4RtpwGrY9tOjx13KvBUSKcmuTcPqrWtt83FZt7D2NU2ojGii2OpW3AD7gD2ByYnCKrnA9dl2XZwiANbgKeJFoz5GzAPeBGYnevc3c6XaCeYBvw5BJZDiXVYyudiWc69N9HCsAvTguo5wA3h9ShgCdFMFyOIJvTfKWy7CTg0vB5L1L6eqvKPSnD9JcDBtjWoXplhnx2JniPvGP7jWZG6Rq7kQdVl0lMgLEQPY6/BukpGDzVV4BjgB+F1zqBK1Ml2OTAiy/YOomeng4k6J70FjA/bxgErc5Wl2/kS7dS1k1KXDkv5XKyHa6QH1auAz8TeLwAOABqA+bH8z8a+ua3x2mWCa44FXmBrlT9bUJ1CtCpP6v01wJSezu9B1eWrN03GXoN11SZBUL0UWEU0DOalEAhvzrDfPsAzhPXAs5zrjdjrOmBD2vZ1ucrS7Xz57FzMlCGoTiWazaKeaFLj14ETQ01xFdH4oXqiGTDuDsfcGQLrA8BDwJE9XHM68K3Y+9OAvxM1Ld8B7BLyzyfWrRpoIcuEGPHkQdXlq1A9jL3Dk6tkPQVV6/o5nrGmCuxK1Jz74R6O/0fa+9dybe8p1dMPJM0negic7iIz+2WWw64nahpeTLR47IPAJjNbK6kZuJ2oDfxBtq4iUE9U+5xM9BB7kaT3m9nrWa5xMlFNN+Vu4FYze0fSWURNy4cQPcROZ1nudSrRPwQMHDgwy2WdyyzTuq+NjY1ZV+bp6OgAYO7cubS0tDB37lwaGxs7t6evOQt0rt7T2NjYuWJPaptzlSx8bmNmVxNVmkYAP5AEUfyYlOGwd0maFXs/OO19fh/kCf8T2B74DtHzx+eInmk+T54zTfRwjYXEaqoZtj9IaOdOy59K6KwEXA2cFtu2AGjIcr4JwJM5rldHaBbAm39dmerpmap3eHKViH6c/IGoQ2zOlNf5El705hD0jgPWha/3A18u4I11CarAEGDb8Pow4PexbaPC1+FEvXT3Cu+PJEx+DIwkel6a7eH0N4GL0/J2jr0+HngovN6R6Bny8JCeBXbs6Z48qLpiK9QYWe/w5MpJfwbVQqekAe+VVHACXg9f3wv8uc8FiILXKqIu1C8D94T8MURTRC0nmtlit9gxtxINm1kGnBzLV6hRLyMaihPftjTtuiuAcWl5lxINxXkEaI9vJ5pp4+mQPpfk3jyoulLxDk+ukvVXUCXBCJGw37sTnzPhCdcA9eH1KmAY0fCWvB7g1lryoOpKpdgdnrwG64qpH4PqE8APgA+Rtih5iHEHhu2PJz1naihJTpIWAJeY2QJJtxJ1EHoT2N8yP/h1wLbbbmvr168vdTEKZ+bMKLmK1draSkNDQ2eHJ4D29nY6OjpoaGjo1rEp3pkpPd+5YpH0lplt2w/XGUjUL+dMog6vK4gecQ4N758i6kPzIzP7Z6KTJozmewD/Gl7vBFxH1Pu2W8chT1VcU4VSl8AVSb4dnrxZ2BUTpVmlZhfg48B/EM0V/N7enCdRTdX1TtXVVCXw35eqlG8NFrIPzfFarOur/qqpFkPOoCrpkJ5OYGZtBS1RFamKoDpzJlx8cff8GTO8KbgGpAfL9HGt3izsiqGag+qzaVm7EA1TSTEz2wOXUVUE1TivqdacXDXYadOmMX36dGbPnk1LSwuzZs3KcSbnkqvaoNptZ2mtmQ0vYnmqigdVV80ydWBKNRdnC8LOJVHJQXVAnvv7J2otmzGj1CVwZSLeDDxr1qzOKRHr6+tpamqivb29y34NDQ0lLrFzyUgaIGnnXh+fZ031NTPbsbcXqzVVV1N1LujN0BznkipFTVXSMKIxqScBG81sW0nHAgeY2dcTn8eDavF4UHW1yp+1ur4oUVC9DVgLzAKWmdlwSTsBD5rZ2KTnyblKjaRFdG3yHSrp9/F9zOzg5MV2zlW79vb2jKvmOFfmDgX+xcw2SjIAM1staVQ+J+lp6bfr0t7/KJ+TO+dqS/oQnMbGRh+/6irFG0QLsfw9lSFp1/j7JHIGVTO7qVdFc87VpFzrvnpQdWXuOuBnki4CBkj6EHAJ0ZKiiWV9pirpWDO7q8cTJNyvFvkzVeecy1+JnqkK+BLRXMC7Ea0Zfg3wXcuj81GuoHoLsA/wv8B9RMuwpSYa3gv4KPAZoiXVPtPrO6liHlSdcy5/VTlO1cz+A5hCtG7qT4DVwAaitVVvAt4DfMoDqnOuJ62trZ1jV1Pa29tpbW0tUYmc60rSBZIa0vIOkJTXrCU5J38ws8fM7Atm9q/AdkTTFA41s73M7L/M7Im8S+6cqzmpsas+KYQrY18ClqXlLQP+K5+T+Co1ReTNv85t5euyuqRK9Ez1VWBni62bGtZbfSmf+RnynabQOed6pbGxkebmZmbPnk1zc7MHVFdulgBnp+WdBfw5n5N4TbWIvKbq3FZeU3VJlaim+j7gXqJxqc8AewLvBg4zs/Rm4ezn8aBaPB5UnYvkWpfVA6tLV6rev5K2A/4fMJpomdNfmdmb+ZwjUfOvpHMljcy/iM45l3tSCOfKhZm9aWa3mtllZnZbvgEVEtZUJd0FHAIsJBpec6eZvZPvxWqN11Sdcy5/JWr+3R34BjCRaLRLJzPbNel5EhPR/2kAABvrSURBVNVUzexYohkmfkPUvfglSddJ8sn0nXO95uNXXRm5BdgCfAX4bFpKLHHvXzN71cyuMrMPEc2m1AC0S1op6aLQFu2cc4n5+FVXRt4HnGJmvzGz++Ipn5PkNaRG0qGSbiBqBn4ZOIUoiu9LVIt1zrnEUs9Wm5qamD59undecqX0e6JY1ic9Lf0GgKRvAScTLY3zY+DrZva32PaHiBZ3dc65vMTHr7a0tHhAdaWyErhH0s+Bl+IbzGx60pMkCqrAIOB4M8vYVS8s6jop6UWdcy7FFzV3ZWJb4G5gG6IpeXvFx6kWkff+dS63bONXTzjhBE4++eQuwbW9vZ2Ojg6mTctrfnNXgapylZo4SYsk/T5DulfSDZL+X7EL6pyrPtnGrwLegcmVhKShknaXtEcq5XV8wnGqs4FTiZZ8e4GoanwKURdkAf8JXGZm3g8+xmuqzvWeT2tYu0o0TnU80frhEwAjim0GYGZ1Sc+T9Jnq4cARZrY8VoD/BW4ysw+GB7u3AR5UnXMF4R2YXD/7AdAONALPAmOAS4EH8zlJ0prqG8Co+CxKkgYDfzezYeH9m2bmY1VjvKbqXO95TbV2laimupYozm2U9LqZDZO0LfC4me2e9DxJx6n+HrhB0p6SBknaE/ghcH8ozAeIZvbPm6RPSnpC0pZ4D2JJA8Pz2sckPSJpcmzbpyQ9Go5rTTtfk6RlYdstGa43VNLSWFoj6Yqw7TRJq2PbTo8dtzmWf1dv7tU5l0y8A9OsWbM6x7Kmz77kXAG9TdTzF2CNpF2JYuSIfE6SNKieGvZdBqwHngDqgNPC9n8CU/K5cMzjwAlEgTvuDAAz+wBwGPBtSQMkjQAuAw41s/cB75Z0KICkscCFwEFhW7cV281snZlNTCXgOeDnsV1uj22/Lpa/IZZ/bC/v1TmXgE/A75KQVCfpYUm/yrBNkr4n6elQCduvh9MtAprC6zuIJjS6D2jLp0w9PlOVVEcUnE4D/gPYCVhtZltS+5jZX/O5aFzqOa2k9E3jgQVhn1ckvQ5MInpw/KSZrQ77zQdODPueAVxlZmtTx/Vwb2OBUUTfTOdcmcg0bMbHr7oMvgQsB7bPsO0oYGxIHwTmhq8ZmVlT7O1/E1X4hhJ10E2sx5qqmW0GzgH+aWZbzOzleEAtokeA4yTVh9UD9ifqdfw0ME7SGEn1wCfYOlB3L2AvSQ9IekjSkT1cYwpRzTT+YPnE8F/NHZLiA4AHSVoczvuJgtyhcy5vPgm/A5A0GjgauC7LLscBP7bIQ8AwSTvnON/5qdch1t1sZnOBs/IpV9Lm35vyPXGcpPmSHs+Qjstx2PXAKmAxcAVRD6xNoRbaDNxOVMNcCWwKx9QT/VcymShgXidpWI5rnAzcGnt/NzDGzPYhqgHH/0PZ1cwmEdXWr5D0r1nudWoIvos3bdqUaRfnXB/4JPw1oT71ORrS1Az7XAFMI1pZJpP3Eg0BTVkV8rLJNhXh13ssbUzSITUHAF+UNI2okJ01OzPrcfk3M/tYPoUKx2wCvpx6L+lB4Kmw7W6iAEj4Zm8Ou60CHjKzjcCzkv5KFGS7PYiRNAGoN7MlsWu+Gtvlh8Cc2LYXw9cVkhYSTbz8TIZyXwtcC1Hv33zv2zmXW3wSfu8ZXLU2hUpMRpKOAV4xsyXxTqzpu2XI6/aZLOmQ8LJOUmPacXsA65IVOZI0qP4wpH4jaQjRkJ/1kg4j+iYvC9tGheesw4Gz2fpw+U6iGuqNkkYSNQevyHKJKXStpSJpZzNL9WI+lqitnnCdt8zsnXDeg/Axuc6VjI9hrXkHAcdK+jjR3PTbS7rZzD4T22cVXefwHQ28mOFcPwpfBxG1kKYY0WpsX8yrZGZW0gQcH27+nXAD94T8McBfiQLbfGC32DG3EvVEXgacHMsX8J2Q/1jatqVp110BjEvLu5SoZ/MjRIOAx4X8D4fzPRK+/meSexsyZIg55wqvra3NRo4caS0tLTZy5Ehra2uzOXPmWFtbW7f95syZU6JSut4C1lvyGDIZ+FWG/KOJevAKOBD4Uw/n+XHSa+Y8T8JCi6hnbRvwaMg7GGgqRCGqNXlQda7wUgE1FUBT77/97W9nzE8PtK789TaoEvX9Ocu2xq2riB7TPQZMSnrOcHwj8JF8jjGzvOb+PYzowfDVFs00sQfwUzPbv8cT1CifUcm5wmttbaWhoSHjCjapTkz+rLWylWhGpfuA/zazByR9DTiPqBPsVWZ2SeLzJAyqLwD7mtkaSWvNbLiigaWvmdnwXt5D1fOg6lz/mz59euez1lmzZpW6OK4XShRUXyWapnCzpKeB/we8CTxgZrsmPU/SITV14eSwtffUdrE8V8tmzix1CZwDui947tMaujwMACwMl5SZLTezF4C8Ko5Ja6rXEU1F+GWiOX5HAJcDA83s7HxLXitqpqYqQYLfI+eKKduC594EXHlKVFO9m2jI6M7AM2Z2fgiw860IE+qfB/wL8AawA1ENdTfga3mV2jnnisTnC3Z9dBrwOvAoMDPkjQO+m89JEtVUO3eWRhEF0xfM7KV8LlSLqrqmOnMmXHxx9/wZM7w52DnXJ6WoqRZKb4JqlzVTzSzb5Ao1r6qDapw3/7oylqu3cKaJ+13p9VdQlXSRmX0jvM7aq83Msk1h2E2i5l9JR0r6G/AS0YT2qfRU0gs551wp+FzBLofRsde75EiJJe2o9AzRGqY3mdmGfC5Qy2qmpjpzpjf5urKWCqQ+frUyVH3zr6TXgBGWT1uxq52g6lwF8PGrlaNEvX/HAx8BdgReAxZZmG8+H0l7//4I+Fy+J3fOuXLg41ddNopcTzSV4X8TLaZyEfCopBvCREeJJV2l5kDgXEkXED1X7WQJln5zzrlSSR+v2tjY6ONXXdxUovmDDzSzzvFXkhqIFm85E7g66cmSBtXryL66unPOla1s41cvu+yyzvcp3iu4Jn0WODceUAHMrEPSfwEXkkdQzWtIjcuPP1N1rnz5DEzlqz+fqYY+Q7uZWbfFyCUNBZ7PZ477nM9UJX0v7f1/pr3/WdILOedcwfWh13mqxtrU1MT06dM7A2pHR0e3Z67t7e20trb2sbCuTNVlCqgAIT9p3yNIsPNpae8vS3t/WD4Xc865gso0q1ceGhsbaW5uZvbs2TQ3N9PY2OjjWmvPNpIaJR2SKZH8MSkk2Dm911NevaCcc66cpfcKTnVkStVgfVxrTXgFuL6H7cn1sPL5P9Lev5Zru6euaciQIeacK7AZM8yiiTG7phkzuu7Tg7a2Nhs5cqS1tbVlfN/S0mKAtbS0FP4eXE7AeiuDz/DepJwdlSS9BRzN1hrqncBxsfd3W4XOetEfvKOSc0WWbd7pBPNR55oTONUEHK+ppvK9t3DxVe2MSpJWsnVR8owsj3Xmao0HVeeKrA9BNZtsvYIvvPBCLr30Uu8t3A8qOajm7KhkZmPMbPdcqb8K6pxz3cyYsfX1zJlRME1NgJN6nWcP4WzjWjdt2uS9hV3PSt3+XM3Jn6k6VyJQtFOnP2vt6dmsyx8V/Ew1r/E3zjlX0fq4mlKmOYSzjXf1JuEaVeqoXs3Ja6rOFUCCnryJj+lDDdZ7C/cfKrim6tMUFpF3VHKuAPrQ6aiQ58q3t7DXVHuvajsqOedcxStQB6Zp06Z1C5TxGZjmzZvHrFmzOpuCzzzzTO/AVIM8qDrnyk+BAmHnucy21lBTr/v4fDUlW29hwKc7rEHe/FtE3vzrXAGUSfNvb6QCqTcL58ebf51zrhLEx7X2g0wT9rvq5kHVOVfeChkIC9Tkm1SmITiuunnzbxF5869ztcsXQe89b/51zjnXRbYOTB0dHSUumSumkgdVSZ+U9ISkLZImxfIHSrpB0mOSHpE0ObbtU5IeDce1pp2vSdKysO2WDNcbKmlpLK2RdEVPx0s6VdJTIZ1a8G+Ec66qZBuCA/Q81Kafm6ldAZV69glgb+DfgIXApFj+OcAN4fUoYAnRPwEjgOeBncK2m4BDw+uxwMPA8NRxCa6/BDg41/HAjsCK8HV4eD28p3P7jErOuXSJ5gou4tzFlYAKnlGp5DVVM1tuZn/NsGk8sCDs8wrwOjAJ2AN40sxWh/3mAyeG12cAV5nZ2thxWUkaSxSwF/Vw/BHAvWb2Wth2L3BkvvfqAv8v3FW6PvwO+1zB1a3kQTWHR4DjJNVL2h3YH9gFeBoYJ2mMpHrgEyEfYC9gL0kPSHpIUk+Bbwpwe/jPKNfx7wVeiB23KuS53rj44uzbPOC6SpDrdziBjENtCjnhhSuZfgmqkuZLejxDOi7HYdcTBa/FwBXAg8CmUFNsBm4nqmGuBDaFY+qJmnAnEwXM6yQNy3GNk4FbY++zHa8Mx2bsNi1pqqTFkhZv2rQp0y4ulz5+WDlXUgkDYMahNkWe+amaSBok6U+hv80Tkrp9cEjaQdLdsX0+1x9l65egamYfM7P3Z0i/zHHMJjP7splNNLPjgGHAU2Hb3Wb2QTP7EPDXVD5REP6lmW00s2fDtrGZzi9pAlBvZkti2dmOX8XW2jDAaODFLOW+1swmmdmk+vr6Hr83NcP/C3eVLsnvcIJ/CuNDa2YNGNDZFJxoDKv/vaS8AxxiZhOAicCRkg5M2+ccYFnYZzLwbUkDi12wsm3+lTRE0rbh9WFEtdRl4f2o8HU4cDZwXTjsTqAxbBtJ1Jy7IsslptC1lprr+HuAwyUND9c8POS5pHL9F57kw8o/TKpfKX/GSa5doJpkl6E2F1+ceahNtgkvvCUH6OzF9WZ4u01I6a2HBgyVJGA74DW2tmoWtXCl7v17PFFN8B3gZeCekD+GqKa4nKgz0m6xY24FloV0cixfwHdC/mNp25amXXcFMC4tL9fxnyd6nvs08Lkk9+a9f7PI1bMx27Ya7w1ZE0r5M8732vH9Z8xIhdeuKb6ma4L1XefMmdO1B7BFPYPnzJmTrJy9WXe2TIV4sDiWplr32FEHLAXeBOZk2D4UaAf+HvY5On2fYqSSB9VqTh5Us8j1x+9BtXaVY1DN9rua7yLoCYLwilNOyTjUZsUpp2Tcv1sZquhvhDyG1BA9GmwH3p+WfxJweags7Qk8C2yf9Ly9TSUPPNWcPKj2Qvp/90k+TFzlKuXPOMm1+1KD7UV+KpC2tLR0H7vaU3lqNKhGuzMDOD8t7/+Aj8TetwEH5HPe3qSSB55qTh5UC6iKPjBcFuVYU823TPn+U5jh/C0tLQZYS0tLz+Wp0n88ewqqwE7AsPB6MNFIkGPS9pkLzAyv3w38DRiZ67yFSCUPPNWcPKgWkAfV6lcuQbUYgSph83KPNdXePDqpQAmC6j5Es989CjwOTA/5ZwFnhdf/Avwu9I95HPhMrnMWKpU88FRz8qBaQBX+n7dLoJQ/43yfkeYrwXmyTV84derUnjswFbKsZSDf5t9ySmU7pMa5LnxITfUr9yE1fZFgTdhsq9oAXcaxpsa5NjQ0JLuG/+30K19PtYh8PVXnKlxqHHWJpQJpc3Mzc+fOzW+uYGnr2NoK4eupOlcqZfCB56pYmfx+ZZwr2JUlD6qusvkMM5WnTAJVJck4VzDQ2tqacW3WBw47zKcFLZVSP9St5uQdlfpBFXXOqBn+M8tLrvVX+7Q2axl3/sM7KjnXj3xyfldDsnVg6ujo6NvarN7KUxyljurVnLym2g9KVesp4//yiybfKft8dqx+k3PCiGIPFyoCKrimWvICVHPyoNoP+jLjTSVet5T6OAVfj/muV3qcMCKuQv658aDqKWPyoNoPCjXheb4fKvmepxwDSaHu2YNqySR6pppNGT9rreSg6s9UXWXL9hw13+dF2fZPX9O1AItUl42+3PPkyfnlx8+ZYCIEl0yu5629Vkm/w2XIJ38oIp/8oYTiA95nzsz8QTFjxtYP+2wD5PuSn+S6pZxcoBj3nCTflYdsv3tl8HPzyR+cKwfZalWw9elR+utC9CLO97pJa7bZytGb8vVUVu85XXvy/b3w35FkSt3+XM3Jn6mWUF+e8yXpzNGXRaqT5PfmmCIsnJ3onvPNd0U3Z86cZJPwZ1IGz8Kp4GeqJS9ANScPqiVUjMBTqOsm7YFZ7I5BZfDh6YqjTyvelMHvhQdVTxmTB9USyreWVKheu33tUZst4H70o5nzsy14XYzatqsomYbaJOotXAbjiz2oesqYPKhWgWI3Yfa1+TfJh14tjql1ZpZ5Uoi8xrXG9eMQHA+qnjKmmg+q/uHcs1zfI2/OdX2QK3imB9tEz2D78ffLg6qnjKnmg6p/mPdNqSa2cBUvyST8SZqFuzyDDb9HiYNtH3hQ9ZQxeVCl1CWoLR48XZCt5jl16tSCBNsVp5xixXzW6kHVU8ZUk0G1QuYWda4W9dTMm/QZbJfzhH+eEw/ZScCDqqeMqSaDapzXVJ2rGPk8g+1Sg4X8mosT8KDqKWPyoEqpS+CcSyDfZ7Bd9jn44MzNxbFgm7iHceBB1VPGVPNB1Zt8nasIvXkGa5a7udigVwHVzCo6qPqE+kXkE+o75ypZa2srDQ0NnavgALS3t9PR0UFDQwNNTU00Nzczd+7caLWc++7reRGJBCp5Qn0PqkXkQdU5V43a29tpamrqXHYu/h6gqamJ1WvWsNPIkV2WpkuqkoOqr1LjnHMuL9nWcb3tttu6BNd58+bR1NREe3t7KYvbr7ymWkReU3XO1ZIuzcVhvdZUc/G0adMSn6eSa6oeVIvIg6pzzuWvkoNqyZt/JX1S0hOStkiaFMsfKOkGSY9JekTS5Ni2T0l6NBzXmna+JknLwrZbMlxvqKSlsbRG0hU9HS9pc+yYuwr+jXDOOVfx6ktdAOBx4ATgmrT8MwDM7AOSRgG/kdQADAcuA/Y3s9WSbpJ0qJktkDQWuBA4yMzWhuO6MLN1wMTUe0lLgJ+H17mO32BmE3HOOeeyKHlQNbPlAJLSN40HFoR9XpH0OjAJMOBJM1sd9psPnBj2PQO4yszWpo7Lde0QREcBi0JWXsc755xzcSVv/s3hEeA4SfWSdgf2B3YBngbGSRojqR74RMgH2AvYS9IDkh6SdGQP15gC3G5bHyznOn6QpMUh/xOFuknnnHPVo19qqpLmA+/JsOkiM/tllsOuB/YGFgPPAQ8Cm0KzbDNwO7Al5O8RjqkHxgKTgdHAIknvN7PXs1zjZOCzsfe5jt/VzF6UtAfQJukxM3smw71OBaYCDBw4MMtlnXPOVaN+Capm9rFeHLMJ+HLqvaQHgafCtruBu0P+VGBz2G0V8JCZbQSelfRXoiDZkX5+SROAejNbEsvOeryZvRiuvULSQmBfoFtQNbNrgWvDNbZI2pDvvcfUA5v6cHwlqrV7rrX7Bb/nWtGXex5cyIL0p5I/U81G0hCiIT/rJR1GVEtdFraNCs9ZhwNnA03hsDuJmnRvlDSSqDl3RZZLTAFuTcvLeHy4zltm9k7IPwhopQdm1qfmdUmLzWxSz3tWj1q751q7X/B7rhW1eM9QBkFV0vHA94GdgP+TtNTMjiDqQHSPpC3A3+jaTPvdUNMEmGVmT4bX9wCHS1pGVHv9qpm9Gq6zNK33bhPw8bTiZDxe0oeBa0JZBgDfTAV455xzLsUnfyhjtfifXq3dc63dL/g914pavGco796/LjybrTG1ds+1dr/g91wravGevabqnHPOFYrXVJ1zzrkC8aBahiQdKemvkp6WdEGpy1MMkq6X9Iqkx2N5O0q6V9JT4evwUpax0CTtIqld0vIwt/SXQn7V3rekQZL+FObvfkLSxSG/au8ZQFKdpIcl/Sq8r/b7XalonvalkhaHvKq+52w8qJYZSXXAVcBRRFM1TpE0vrSlKoobgfQZry4AFpjZWKJpJ6vtH4pNwFfMbG/gQOCc8LOt5vt+BzjEzCYQzbl9pKQDqe57BvgSsDz2vtrvF6DRzCbGOifVwj1340G1/BwAPG1mK8zsn8BtwHElLlPBmdnvgdfSso8DbgqvbyKagrJqmNnfzezP4fU6og/d91LF922RN8PbbUIyqvieJY0Gjgaui2VX7f3mUIv37EG1DL0XeCH2flXIqwXvNrO/QxSAiMYqVyVJY4hm5fojVX7foSl0KfAKcK+ZVfs9XwFMI5pGNaWa7xeif5R+J2lJmOUOqv+eMyr55A+um27L9RD9wroqIWk74GfAf5nZPzKs0FRVzGwzMFHSMOAXkt5f6jIVi6RjgFfMbIlia0DXgIPC3OijgHsl/aXUBSoVr6mWn1VsXXUHoon9XyxRWfrby5J2Bghfq27pPUnbEAXU/zWzn4fsqr9vgLAwxUKiZ+nVes8HAcdKWkn06OYQSTdTvfcLQGxu9FeAXxA9xqrqe87Gg2r56QDGStpd0kCilXTuKnGZ+stdwKnh9alAthWMKpKiKumPgOVm9p3Ypqq9b0k7hRoqkgYDHwP+QpXes5ldaGajzWwM0d9um5l9hiq9XwBJ20oamnoNHA48ThXfcy4++UMZkvRxoucydcD1ZvaNEhep4CTdSrTE3kjgZWAG0YIG84BdgeeBT5pZememiiXp34FFwGNsfd7230TPVavyviXtQ9RJpY7on/h5ZjZL0giq9J5TQvPv+WZ2TDXfr6LlMH8R3tYDt5jZN6r5nnPxoOqcc84ViDf/OueccwXiQdU555wrEA+qzjnnXIF4UHXOOecKxIOqc845VyAeVJ1z3Uh6MwyVcM7lwYOqc2UoLKX1MUmnSbq/yNdaKOn0eJ6ZbWdmK4p5XeeqkQdV56qYJJ/f27l+5EHVufK1N3A18KHQHPs6gKR3SfqWpOclvSzp6jAFIJImS1ol6WuSXgJukDRc0q8krZa0NrweHfb/BvAR4MpwjStDvknaM7zeQdKPw/HPSfq6pAFh22mS7g/lWSvpWUlH9ft3yrky4UHVufK1HDgL+ENojh0W8ucAexEt+r0n0dKA02PHvQfYEdgNmEr0d35DeL8rsAG4EsDMLiKaOvEL4RpfyFCO7wM7AHsAHwVOAT4X2/5B4K9EU062Aj9StS+941wWHlSdqyAhWJ0BfNnMXguLnV9CNHl7yhZghpm9Y2YbzOxVM/uZmb0V9v8GUXBMcr064FPAhWa2zsxWAt8GPhvb7Tkz+2FY4u0mYGfg3X28Vecqkj9vca6y7AQMAZbEKoMimrA+ZbWZvd25URoCXE605NrwkD1UUl0IhLmMBAYCz8XyniOqHae8lHphZm+Fcm2X9IacqyZeU3WuvKWveLGGqPn2fWY2LKQdzGy7HMd8Bfg34INmtj1wcMhXlv3Tr7eRqOk4ZVfgb3ncg3M1w4Oqc+XtZWB0WFsXM9sC/BC4XNIoAEnvlXREjnMMJQrEr0vakWiZvfRrZByTGmqy84BvSBoqaTfgPODmPtyTc1XLg6pz5a0NeAJ4SdKakPc14GngIUn/AOYT1USzuQIYTFTrfAj4bdr27wInhd6738tw/BeB9cAK4H7gFuD63t2Oc9XN11N1zjnnCsRrqs4551yBeFB1zjnnCsSDqnPOOVcgHlSdc865AvGg6pxzzhWIB1XnnHOuQDyoOueccwXiQdU555wrEA+qzjnnXIH8f5NjuX4V7nBpAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig, axs = plt.subplots(1,1)\n", "axs.plot(energy, 'kx', label=\"energy\")\n", "axs.set_ylabel(\"Energy (Hartree)\", fontsize=12)\n", "axs2 = axs.twinx()\n", "axs2.plot(distance, 'r+', label=\"Distance\")\n", "axs2.set_ylabel(\"Distance (Bohr)\", fontsize=12)\n", "axs.set_xlabel(\"Iteration\", fontsize=12)\n", "axs.legend(loc=\"upper center\")\n", "axs2.legend(loc=\"upper right\")\n", "axs.ticklabel_format(useOffset=False)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.10" } }, "nbformat": 4, "nbformat_minor": 4 }