{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Running a wavelet computation on a methane molecule, with AiiDa\n",
    "\n",
    "This is a copy of the CH4 tutorial, please check the CH4 tutorial first for details about computations and datasets.\n",
    "The purpose of this tutorial is to outline the integration of AiiDa with BigDFT workflows and datasets. This allows jobs to be offloaded to remote supercomputers and database storage of computed results, as well as transparent integration with PyBigDFT analysis tools.\n",
    "\n",
    "## Introduction: AiiDa installation\n",
    "\n",
    "AiiDa must be functional on your system for this tutorial. Please install the prerequisites from [here](https://aiida.readthedocs.io/projects/aiida-core/en/latest/install/prerequisites.html), follow the steps from the guide [here](https://aiida-core.readthedocs.io/en/latest/install/installation.html) to setup the latest version of AiiDA, and take some time to familiarize with this environment and verdi command line tools.\n",
    "\n",
    "This tutorial should be run from Aiida's python virtual environment. Jupyter should also be setup in a way that either \"import aiida\", \"%aiida\", or with a compatible aiida-aware kernel. Check [this page](https://aiida.readthedocs.io/projects/aiida-core/en/latest/intro/installation.html#intro-install-jupyter) for details.\n",
    "\n",
    "Once aiida is installed, you need to install the BigDFT plugin, with 'pip install aiida-bigdft', or locally, using the one available on https://github.com/aiidateam/aiida-bigdft-plugin.\n",
    "\n",
    "Once installed, you need to start the verdi daemon with PyBigDFT's path in the PYTHONPATH. If BigDFT is installed locally, this can be done by executing 'source BIGDFT_INSTALL_PATH/bin/bigdftvars.sh' before running. \n",
    "\n",
    "If BigDFT is not installed locally, you may need to execute 'export PYTHONPATH=BIGDFT_SOURCE_PATH/PyBigDFT/:BIGDFT_SOURCE_PATH/futile/src/python'\n",
    "\n",
    "Next you need to setup machines and codes you will use.\n",
    "In this tutorial we will use a 'localhost' machine, without scheduler. BigDFT needs to be installed on the machine."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aiida setup for the tutorial\n",
    "\n",
    "#### Computer\n",
    "`verdi computer setup`\n",
    "\n",
    "- Computer label: localhost \n",
    "- Hostname: localhost \n",
    "- Description []:  \n",
    "- Transport plugin: local \n",
    "- Scheduler plugin: direct \n",
    "- Shebang line (first line of each script, starting with #!) [#!/bin/bash]:  \n",
    "- Work directory on the computer [/scratch/{username}/aiida/]: path/to/a/folder/to/store/aiida/generated/files \n",
    "- Mpirun command [mpirun -np {tot_num_mpiprocs}]:  \n",
    "- Default number of CPUs per machine: number of cores of your computer (or less)\n",
    " \n",
    "Then, no need for a specific script, just hit escape and then :wq to exit.\n",
    "              \n",
    "Once this is done, just launch\n",
    "\n",
    "`verdi computer configure local localhost`\n",
    "\n",
    "and keep the default setting to allow usage of the computer. \n",
    "\n",
    "\n",
    "#### Code\n",
    "\n",
    "`verdi code setup`\n",
    "\n",
    "- Label: bigdft\n",
    "- Description []: bigdft on localhost\n",
    "- Default calculation input plugin: bigdft\n",
    "- Installed on target computer? [True]:    \n",
    "- Computer: localhost\n",
    "- Remote absolute path: /path/to/bigdft/install/bin/bigdft\n",
    "\n",
    "And in the \"pre execution script\" section, add \n",
    "\n",
    "`source /path/to/bigdft/install/bin/bigdftvars.sh`\n",
    "\n",
    "Be careful, an incorrect path here necessitates a re-creation of the code.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Running the code\n",
    "If not already activated, activate the aiida env with\n",
    "\n",
    "`source ~/envs/aiida/bin/activate`\n",
    "\n",
    "Check verdi status with `verdi status`\n",
    "\n",
    "Start postgres service with `sudo service postgresql start`\n",
    "\n",
    "Start rabbitmq server with `sudo service rabbitmq-server start`\n",
    "\n",
    "And verdi daemon with `verdi daemon start 2`\n",
    "\n",
    "Be careful not to start more processes than there are physical processors on your machine.\n",
    "\n",
    "These steps will likely be neccessary after each reboot,\n",
    "consider adding to an `aiida_service_start.sh` file for easier startup."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#prepare computation the same way as in normal computation:\n",
    "inp = { 'dft': { 'hgrids': 0.55, 'rmult': [3.5, 9.0]} }\n",
    "# becomes\n",
    "from BigDFT import Inputfiles as I\n",
    "inp=I.Inputfile()\n",
    "inp.set_hgrid(0.55)\n",
    "inp.set_rmult([3.5,9.0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "inp.set_atomic_positions('CH4_posinp.xyz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A simple calculation can now be run on the selected platform using an AiidaCalculator instead of SystemCalculator.\n",
    "- code is a mandatory option, to select the code and the target computer\n",
    "- num_machines, mpiprocs_per_machine, omp (number of OpenMP threads per MPI process), walltime are optional parameters "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initialize an Aiida Calculator for Remote code 'bigdft' on localhost, pk: 256, uuid: 3f301d19-e59e-4b8f-92b6-74cd7ad689d8 with 1 machine(s) , 1 processes per machine, and 1 cores per process\n",
      "Creating the yaml input file \"./input.yaml\"\n",
      "setting data dir to /home/pi/.aiida/repository/pi/repository/node/14/51/1d86-1638-415d-ae9a-e01f77b35009/path\n"
     ]
    }
   ],
   "source": [
    "%aiida\n",
    "from BigDFT import AiidaCalculator as calc\n",
    "#reload(calc)\n",
    "study = calc.AiidaCalculator(code=\"bigdft@localhost\", \n",
    "                             num_machines=1, \n",
    "                             mpiprocs_per_machine=1, \n",
    "                             omp=1,\n",
    "                             walltime=3600) #Create a calculator\n",
    "ch4 = study.run(input=inp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*ch4* is the instance of the class Logfile which can handle easily all information coming from the output file **log.yaml**.\n",
    "Then we can display some information as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(ch4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Construction of the input dataset dictionaries\n",
    "As in plain BigDFT notebook, you can group computations in Datasets.\n",
    "Aiida's main feature in this case is the possibility to launch all the computations from the dataset asynchronously, creating and tracking multiple parallel jobs to execute on the target computer.\n",
    "\n",
    "metadata can be assigned to Datasets to tune runtime parameters or scheduler options for all internal computation tasks.\n",
    "\n",
    "Let us build two different dataset of inputfiles, the first one synchronous, the other asynchronous.\n",
    "We also label each of the input dictionaries by a unique name identifying the run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from BigDFT import Datasets as D\n",
    "import copy\n",
    "hgrids = [ 0.55 - i*0.05 for i in range(8)]\n",
    "crmult = [ 3.5 + i*0.5 for i in range(8)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "metadata={\n",
    "    #'queue_name':'my_queue', \n",
    "    #'account':'my_account', \n",
    "    #'qos':'myqos'\n",
    "}\n",
    "h_only_dataset=D.Dataset('h_only',run_dir=\".\",metadata=metadata)\n",
    "for h in hgrids:\n",
    "    inp_run=copy.deepcopy(inp)\n",
    "    inp_run['dft']['hgrids'] = h\n",
    "    label_run='h='+str(h)\n",
    "    h_only_dataset.append_run(id={'h':h} , runner=study , input=inp_run)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "c_only_dataset=D.Dataset('c_only',run_dir=\".\", asynchronous=True,metadata=metadata)\n",
    "for c in crmult:\n",
    "    inp_run=copy.deepcopy(inp)\n",
    "    inp_run['dft']['rmult'] = [c,9.0] #only change crmult\n",
    "    label_run='c='+str(c)\n",
    "    c_only_dataset.append_run(id={'c':c} , runner=study , input=inp_run)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "h_and_c_dataset=D.Dataset('h_and_c',run_dir=\".\", asynchronous=True,metadata=metadata)\n",
    "for h,c in zip(hgrids,crmult):\n",
    "    inp_run=copy.deepcopy(inp)\n",
    "    inp_run['dft']['hgrids'] = h\n",
    "    inp_run['dft']['rmult'] = [c,9.0] #only change crmult\n",
    "    h_and_c_dataset.append_run(id={'h':h,'c':c} , runner=study , input=inp_run)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then run the first dataset synchronously, one job at a time, which will block us a while:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h_only_dataset.run()\n",
    "h_only_results = h_only_dataset.results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Launch  the second and third datasets asynchronously.\n",
    "\n",
    "- WARNING : if you are running locally, this will launch several instances of bigdft at the same time, which could use a lot of memory. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c_only_dataset.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h_and_c_dataset.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As the dataset has been launched asynchronously, we now must wait for the end of execution, and gather the results (we can do other work while the dataset is running).\n",
    "The wait() call is an optional commodity which will return only when AiiDa reports that no process from the dataset is running anymore, allowing to postprocess the data.\n",
    "- Hint : use jupyter-notify to get a notification in your browser when a cell finishes : pip install jupyternotify, then uncomment next cell and the %notify magic which activates it for a given cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%load_ext jupyternotify"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'0 processes still running'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#%notify\n",
    "c_only_dataset.wait()\n",
    "h_and_c_dataset.wait()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Aiida retrieves logfiles and other output files from BigDFT automatically. To load them in this notebook, get_logfiles() can be used"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c_only_results=c_only_dataset.get_logfiles()\n",
    "h_and_c_results=h_and_c_dataset.get_logfiles()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now store the energies of each of the dataset runs, and identify the minimum as the minimum value from the `h_and_c` dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "energies_h=np.array([h_only_results[c].energy for c in h_only_results])\n",
    "energies_c=np.array([c_only_results[c].energy for c in c_only_results])\n",
    "energies_hc=np.array([h_and_c_results[c].energy for c in h_and_c_results])\n",
    "#find the minimum\n",
    "emin=np.min(energies_hc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We plot the energy values varying the grid spacing or the extension"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fe018d0d710>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEKCAYAAAAb7IIBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi40LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcv7US4rQAAIABJREFUeJzt3Xl4VdXZ9/HvnRASIYAMigxiUJmJgERQUaFgFRWVWBAoioiK+mi1UqeqfaWtLVR87GPrUKkDYlVUKAiCs0yKA0FQGUQZJaAyzwIZ1vvHPgknISHDOcneyfl9riuX5+yzh5ttsu6z19r7XuacQ0REYk+c3wGIiIg/lABERGKUEoCISIxSAhARiVFKACIiMUoJQEQkRikBiIjEKCUAEZEYpQQgIhKjavgdwNE0atTIpaSk+B2GiEiVsmjRoq3OueNKWi/QCSAlJYWMjAy/wxARqVLMbH1p1qvUBGBm/YFLgLrAs865dyvz+CIiclipxwDM7Dkz22xmSwst72tmK81slZnde7R9OOemOeduAG4CBpUvZBERiYayXAFMAB4HJuYtMLN44Angl0AmsNDMpgPxwJhC249wzm0OvX4gtJ2IiPik1AnAOTfPzFIKLe4GrHLOrQEws0nA5c65MUC/wvswMwPGAm85574o6jhmNhIYCdCiRYvShiciImUU6W2gzYANYe8zQ8uK8xvgfGCAmd1U1ArOufHOuTTnXNpxx5U4iC0iIuVUqYPAzrl/AP+ozGOKiEjRIr0C2AicGPa+eWhZRMzsUjMbv2vXrkh3JSIixYg0ASwEWplZSzOrCQwGpkcalHNuhnNuZL169SLdlYhIleGcY9nq9Xz+5PXs3r2zwo9XlttAXwE+AdqYWaaZXeecywZuBd4BVgCvOeeWVUyoIiLV0/fb9vPPD75jxLgXqfXCL+n8039Zv2RuhR+3LHcBDSlm+SxgVtQiwusCAi499dRTo7lbEZHA2L7vEDO/2sS0JZtYtH4H58ct4snEJyGpFocGziC1VY8KjyGQpSCcczOAGWlpaTf4HYuISLT8fCiH91f8xLTFG5n77Raycx2tj6/Nq23n033dU9CkCwx6Ceod7WbK6AlkAhARqS5ych0LVm9l2uJNvL30B/YdyuGEuklcd05L+nc4lraf3YMtfwNOGwSXPgYJx1RabIFMAOoCEpGqzDnHsk27mbp4IzO+3MTmPQepk1iDS05rQv/Ozeh+ckPid62HSQNg83K44CE461Ywq9Q4A5kA1AUkIlXRhu37eWPJRqYu3sjqLftIiDd6tTme9C7N6N32eJIS4r0V186D164BlwNDX4dTz/cl3kAmABGRqmLHvkPM/PoHpi3eSMb6HQB0S2nAiHNacklqE46tVfPwys7B5+Ph7d9Dw1NhyCvQ8BSfIlcCEBEpswNZeYO5m5j77Waychytjk/mrgvbcFmnppzYoNaRG2UfhJmjYPF/oM3FkP40JNWt/ODDBDIBaAxARIImJ9fx6ZptTF28kbeX/sjeg9k0rpvI8LNT6N+lGe2b1MWK68Pf8yO8ejVkfg7n3Q29fg9x/s/IG8gEoDEAEQmCvMHcN5ZsZPqXm/hp90GSE2vQt+MJpHdpxpknNyQ+roSB242LYNJVcGAnDHwBOvSvnOBLIZAJQETET5k79vPGkk1MW7yR7zbvpUacN5j7h35NOb9d48ODuSX5chJMvw3qNIbr3oUTUis28DJSAhARAXbu9wZz31i8ic/XbQcg7aT6PNS/I5ekNqF+7Zol7CFMTja8/yB88jiknOt986/dsIIiLz8lABGJWQeycvjwm81MXbyROSu9wdxTjqvNnRe05vLOzYoezC3Jzztg8ghY/SF0uxEu/AvEJ0Q/+CgIZALQILCIVJTcXMena7cxbfFG3vr6R/YczOa4OokMOyuF9C7N6ND0KIO5Jdm8Al4ZArsy4bJ/wunDoht8lAUyAWgQWESiyTnHih/28MaSjbyxZBM/7j5A7Zrx9O3YhP5dmnL2KY1KHswtyTcz4b8jIaEWDJ8JLbpHJ/gKFMgEICJyNLm5jt0Hsti+7xA79h9i+74sduw/xI59h9ie99+wZdv2HWLXz1nUiDN6tj6O+y9px/ntGnNMzVIO5h49GJj/CMz+CzSt3GJukVICEBFfOefYfSCbnfsPFWzQCzTmh9i5Pyv//Y79h8h1Re+vZnwcDWrXpH7tmtSvlUC7pnVpUKsmrU+owyWpTWhQlsHckhzcC9NuhhXTfSnmFiklABGJGucc+w7l5DfaxTXg4Q39zv2HyC6mNU+IN46tVZMGtWpSv3YCrRsnU79WTa+BDy0Lf9+gdk1q1Ywvfx9+WexYB6/8GrasgAv+AmfdUunF3CIVyASgQWARfzjn+Dkrh70Hstl7MPRzIJs9of/mLdtzIJs9B7K8Rj3UmHvdLVkcysktct/xcUb9WgmhhrsmJzdKputJ3rf08Aa8fu3DDX5yYo3KaczLas1ceH14qJjbZDi1j98RlYs5V8x1VACkpaW5jIwMv8MQCbycXFegwd57MIs9B8LfF3xdoEHPW3Ygi70Hs4vtWglXMz6O5KQa+Y334W/pNWkQ/q08r2GvVZM6STWIi3Sg1W/hxdwatYLBL/tazK04ZrbIOZdW0nqBvAKI1L/nrWH+qq3EG8SZERdnxJsRF+e9j48zb7kZ8XGF1jHC1i+0Tmjb+DjDDOIt77URbxx+HdreQssOr+Pt6/Drw+sY3msDMDC8WMwOL7fQh4ffG4YXW96XJAtta1bwdVwx24avA96/vfDy8Jgs75zmvQ+dm/A4Dq8Ttq8gfosrJ+ccuQ5ynSMn1+Ec5DhHrnPk5nqfectdaDmh5d76edsW2D70+cHs3CIa6awC38L3HTryW/n+Qzmlir12zXiSk2qQnBj6SapBo+RaJCcmUCfp8LLkxBqH34eW1UlMIDmpBrUT40msEYXB06qmQDG3S+CKpyGxjt9RRaRaJoCfs3LY9XOW9weYW+iPNO8P0Tlycwn7ozz8h5v3hxz+x5vjvH1I+eQlhZISRV4C8tY5nLwOr+Ml8rwEGRfKUHlJKS6UaOLC1itvQ11wffJ/byr79yA+zvIb4rxGuX7tmpzYoBZ1kmpQu2bhRjvhyEY8tF7EtzrGqj0/wqtXQebCQBVzi1S1TAC39WnFbX1aRX2/rnCSCP9Gl58wSkgsecvDGpTcXIfDu7p0rpjXuNB/w5eHlh2xHMA7blHbElon1x25T0LrF7UtoXXy4s11h2PIzfssP6bwb7retrklrBf+/nAMh/efH1PoPObFkBsKMO/8Fl4v/AosrtDVXN6VXN5VYFzoiixvnbz3VuhKMO/KMO+KL84K7jfvSjLOCl4NFjiOWcHt44zEUNdK+LfupIS4anUFVeVkLoJXh8KBXYEr5hapapkAKkp4V4+IxIAlr8CM20PF3N6DEzr6HVFUKQGIiBRWRYq5RUoJQEQk3P7tXjG3NbMDX8wtUoFMAHoOQER8kVfMbfdGuOxxOP1qvyOqUIEcxnbOzXDOjaxXr57foYhIrPhmJjxzPmTt94q5VfPGHwJ6BSAiUmlyc2HeOJjzV2h6Ogx+Ceo29TuqSqEEICKx6+BemHYTrJgBpw0OFXNL8juqSqMEICKxaftamDTUK+Z24V/hzP+pcsXcIqUEICKxZ81ceP0a74nBq6bAKb39jsgXgRwEFhGpEM7BZ0/Di+mQ3Bhu+DBmG3/QFYCIxIrsg/DmKFhSfYq5RUoJQESqv90/wGtXe8Xcet4DPe+tFsXcIqUEICLV1887vXIOnz7ldf9cORHaX+53VIERyASgJ4FFJCIHdsNn/4IFj8PBXdC+P/R+wJvERfIFMgE452YAM9LS0m7wOxYRqUIO7fNm7Pr4Mfh5h9fX/4vfwwmpfkcWSIFMACIiZZL1M2Q8Bx/9HfZtgVN/Cb+4D5qd7ndkgaYEICJVV/ZB+GIizHsE9v4ILXt6XT0ndvM7sipBCUBEqp6cLG9u3nmPwO5MaHE2DHgWUs7xO7IqRQlARKqOnGz46lWY+zfYuR6apcHl/4STfxFzZRyiQQlARIIvNweW/hfmjoVtq6BJJ7h4HLS6QA1/BJQARCS4cnNhxXSYMwa2fAPHd4BBL0HbS9TwR4ESgIgEj3Ow8i2Y/Vf46Wto1BoGPO/dz68neKNGCUBEgsM5WPUBzP4LbPoCGpwM6eMhdQDExfsdXbWjBCAiwbBmrtfwb/gM6rXw5uTtNATi1UxVFJ1ZEfHX+k+8hn/dfKjTFC55FLpcDTVq+h1ZtacEICL+yMzwGv7VH0Lt46Hv36Dr8JiaktFvlZYAzKwdcDvQCPjAOfdUZR1bRAJk0xJvcPe7d6BWQ7jgIUi7DmrW8juymFOqBGBmzwH9gM3OuY5hy/sCjwHxwDPOubHF7cM5twK4yczigImAEoBILPlpmdfwf/MmJB0Lff4fdLsREpP9jixmlfYKYALwOF7DDYCZxQNPAL8EMoGFZjYdLxmMKbT9COfcZjO7DLgZeDHCuEWkqtjyrXcf/7Kp3gxcvX4PZ94MSfX8jizmlSoBOOfmmVlKocXdgFXOuTUAZjYJuNw5NwbvaqGo/UwHppvZTODl8gYtIlXAttUw92H4+jWocQycOwrOuhVqNfA7MgmJZAygGbAh7H0m0L24lc2sF3AFkAjMOsp6I4GRAC1atIggPBHxxc7vvYZ/ycsQXxPOugV6/BZqN/I7Mimk0gaBnXNzgDmlWG88MB4gLS3NVWxUIhI1uzd51Tm/mAgWB91ugHNGQZ3GfkcmxYgkAWwETgx73zy0TERiyZ6fvIlYMp4DlwunD4Nzfwf1mvkdmZQgkgSwEGhlZi3xGv7BwK+jEZTmBBapAvZthY//Dz5/BnIOQedfw3l3Qf2T/I5MSqm0t4G+AvQCGplZJvCgc+5ZM7sVeAfvzp/nnHPLohGU5gQWCbivJ8OM2yFrP6ReCT3vhoan+B2VlFFp7wIaUszyWRxlQFdEqpmcLHjvQfj0CWhxFlz6GBzXxu+opJwCWQpCXUAiAbR3M7x+Laz/CLrf5D3BG5/gd1QSgUAW1nbOzXDOjaxXTw+KiARCZgY83RM2LvLKM1/0NzX+1UAgrwBEJEAWTYBZd0GdE+C6d6HJaX5HJFESyASgLiCRAMg+6DX8X7wAp/SGXz2rp3irGXUBiciRdm2E5y/yGv9zfwdDJ6vxr4YCeQUgIj5aOx9eH+5dAQz6D7S71O+IpIIoAYiIxzn49El49w/ePf2DXoLjWvsdlVSgQCYAjQGIVLJD+2D6bbB0MrTtB/2fgqS6fkclFUxjACKxbvsaeOaXsHSKN0nLlS+q8Y8RgbwCEJFK8t17MOU6wOCqyXDq+X5HJJVICUAkFuXmwvxHvCkaT+joDfbWT/E7KqlkSgAisebALph6E6ycBacNgn7/pwnZY1QgE4AGgUUqyOYVMGko7FwPFz0M3UaCmd9RiU80CCwSK5ZNhX/3gYN74JoZ0P1GNf4xLpBXACISRTnZ8OGf4OPHoPkZcOVEqNvU76gkAJQARKqzfdtgyghYMwfSRkDfsVAj0e+oJCCUAESqq02L4dWrvTr+lz0Op1/td0QSMIEcAzCzS81s/K5du/wORaRqWvIyPHuhV95hxNtq/KVIgUwAGgQWKafsQzDzdzDtZmjRHW6cC81O9zsqCSh1AYlUF7t/gNeGQebncPZvoM9oiNefuBRPvx0i1cH6T+D1a+DgXhjwPHS8wu+IpApQAhCpypyDhc/A2/fCsS3g6mnQuL3fUUkVoQQgUlVl/Qxv3gFfvgKt+0L603DMsX5HJVWIEoBIVbRjPbx6Ffz4FfT6PZx3N8QF8p4OCbBAJgDVAhI5itUfwuQRXkXPIa9Cm75+RyRVVCC/Mug2UJEiOAfzH4X//ArqNIGRs9X4S0QCeQUgIoUc3OPd279iBnS4Ai77JyQm+x2VVHFKACJBt/U7r4TztlVwwV/grFtUxVOiQglAJMhWvOlN3lIjEYZNg5bn+R2RVCNKACJBlJvjTdc4/xFo2sWbsrFec7+jkmpGCUAkaPZvh//eAKvehy5XwcX/CwlJfkcl1ZASgEiQ/Pi119+/e5M3V2/X4ervlwqjBCASBLk58MUL8PZ93tO8174FJ57hd1RSzSkBiPjJOfj2HXh/NGxZASnnwoDnIPl4vyOTGBDIBKAngSUmZC6C9/4A6z+GBifDwBeg/eXq8pFKE8gE4JybAcxIS0u7we9YRKJu22r44E+wfBrUagQXP+L19ccn+B2ZxJhAJgCRamnvFpj3MGQ8B/E1oec93sQtiXX8jkxilBKASEU7tA8+eRI+fgyy9sPpw6DXvVDnBL8jkxinBCBSUXKyYcl/YPYY2PsjtO0HfR6E41r7HZkIoAQgEn3Owcq3vDt7tq6E5t3gyhegxZl+RyZSgBKASDRlZsC7f4DvF0DDU70SDm376c4eCSQlAJFo2LYaPvgjLH8Dah8Plzzq9fXrzh4JMCUAkUjs3QJz/waLnof4RG96xrNuVa1+qRKUAETK4+Be+OQJWPAPb3L2rsO92zrrNPY7MpFSUwIQKYucbFg8EeaMhb0/QbtLvTt7GrXyOzKRMlMCECkN52DlrNCdPd/CiWfClS9Ci+5+RyZSbkoAIiXZ8Ll3Z8+GT6FhKxj8MrS5WHf2SJWnBCBSnK2r4IPR3kTsyY29+vxdroZ4/dlI9VCpv8lmVhuYC4x2zr1ZmccWKbW9m70+/kUTIOEY6HWfNxG77uyRaqZUCcDMngP6AZudcx3DlvcFHgPigWecc2NL2NU9wGvljFWkYh3cC588Dh//A3IOQtoI6Hm3avNLtVXaK4AJwOPAxLwFZhYPPAH8EsgEFprZdLxkMKbQ9iOATsByQJObSrDkZMEXoTt79m32avL3eRAanuJ3ZCIVqlQJwDk3z8xSCi3uBqxyzq0BMLNJwOXOuTF4VwsFmFkvoDbQHvjZzGY553LLH7pIhJyDb9707uzZtgpanO0N8GoqRokRkYwBNAM2hL3PBIq9J845dz+AmQ0HthbX+JvZSGAkQIsWLSIIT+Qovv8U3vt/sOEzaNQGBr8CbS7SnT0SUyr9dgbn3IQSPh8PjAdIS0tzlRGTxJAt33o1e755E5JPgEsfg85X6c4eiUmR/NZvBE4Me988tCximhNYom7PTzBnjNfXn1ALej8AZ/4P1Kztd2QivokkASwEWplZS7yGfzDw62gEpTmBJWoO7oEF/4QFj3t39pxxvXdnT+1Gfkcm4rvS3gb6CtALaGRmmcCDzrlnzexW4B28O3+ec84tq7BIRcpq7TyYfJ13Z0+HdOj9B93ZIxKmtHcBDSlm+SxgVlQjQl1AEiHnYOEz8NY93qQsQyZB865+RyUSOHF+B1AU59wM59zIevXq+R2KVDXZh+DNO2DWnXDq+XD9+2r8RYqhWx+k+ti3FV4bBus/hnPu8Lp84uL9jkoksJQApHr4cSm8MsTr77/iGThtoN8RiQReILuAzOxSMxu/a9cuv0ORqmDFDHj2AsjNgmtnqfEXKaVAJgCNAUipOAdzH4ZXr4Lj28HIOdBM/f0ipaUuIKmaDu2Daf8Dy6fBaYO9J3oTVGdQpCyUAKTq2bkBJg2Bn5bBBQ/BWbeqho9IOQQyAeg5ACnW+k+8Lp+cQ/Dr16DVL/2OSKTK0hiAVB1fTIQXLoWkenD9B2r8RSIUyCsAkQJysuHd++Gzf8EpvWHAc3BMfb+jEqnylAAk2PZvh8nXwpo5Xl//+X9U6WaRKAnkX5LGAASAzd/AK4Nh90a4/AnocpXfEYlUKxoDkGBa+TY8c753u+c1b6rxF6kAgUwAEsOcg4/+7n3zb3gyjJwNLYqdaVREIhDILiCJUVk/w/TfwNevQ4crvG6fmrX8jkqk2lICkGDYvQkm/Ro2Lfamazz3Tj3cJVLBlADEf5kZMGkoHNoLg1+Gtpf4HZFITAjkGICqgcaQL1+F5y+GGolw3btq/EUqUSATgO4CigG5OfDuH2DqSDixG9wwGxp38DsqkZiiLiCpfAd2wZTr4bt34Yzroe9YiE/wOyqRmKMEIJVr22rvFs/ta+CSR+GM6/yOSCRmKQFI5Vn9Ibw+HCwehr0BKef4HZFITAvkGIBUM87Bp0/Bf34FdZt7D3ep8Rfxna4ApGJlH4SZo2Dxf6BtP0h/GhKT/Y5KRAhoAlAxuGpi72Zv8pYNn0HPe6DnvRCni06RoAjkX6NuA60GNi2B8b3gh69g4AT4xX1q/EUCJpBXAFLFLf2vN2F7rYZw3TvQpJPfEYlIEZQAJHpyc2HOX2HeODjxTBj0IiQf73dUIlIMJQCJjoN74L83wsqZ0OVquOR/vfIOIhJYSgASuR3r4JUhsGUlXPQwdBupSp4iVYASgERm7Xx4bRi4HLhqCpzyC78jEpFS0m0ZUn4Ln4EX+0PtRl4xNzX+IlWKrgCk7HKy4K27IeM5aHUB/OoZSNItuyJVjRKAlM2+bV6Xz/qPoMft0OdBiIv3OyoRKYdAJgA9CRxQ+7fDhEu8Sp7p46HTIL8jEpEIBHIMQE8CB9DBPV4xt+1r4KrJavxFqoFAXgFIwGT97N3m+cOXMPglaHme3xGJSBQoAcjR5WTB69fCuo/gin9Dm4v8jkhEokQJQIqXmwvTboZv3/Ke7D1toN8RiUgUBXIMQALAOZh1J3z9unenzxnX+x2RiESZEoAU7YM/Qcaz0OO3cO4ov6MRkQqgBCBH+ujv8NGj0PVaOH+039GISAVRApCCMp6D90dDxwFev7+KuolUW0oActjXk+HNUdDqQkj/l57wFanmlADEs/JtmHojnNQDrnwB4hP8jkhEKpgSgHglnV+/Bk5IhSGvQMIxfkckIpVACSDWbVwErwyG+ikwdAok1fU7IhGpJEoAsWzzCq++T62GcPVUqN3Q74hEpBJVWgIws15mNt/M/mVmvSrruFKM7WthYn+IT4Rhb0Ddpn5HJCKVrFQJwMyeM7PNZra00PK+ZrbSzFaZ2b0l7MYBe4EkILN84UpU7P4BJl4OOQe9b/4NWvodkYj4oLS1gCYAjwMT8xaYWTzwBPBLvAZ9oZlNB+KBMYW2HwHMd87NNbPGwKPA0MhCl3LZv92bxnH/Nhg2HRq39zsiEfFJqRKAc26emaUUWtwNWOWcWwNgZpOAy51zY4B+R9ndDiCx7KFKxA7shv9c4XX/XDUZmnf1OyIR8VEk1UCbARvC3mcC3Ytb2cyuAC4EjsW7mihuvZHASIAWLVpEEJ4UkF/T/yvV9BcRoBLLQTvn/gv8txTrjQfGA6SlpbmKjism5GTB68Nh/ceq6e+zrKwsMjMzOXDggN+hSDWQlJRE8+bNSUgo34ObkSSAjcCJYe+bh5ZFTHMCR1FuDky9Cb59WzX9AyAzM5M6deqQkpKCqc6SRMA5x7Zt28jMzKRly/LdyBHJbaALgVZm1tLMagKDgekR7C+f5gSOkrya/ksnq6Z/QBw4cICGDRuq8ZeImRkNGzaM6GqytLeBvgJ8ArQxs0wzu845lw3cCrwDrABec84tK3ckEn0f/NGr7qma/oGixl+iJdLfpVIlAOfcEOdcE+dcgnOuuXPu2dDyWc651s65U5xzf4kokjBmdqmZjd+1a1e0dhl75j/q1fVXTX+pZBMmTODWW28FYNq0aSxfvrzU237++ed07tyZzp0706lTJ6ZOnVrkesOHD6dly5b56y5ZsiQqsceaQJaCUBdQhBY+6337V01/KaPs7Oyo7q+sCaBjx45kZGSwZMkS3n77bW688cZiYxo3bhxLlixhyZIldO7cOVohxxRNCl/dfPU6zPwdtO6rmv5SpIkTJ/LII49gZpx22mnEx8eTlJTE4sWL6dGjB3Xr1mXt2rWsWbOG77//nr///e98+umnvPXWWzRr1owZM2aQkJBASkoKGRkZNGrUiIyMDO68807mzJmTf5wFCxYwffp05s6dy0MPPcSUKVM45ZRTjhpbrVq18l8fOHBA3WUVLJAJQHcBldPKt7ya/innwMAJqukfcH+csYzlm3ZHdZ/tm9blwUs7FPv5smXLeOihh1iwYAGNGjVi+/btjBo1iszMTBYsWEB8fDyjR49m9erVzJ49m+XLl3PWWWcxZcoUHn74YdLT05k5cyb9+/cvMZazzz6byy67jH79+jFgwADA+9b+0ksvHbHueeedxz/+8Q8APvvsM0aMGMH69et58cUXqVGj6Gbq/vvv509/+hN9+vRh7NixJCbq+dKyUhdQdbF2Prx2DTTppJr+UqwPP/yQgQMH0qhRIwAaNGgAwMCBA4mPP3y1eNFFF5GQkEBqaio5OTn07dsXgNTUVNatW1fu499111353TbhP3mNP0D37t1ZtmwZCxcuZMyYMUXe5TJmzBi++eYbFi5cyPbt2/nb3/5W7phiWSCvAKSM8mr6N2gJV02BxDp+RySlcLRv6pWtdu3aBd7nfZuOi4sjISEhvysmLi4uv0++Ro0a5ObmApT6VsTSXAHkadeuHcnJySxdupS0tLQCnzVp0iQ/zmuvvZZHHnmkVMeXggJ5BSBlUKCm/zSo1cDviCTAevfuzeuvv862bdsA2L59e7n3lZKSwqJFiwCYMmVKkevUqVOHPXv25L8v6Qpg7dq1+Qlm/fr1fPPNN6SkpByx3x9++AHwHoaaNm0aHTt2LPe/I5YFMgHoNtBSOqKmfxO/I5KA69ChA/fffz89e/akU6dOjBpV/udDHnzwQW6//XbS0tIKdB+FGzx4MOPGjaNLly6sXr26xH1+9NFHdOrUic6dO5Oens6TTz6Z31118cUXs2nTJgCGDh1KamoqqampbN26lQceeKDc/45YZs4Ft9xOWlqay8jI8DuMYNr9Azx3IRzcDde+Bce38zsiKYUVK1bQrp3+X0n0FPU7ZWaLnHNpxWyST2MAVVF4Tf9rpqvxF5FyUQKoagrU9J8CzVTTX0TKRwmgKsmr6f/j1zDoJWh5rt8RiUgVpkHgqiK8pn/609Cmr98RiUgVF8gEoAfBCsnN8Z7wzavpnzrA74hEpBoIZAKQMM55tX2WTvGqep5xnd8RiUg1oQQQdO+PhkXPwzl3eD8iAXX22WdHZT/h5aTLa/o6DtUAAAAPDklEQVT06YwdO7bIz5KTkyPad2VJSUlh69at7Ny5kyeffLJCjqEEEGTzH4WP/w/SRngzeon4yDmXX/qhKAsWLKjEaIqXnZ3NZZddxr333hvRftatW0evXr3KvF1J56msYi4BaBAYWPiMV9M/dSBcrJr+Eh333nsvTzzxRP770aNH88gjj7B371769OnD6aefTmpqKm+88QbgNYJt2rRh2LBhdOzYkT//+c/89re/zd/+3//+N3fc4V2Z5n2znjNnDr169WLAgAG0bduWoUOHkvfA6axZs2jbti1du3bltttuo1+/fkXGuWnTJvr27UurVq24++6785c/++yztG7dmm7dunHDDTfkXykMHz6cm266ie7du3P33XcXuIpYu3YtZ511FqmpqQWeGP7hhx8477zz6Ny5Mx07dmT+/PmlOoc//fQT6enpdOrUiU6dOrFgwYIjztOGDRtITk7mrrvuokOHDpx//vl8/vnn9OrVi5NPPpnp073Zcwtf7fTr169ASe28/2erV6+mc+fO3HXXXaWKsbQCeRuoc24GMCMtLe0Gv2PxxVevwcw7vZr+/Z+CuEDmaYnUW/d6t/RG0wmpcFHRXR8AgwYN4re//S233HILAK+99hrvvPMOSUlJTJ06lbp167J161bOPPNMLrvsMgC+++47XnjhBc4880z27t1Lp06dGDduHAkJCTz//PM8/fTTRxxn8eLFLFu2jKZNm9KjRw8+/vhj0tLSuPHGG5k3bx4tW7ZkyJAhxca5ZMkSFi9eTGJiIm3atOE3v/kN8fHx/PnPf+aLL76gTp069O7dm06dOuVvE17SesKECfnLb7/9dm6++WaGDRtWIPm9/PLLXHjhhdx///3k5OSwf//+Up3i2267jZ49ezJ16lRycnLYu3cvO3bsKHCeAPbt20fv3r0ZN24c6enpPPDAA7z33nssX76ca665Jv/8lmTs2LEsXbq0QmY9C2QCiGkr34KpN6mmv1SILl26sHnzZjZt2sSWLVuoX78+J554IllZWdx3333MmzePuLg4Nm7cyE8//QTASSedlN+oJScn07t3b958803atWtHVlYWqampRxynW7duNG/eHIDOnTuzbt06kpOTOfnkk2nZsiUAQ4YMYfz48UXG2adPH/LuAmzfvj3r169n69at9OzZs0AJ62+//TZ/m8IlrfN8/PHH+cXqrr76au655x4AzjjjDEaMGEFWVhb9+/fPn1UsPT2dtWvXcujQIb7//vv85bfffjvXXnstH374IRMnTgQgPj6eevXqsWPHjgLnCaBmzZoFymgnJibml9iOpKR2NCkBBMnaearpH0uO8k29Ig0cOJDJkyfz448/MmjQIABeeukltmzZwqJFi/Jn+8or8Vy4VPT111/PX//6V9q2bcu1115b5DHCJ2eJj48v81ST5dm+cJzhippZ7LzzzmPevHnMnDmT4cOHM2rUKIYNG5Y/D/G6desYPnz4EV0ypT1+4TLa4SW2iyqpDaUvqx0t6lsIisxF3lO+qukvFWzQoEFMmjSJyZMnM3DgQAB27drF8ccfT0JCArNnz2b9+vXFbt+9e3c2bNjAyy+/fNRunMLatGnDmjVr8r/9vvrqq2WK+4wzzmDu3Lns2LGD7OzsYktQF9ajRw8mTZoEUGAugvXr19O4cWNuuOEGrr/+er744otS7a9Pnz489dRTAOTk5BDJWGVKSgpLliwhNzeXDRs28Pnnnx+xTuGS2tGkBBAEPy2Hl1TTXypHhw4d2LNnD82aNcufWGXo0KFkZGSQmprKxIkTadu27VH3ceWVV9KjRw/q169f6uMec8wxPPnkk/Tt25euXbtSp04dyvKwZ7Nmzbjvvvvo1q0bPXr0ICUlpVTbP/bYYzzxxBOkpqaycePG/OVz5syhU6dOdOnShVdffZXbb7+9VHE89thjzJ49m9TUVLp27VqmSe8L69GjBy1btqR9+/bcdtttnH766Ues07BhQ3r06EHHjh2jPgisctB+O7Qf/nm698DXiLe9KwCptqpLOeh+/fpxxx130KdPnzJtt3fvXpKTk3HOccstt9CqVav8u4jKsn12djbp6emMGDGC9PT0soZfrURSDlpXAH7b84P30/sBNf4SeDt37qR169Ycc8wxZW78wbtttHPnznTo0IFdu3Zx4403lmn70aNH59+22bJly1JNTi/FC+QgsJldClx66qmn+h1K5amRWPI6Ij479thjC9x5U1Z33HFHmb7xF6a5f6MrkFcAKgYnIlLxApkARKqzII+7SdUS6e+SEoDf1BjElKSkJLZt26YkIBFzzrFt2zaSkpLKvY9AjgGIVFfNmzcnMzOTLVu2+B2KVANJSUn5T1yXhxJAYKjYWyxISEjIL4Ug4jd1AYmIxCglABGRGKUEICISowJdCsLMtgDFV6U6ukbA1iiGU9GqUrxVKVaoWvFWpVihasUbS7Ge5Jw7rqSVAp0AImFmGaWphREUVSneqhQrVK14q1KsULXiVaxHUheQiEiMUgIQEYlR1TkBFD3XXHBVpXirUqxQteKtSrFC1YpXsRZSbccARETk6KrzFYCIiBxFlUwAZtbXzFaa2Sozu7eIz0eZ2XIz+8rMPjCzk8I+u8bMvgv9XBPwWHPMbEnoZ3pFx1rKeG8ys69DMX1kZu3DPvt9aLuVZnZhUGM1sxQz+zns3P6romMtTbxh6/3KzJyZpYUtC9S5LS7WoJ5bMxtuZlvC4ro+7LOgtQlHizW6bYJzrkr9APHAauBkoCbwJdC+0Dq/AGqFXt8MvBp63QBYE/pv/dDr+kGMNfR+bwDPbd2w15cBb4detw+tnwi0DO0nPqCxpgBLg3ZuQ+vVAeYBnwJpQT23R4k1kOcWGA48XsS2QWwTiow19FlU24SqeAXQDVjlnFvjnDsETAIuD1/BOTfbObc/9PZTIK9c3oXAe8657c65HcB7QN+AxuqH0sS7O+xtbSBvEOlyYJJz7qBzbi2wKrS/IMbqhxLjDfkz8DfgQNiywJ3bo8Tqh9LGW5TAtQmVqSomgGbAhrD3maFlxbkOeKuc20YqklgBkswsw8w+NbPKmPy0VPGa2S1mthp4GLitLNtGUSSxArQ0s8VmNtfMzq3AOPOUGK+ZnQ6c6JybWdZtoyySWCGA5zbkV6Gu1slmdmIZt42WSGKFKLcJVTEBlJqZXQWkAeP8jqUkxcR6kvOeBvw18H9mdoovwRXinHvCOXcKcA/wgN/xHE0xsf4AtHDOdQFGAS+bWV2/YgQwszjgUeB3fsZRGiXEGrhzGzIDSHHOnYb3Lf8Fn+M5mqPFGtU2oSomgI1AeEZsHlpWgJmdD9wPXOacO1iWbaMoklhxzm0M/XcNMAfoUoGxQtnPzyQg71tIIM9tmPxYQ10p20KvF+H1ybauoDjzlBRvHaAjMMfM1gFnAtNDg6tBO7fFxhrQc4tzblvY39YzQNfSbhtlkcQa/TahogY7KuoHbxKbNXiDYXmDKB0KrdMF7xevVaHlDYC1eIM99UOvGwQ01vpAYuh1I+A7ihiI8yHeVmGvLwUyQq87UHCgcg0VO1AZSazH5cWGNxi3sSJ/D0obb6H153B4YDVw5/YosQby3AJNwl6nA5+GXgexTSgu1qi3CRX2P6WC/4dfDHwbajjvDy37E943aID3gZ+AJaGf6WHbjsAbRFsFXBvUWIGzga9DvyBfA9cF5Nw+BiwLxTo7/JcX7ypmNbASuCiosQK/Clv+BXBpEM5toXXnEGpUg3hui4s1qOcWGBOK68vQ70LbsG2D1iYUGWtFtAl6ElhEJEZVxTEAERGJAiUAEZEYpQQgIhKjlABERGKUEoCISIxSApDAMrPGZvayma0xs0Vm9omZpRezblMzm1zMZ3PCK2uWcMz+4RVOI2VmTczszdDrXma2K1TJ8Ssze9/Mji9h+9Fmdmcpj/W+mdWPRtwSG5QAJJDMzIBpwDzn3MnOua7AYIoolmdmNZxzm5xzA6Jw6P541TejZRTw77D3851znZ33mP9C4JZID2CeOOBF4H8i3Z/EDiUACarewCHnXH49eefceufcPyG/Zvp0M/sQ+CBUh35p6LNjzGySma0ws6nAMUUdwMzG2uG5GB4xs7PxykaPC31LPyX083boCmS+mbUNbTvBzP4VKsz1rZn1K+bf8Svg7SKObXglFXaE3jcws2mhWD41s9PCVm8fuopZY2a3hdZPCdWUnwgsxSsvMB0YUtoTLFLD7wBEitEB70nSozkdOM05t93MUsKW3wzsd861CzWkR+zHzBriPWbf1jnnzOxY59zO0CQbbzrnJofW+wC4yTn3nZl1B57ES07g1b7vBpwCzDazU51zB8KO0RLY4cLqOwHnmtkSoCGwD7gvtPyPwGLnXH8z6w1MBDqHPmuLN29EHWClmT0VWt4KuMY592nYMRPNrKEL1eMRORpdAUiVYGZPmNmXZrYwbPF7zrntRax+HvAfAOfcV8BXRayzC6+O/bNmdgWwv/AKZpaM9/j966FG+2mgSdgqrznncp1z3+HVd2lbaBdNgC2FluV1AZ0IPI9XphrgHLwuHJxzHwINw6poznRekbWtwGagcWj5+vDGP2Qz0LSIf6/IEXQFIEG1DK/7BADn3C1m1gjICFtnX3l37pzLNrNuQB9gAHArh7/Z54kDdjrnOhfePm83Jbz/GUg6ShjTgSmlCDf8CiKHw3+3Rf37k0LHFSmRrgAkqD7Em/zi5rBltUq57Ty8eumYWUfgtMIrhL7d13POzQLuADqFPtqD19WC82YUW2tmA0PbmJl1CtvNQDOLC9VkPxmvUFu4b/G6iYpzDl5BMID5wNDQcXoBW13BGc1KFBpXOAFYV5btJHbpCkACKdQv3x/4u5ndjdeVsg9vYpeSPAU8b2YrgBXAoiLWqQO8YWZJgOHdrQPevAH/Dg22DsBrlJ8ysweAhNDnX4bW/R74HKiLN05QYGpE59w+M1sdGhtYFVqcNwZgeN1QeRN+jwaeM7Ov8LqjyjM5eVe80sHZ5dhWYpCqgYqUg5lNIGyw+CjrpQNdnXMVPnOamT2GV078g4o+llQPugIQqUDOuamhO44qw1I1/lIWugIQEYlRGgQWEYlRSgAiIjFKCUBEJEYpAYiIxCglABGRGKUEICISo/4/ozm4mUhyAh4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "plt.xlabel('Grid step (Bohr)')\n",
    "plt.plot(hgrids,energies_h-emin,label='crmult=3.5')\n",
    "plt.plot(hgrids,energies_hc-emin,label='varying hgrids+crmult')\n",
    "plt.yscale('log')\n",
    "plt.legend(loc='best')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fdff10bb6d0>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEKCAYAAAAb7IIBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi40LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcv7US4rQAAIABJREFUeJzt3Xd8FVX+//HXJ8mFUJJQBIQESRSlBqJEEFFAAoqCaJSqVCWUVUHW1XXVdXH9rbLr6ndhwQIiRWkKolRBBQSxYICINAU1LgEURAwdUs7vj7kJ6blpd+bmfp6PRx7kzj0z884Y87lnyjlijEEppZT/CbA7gFJKKXtoAVBKKT+lBUAppfyUFgCllPJTWgCUUspPaQFQSik/pQVAKaX8lBYApZTyU1oAlFLKTwXZHaAol1xyiYmMjLQ7hlJK+ZStW7f+aoypV1w7RxeAyMhIEhMT7Y6hlFI+RUR+8qSdVwuAiNwJ9AJCgZnGmLXe3L9SSqmLPL4GICJviMgREdmZZ3lPEflWRPaLyONFbcMY854xJgEYAwwoXWSllFLloSQ9gNnAVGBu1gIRCQSmAT2AFOArEVkGBALP51n/PmPMEff3T7nXU0opZROPC4AxZqOIROZZ3B7Yb4z5AUBEFgJ3GGOeB3rn3YaICDAJWG2M2VbQfkRkFDAK4LLLLvM0nlJKqRIq622g4cCBHK9T3MsK8xDQHegrImMKamCMmW6MiTXGxNarV+xFbKWUUqXk1YvAxpgpwBRv7lMppVTBytoDOAg0zvE6wr2sTETkdhGZnpqaWtZNKaWUKkRZC8BXwJUiEiUiVYCBwLKyhjLGLDfGjAoLCyvV+lvWvcfmef/gs31HST2TVtY4SilVKXl8CkhEFgBdgUtEJAX4mzFmpog8CKzBuvPnDWPMrgpJWgKyYyGdfl/Noj1bGJ4+ggZ1QmnVMIzW4aG0Cg+jVaNQ6ocE2x1TKaVsJU6cFF5Ebgdub9q0acK+fftKvoHMTM6sfZbqX7zEz6FtmHzJ03z+SxDJx85kN6kfUpXW7mLQqpFVHMJrVcO6UUkppXyXiGw1xsQW286JBSBLbGysKdNQELveg/fGQnAtGPgWJ+u2YfehE+w8dIJdh1LZdfAE+46cJNN9CMKquWgdHkrrRmG0bBRK6/AwourWICBAi4JSyndoAcjy8zew8B44+Qv0mQJtB+Z6+1xaBnt/PsnOg6lWUTh0gr2HT3IhIxOAGlUCadEwNLu30Do8jKb1a+IK1IFUlVLO5NMFoMyngPI6fQzeGQbJm6Djg9D9GQgs/PJHWkYm+345xc5DqVaP4WAquw+f4MyFDACqBAXQ/NIQWjW6WBSaXxpCsCuw7FmVUqqMfLoAZCmXHkCWjDRY8yRseQ2u6AZ934BqtT1fPdOQfOy0u6dgnULaefAEqWetu4wCA4Qr69e0Th01CqN1eBgtGoYQEuwqn/xKKeUhLQCF2TYXVvwRwiJg0AKo36LUmzLGkHL8bPapo50HU9l56ARHT57PbhN1SY1cF5pbNQqjTo0q5fGTKKVUgbQAFOV/X8KiwZB2Bu6aDs17levmj5w4l10Qdh06wc5DqaQcP5v9fqOwYFqFh9HafQqpZaNQ6oVU1esKSqly4dMFoNyvARQk9SAsuhcObYebnoQb/wQBFfcH+PczF9x3IFmnjnYdSuWHX0+T8/DXqBJIaDUXYdVchAa7CK3mIrRaUPbrsGquHO8HXfy+mosaVQL1FlalFODjBSBLhfUAsqSdheXjYcciaNEH7nwFqtasuP3lcfp8OnsOn2DPzyc5fvoCqWfTOHE2zfr3XBqpZ9M54V528nx6kdsKDBBCg4PyFInCi0jettr7UKry8LQAOHpKyArnqgbxr8GlbeDDv8LM72HQfKgd6ZXd16gaRGxkHWIj6xTbNiPTcPJcGifOpucoEFZxuPh97vcO/X42u4hk3dZamGquQHdBKKBY5CkYNaoEERQouAIDcAUKQQEBVAmy/r243P19gNUmMEC0h6KUw/h3DyCn/R/D4hEggdBvNlzexTv79ZJzaRl5ehc5CkZhy92vT54ruvfhKZe7OAQFyMUCkVUsAnIUFPfrKkHWv0GBAVRxty+42Fivs7efY1lQoBAoQkAABIhVhALE+j5AQEQQ3K8DcL9/sY1kt81q70mb3P8W1SZrGUKuPKVRlvpa2r2WaZ9FrFtUnqLXK26fkq9d1vYq0wcU7QGUVNM4SFhvPTT2Zjz0fB7ajyrbb7iDBLsCCXYFUj+05GMgZWQaTp1Lzy4IZy5kkJ6RyYWMTNIzDOmZmVzIMKS7X6dlZpKWnkl6piEtw5CWkUl6RiZpmSbHcndb9/L0jEzre/f20jIMp86nX2yTYa2X63XWvjIMGZnO/SCjfFt2gci1rPBCYi3Pv5IU1a6A/cy5r71HZwfKwpEFIMdFYO/uuO4VMPIjeHc0rH4Mft4BvV6CoKrezeEwgQFCWHUXYdVducb+dpLMTGMViMxM0tKtwpBVLIyBTGPcX9btu5k5lmW9f7Fd4W2ylpHnde72WesbMjPBUESbzJzbAEPpCllZOvKlXbWofRb3c5Q2b1FnLIrbZtbbOdtl5cy9LP9GTf5FuX7GgtYvaD8Usp+snyvnOvVCKv7vjp4CKkhmJmx4Hjb+CyKuhQFvQcil3s+hlFKl4OkpIL31oyABAdDtSeg3B37ZBdO7QspWu1MppVS50gJQlFZ3wv0fQqALZt0KSQvsTqSUUuVGC0BxLm0NCRugcXt4b4w1nlBG+dwVo5RSdnJkAXDcnMA16sKQpdB+NHw+Feb1hTO/2Z1KKaXKxJEFoKxzAleIQBfc9i/oMxV+2gwzusGRPXanUkqpUnNkAXC0a4bA8JXWQHKvd4e9K+1OpJRSpaIFoDQat4dRG+CSq6wHxzb807p1VCmlfIgWgNIKbQQjVkObgbDhOXhnKJw/ZXcqpZTymBaAsnAFQ/yrcMtz1qmgmTfDbz/anUoppTyiBaCsRKDjAzB4CZw4CDNugh8+sTuVUkoVSwtAebmiGySsg5oNrMHkvnytbAO0KKVUBXNkAXDccwCeyhpM7qqe1mByyx6E9PPFr6eUUjZwZAFw5HMAnqoaYg0e1+XPsP0tmN0LTv5sdyqllMrHkQXA5wUEwE1PQP+58MtuHUxOKeVIWgAqUss74P61OpicUsqRtABUtLyDyX3whA4mp5RyBC0A3pA1mFyHMfDFNJh3tw4mp5SynRYAbwl0wa3/hDumwU+fWc8L6GBySikbaQHwtqsHw/BVkHbWGkxuzwq7Eyml/JQWADs0vtYaTK5eM1h0L2yYpIPJKaW8TguAXUIbWT2BtoOsCehfvQH2LNenh5VSXuPIAuCzTwKXlCsY7nwF7p4JGedh0WCY3gW+W6OFQClV4cQ4+A9NbGysSUxMtDuGd2SkwzdvW6eDfv8JwmOh25Nw+U3WgHNKKeUhEdlqjIktrp0jewB+KTAIYu6Bh7bC7ZOt4SPejIdZt0Hyp3anU0pVQloAnCbQBe2Gw7htcNu/4bcfrPGE5vSBA1vsTqeUqkS0ADhVUFVonwDjk6wJZ47shpk94K2+cHCb3emUUpWAFgCnc1WzJpwZ/zV0nwgHE62HyBYMgp+/sTudUsqHaQHwFVVqwA0TYPwOuOlJSN5s3Tr69lA4stfudEopH6QFwNcEh0KXx+DhHdD5Mdi/Dl6+DpaMhF/3251OKeVDtAD4qmq1rNtEH94BncZbk9JPaw/v/QGOJ9udTinlA7QA+LrqdaDHM9Y1gg5jYOcS+G87WD4eUlPsTqeUcjAtAJVFzfrQ8zkYlwTtRsD2eTDlalj1qE5JqZQqkBaAyia0IfT6N4zbbo0zlPgGTG4La56EU0ftTqeUchAtAJVVrcbQZwo8mAit7oIvXrYKwUcTdTIapRTgxQIgIi1E5FURWSwiY721X79XJwriX4EHtkDz2+DT/8B/2sD65+Ds73anU0rZyKMCICJviMgREdmZZ3lPEflWRPaLyONFbcMYs8cYMwboD3QqfWRVKpdcCXe/Dn/4HJp2g0/+CZPbwCcvwPmTdqdTStnA0x7AbKBnzgUiEghMA24FWgKDRKSliESLyIo8X/Xd6/QBVgKryu0nUCVTvwX0nwujN0GTTrD+/1k9gk//AxdO251OKeVFHg8HLSKRwApjTGv3647ARGPMLe7XfwEwxjzvwbZWGmN6FdfOr4aDtsvBrdbpoP0fQY361tPGsfdZcxUopXySN4aDDgcO5Hid4l5WWKCuIjJFRF6jiB6AiIwSkUQRSTx6VO9aqXDh7WDwErhvLdRvDmv+AlNiYMsMSD9vdzqlVAXy2kVgY8wGY8w4Y8xoY8y0ItpNN8bEGmNi69Wr56146rIOMGw5DFsBtSNh1Z+sB8q2zoGMNLvTKaUqQFkKwEGgcY7XEe5lypdF3QgjVsPgd62Hy5aPg6mxkLTAmrVMKVVplKUAfAVcKSJRIlIFGAgsK49QfjMnsFOJQNM4GPkxDFoEVUPhvTHWoHPfLIbMTLsTKqXKgae3gS4APgeaiUiKiNxvjEkHHgTWAHuAt40xu8ojlDFmuTFmVFhYWHlsTpWWCDTrCaM3woC3rNnKltwPr8fB7/+zO51Sqox0UnjlucxM2LkYVv4JAgKh7xtwxU12p1JK5eHTk8LrKSCHCgiANv1h1Hqo2QDeust6fsDBHyKUUoVzZAHQU0AOV/cKGPkRtLwDPvqbNSuZPk2slM9xZAFQPqBqTeg7C3o8C3tXwOvddUYypXyMIwuAngLyESLQaRwMeQ9OH7Umq9+ro3wo5SscWQD0FJCPubwLjPrEOjW0cBCs+wdkZtidSilVDEcWAOWDajWGER9AzGDY+C+YPwDOHrc7lVKqCFoAVPlxBcMdU6H3/8EPG2B6V/h5Z3FrKaVs4sgCoNcAfJiINZroiFXWYHKvd4cd79idSilVAEcWAL0GUAk0bm9dF2h0Nbw7Ej74iw4qp5TDOLIAqEoipAEMWwYdxlhzEs+9E04dsTuVUspNC4CqWIEuuPWfED/dmnzmtS6QosN7KOUEWgCUd7QdAPevhcAgmHUrbJ1tdyKl/J4jC4BeBK6kGraxrgtE3gjLx8Oyh3TWMaVs5MgCoBeBK7HqdeDed+DGR2DbXKs3kJpidyql/JIjC4Cq5AICIe5pa46Bo99Z1wV+3GR3KqX8jhYAZZ8Wt0PCOqtXMPcO+HyaDi2tlBdpAVD2qneVNfVks1thzROwZCRcOG13KqX8giMLgF4E9jPBodbpoLinYecSeL0H/PaD3amUqvQcWQD0IrAfErEuDA9eAicOWuMIfbfW7lRKVWqOLADKjzWNg9GfQK3LYH5/+ORf1lzESqlypwVAOU/tSLhvrTX/8Pp/wKJ74ZyeDlSqvGkBUM5UpTrEvwa3/gv2rYXpN8GRPXanUqpS0QKgnEsEOoyGYcutSednxMGupXanUqrS0AKgnK/J9dZ1gQYt4Z3h8OHTkJFudyqlfJ4jC4DeBqryCW0Ew1dak81sngxv3QWnj9mdSimf5sgCoLeBqgIFVbWmm+wzFf73BUzvAoe2251KKZ/lyAKgVJGuGQL3fWANGzHzFtg+z+5ESvkkLQDKN4VfY10XuKwDvP8HWPkIpF+wO5VSPkULgPJdNS6BwUvh+nHw1eswuxecOGx3KqV8hhYA5dsCg+DmZ6HvLPhll3Vd4KfP7U6llE/QAqAqh9Z3wciPoEoNmNMbvpyuQ0srVQwtAKryaNASEtZD0+6w+lFYOgYunLE7lVKOpQVAVS7VasHABdD1CdixCN64GQ7vsDuVUo6kBUBVPgEB0PXPcM8ia77h1zrDu6Ph9//ZnUwpR3FkAdAngVW5uOoWGLcdOo2zxhD6byysfQrOHrc7mVKOIMbBF8piY2NNYmKi3TFUZfD7AVj/HHy9wJqB7MY/QftR4Aq2O5lS5U5EthpjYotr58gegFLlrlZjiH8FxnwKEdfCh3+FqbHw9UKdcEb5LS0Ayr9c2tqadnLo+1C9DiwdbV0j2P+x3cmU8jotAMo/Xd4VEjbA3TPhfKo1uujcO+Dw1zYHU8p7tAAo/xUQANF94cFEuOV564//a51hSQIc/8nudEpVOC0ASgVVhY5/gHFJcMME2LPMuj6w5kk485vd6ZSqMFoAlMpSrRZ0nwgPbYPo/vD5NJgSA5/+B9LO2p1OqXKnBUCpvMLC4c5pMHYzNO4AH/3NeoYgaT5kZtidTqlyowVAqcI0aAX3vgPDVkDNevDeWHj1Rtj3kQ40pyoFLQBKFSfqRhi5Dvq+AWmnYd7dMLePTkepfJ4WAKU8ERAAre+GB76Cnv90zz3QFZaMhOPJdqdTqlS0AChVEkFV4Lox1hhDNz4Ce1bA1Gvhgyf0jiHlc7QAKFUawWEQ9zSM2wZtBsCXr8DkGNj0kt4xpHyGVwuAiNQQkUQR6e3N/SpVYUIbwR1TYexn0KQjfPwM/LcdbH9L7xhSjudRARCRN0TkiIjszLO8p4h8KyL7ReRxDzb1Z+Dt0gRVytHqt7DmHxi+Emo2gPcfgFdvgO/W6h1DyrE87QHMBnrmXCAigcA04FagJTBIRFqKSLSIrMjzVV9EegC7gSPlmF8pZ4m8ARLWQb/Z1qmg+f1gzu1wcJvdyZTKJ8iTRsaYjSISmWdxe2C/MeYHABFZCNxhjHkeyHeKR0S6AjWwisVZEVlljNFxeFXlIwKt4qFZL9g6Gz6ZBDNusu4i6vZXqBNld0KlAA8LQCHCgQM5XqcAHQprbIx5EkBEhgO/FvbHX0RGAaMALrvssjLEU8pmQVWgwyhoOxA+m2INLbF7GVw7Ejo/CjXq2p1Q+Tmv3wVkjJltjFlRxPvTjTGxxpjYevXqeTOaUhUjOBS6PWWNMRRzD2x5zRpjaNOLcOGM3emUHytLATgINM7xOsK9rMx0TmBVKYU2hD5TYOzn1rWCj/9u3TG07U29Y0jZoiwF4CvgShGJEpEqwEBgWXmEMsYsN8aMCgsLK4/NKeUs9ZvDoAUwYrV1G+myB+GVTvDdGr1jSHmVp7eBLgA+B5qJSIqI3G+MSQceBNYAe4C3jTG7Ki6qUpVMk+th5EfQfy5kXID5/WFeP32iWHmNGAd+4hCR24HbmzZtmrBv3z674yhV8TLS4KvXYe1frQnsBy2Ees3sTqV8lIhsNcbEFtfOkUNB6Ckg5XcCXXDdWBi+As6fhBlx1ikhpSqQIwuAUn7rsusgYb31rMD8AdZsZA7spavKQQuAUk5TqzHctwZa3WnNRrZ0NKSdszuVqoQcWQD0NlDl96pUh76zrOcHdiyC2bfBicN2p1KVjCMLgF4DUAprSInOj8KAeXBkrzUBTcpWu1OpSsSRBUAplUOL3jDyQ2toiVm3wteL7E6kKgktAEr5ggatIGEDRFwLS0fBh0/r08OqzBxZAPQagFIFqFEXhr4HsffD5smwYCCc0/9HVOk5sgDoNQClChHogt4vQa8X4ft18Hp3OPa93amUj3JkAVBKFePakTDkPTj9qzXXwPfr7E6kfJAWAKV8VdSNMGo9hIbDW33hi1f1oTFVIo4sAHoNQCkP1Y6E+9dCs1vhgz/Dsocg/bzdqZSPcGQB0GsASpVA1RDo/yZ0fgy2vwlz+sApnXpbFc+RBUApVUIBAdDtSevp4cNfw/SbrH+VKoIWAKUqk9Z3wX0fAAZm3gK7ltqdSDmYFgClKptGMTBqAzRsA+8Mh3X/gMxMm0MpJ9ICoFRlVLM+DFsOMYNh47/g7SFw/pTdqZTDOLIA6F1ASpWDoKpwx1ToOQm+XQUzb4bjyXanUg7iyAKgdwEpVU5ErJnG7l0MJ1JgRjdI/tTuVMohHFkAlFLlrGmcNdNY9bow9w5IfMPuRMoBtAAo5S/qXgEjP4IrusGKCbDyEWsyeuW3tAAo5U+Cw2DQQug0Hr56Hd6Mh9PH7E6lbKIFQCl/ExAIPf4O8dPhwBZrMLlfdtudStlAC4BS/qrtABix2ho7aGYP2LvS7kTKyxxZAPQ2UKW8JKKdNaLoJVfBwntg4ws6oqgfcWQB0NtAlfKi0EYwYhVE94d1/w8W3wcXztidSnmBIwuAUsrLXNXgrunQ/Rlr/KBZPSE1xe5UqoJpAVBKWUTghofhnkVw7AdrRNH/fWl3KlWBtAAopXK76hZI+Biq1oQ5vWH7W3YnUhVEC4BSKr96zWDkx9Dkenj/AfjgCchItzuVKmdaAJRSBateB+5dAh3GwBfTYH4/OHvc7lSqHGkBUEoVLjAIbv0n3D4FftwEM+Lg6Hd2p1LlRAuAUqp47YZZ8wucS4XX42Dfh3YnUuVAC4BSyjNNOlozjdVuAvP6weYp+tCYj3NkAdAngZVyqFqN4b410LIPfPhXWPes3YlUGTiyAOiTwEo5WJUa0G8OXDMMNr0ImyfbnUiVUpDdAZRSPkgEev8fnD8JHz5tDTPdbrjdqVQJaQFQSpVOQCDEvwYXTsHyh6FqCLS+2+5UqgQceQpIKeUjgqpYp4OaXA/vjoLv1tqdSJWAFgClVNlUqW7NMtagNbw9BJI3251IeUgLgFKq7IJDYfC7UKsJzB8Ah7bbnUh5QAuAUqp81KgLQ9+D6rXhzbvgyF67E6liaAFQSpWf0EYw5D0IdMGbd8LxZLsTqSJoAVBKla+6V8CQpZB2FubeCSd/tjuRKoQWAKVU+WvQCgYvgVNH4M14OPOb3YlUAbQAKKUqRkQsDFoAx76HeX2th8aUo2gBUEpVnMu7QL9ZcCgJFgyCtHN2J1I5aAFQSlWs5r3gzlcgeRMsHgEZaXYnUm5eKwAi0lVENonIqyLS1Vv7VUo5QNsBcNu/4dtV8N4fIDPT7kQKDwuAiLwhIkdEZGee5T1F5FsR2S8ijxezGQOcAoKBlNLFVUr5rPYJ0O2v8M3bsPpRnUvAATwdDG42MBWYm7VARAKBaUAPrD/oX4nIMiAQeD7P+vcBm4wxn4hIA+Al4N6yRVdK+ZwbH7FmFftsijWCaNzTdifyax4VAGPMRhGJzLO4PbDfGPMDgIgsBO4wxjwP9C5ic8eBqiWPqpTyeSLQ4+9WEdj0IlQNhRsetjuV3yrLcNDhwIEcr1OADoU1FpG7gFuAWli9icLajQJGAVx22WVliKeUcqSccwl89DerJxA7wu5Ufslr8wEYY94F3vWg3XRgOkBsbGy+k4RpaWmkpKRw7pzeTqbKLjg4mIiICFwul91R/EvOuQRWTLDmEojua3cqv1OWAnAQaJzjdYR7WZmJyO3A7U2bNs33XkpKCiEhIURGRiIi5bE75aeMMRw7doyUlBSioqLsjuN/suYSmNcXlo62isBVt9idyq+U5TbQr4ArRSRKRKoAA4Fl5RGqqDmBz507R926dfWPvyozEaFu3bram7RTrrkEhkLyp3Yn8iue3ga6APgcaCYiKSJyvzEmHXgQWAPsAd42xuyquKi58nhjN8oP6O+SA+SaS2AgHNxmdyK/4VEBMMYMMsY0NMa4jDERxpiZ7uWrjDFXGWOuMMb8o7xCicjtIjI9NTW1vDZZrpKTk2ndunWp1z906BB9+xZ8vrNr164kJiZ6tJ3z588zYMAAmjZtSocOHUhOTi6wXWRkJNHR0cTExBAbG5u9fOLEiYSHhxMTE0NMTAyrVq0q8c+iVLnIOZfAW3frXAJe4sihIIo6BeTr0tPTadSoEYsXLy7ztmbOnEnt2rXZv38/EyZM4M9//nOhbdevX09SUlK+4jJhwgSSkpJISkritttuK3MmpUpN5xLwOkcWAF+QkZFBQkICrVq14uabb+bs2bN89dVXtGnThpiYGB599NHsXsLs2bPp06cP3bp1Iy4uLlcP4uzZswwcOJAWLVoQHx/P2bNnPc7w/vvvM2zYMAD69u3Lxx9/jNGnK5UvyzWXwB06l0AF89ptoCVR1F1AOT2zfBe7D50o1323bBTK325vVWy7ffv2sWDBAmbMmEH//v1ZsmQJkyZNYsaMGXTs2JHHH889Msa2bdvYsWMHderUyXWq5pVXXqF69ers2bOHHTt2cM0112S/N2DAAL799tt8+/7jH//I0KFDOXjwII0bWzdiBQUFERYWxrFjx7jkkktytRcRbr75ZkSE0aNHM2rUqOz3pk6dyty5c4mNjeXFF1+kdu3aHh0npSpMg1bWNYG5fawJZUasgup17E5VKTmyB+ALp4CioqKIiYkBoF27diQnJ3Py5Ek6duwIwD333JOrfY8ePahTJ/8v8caNGxk8eDAAbdq0oU2bNtnvLVq0KPv0TM6voUOHlijrp59+yrZt21i9ejXTpk1j48aNAIwdO5bvv/+epKQkGjZsyCOPPFKi7SpVYSLaWXMJ/PaDziVQgRzZA/CUJ5/UK0rVqhdHswgMDOTw4cNFtq9Ro0aJ91FcDyA8PJwDBw4QERFBeno6qamp1K1bN1/78PBwAOrXr098fDxbtmyhc+fONGjQILtNQkICvXsXNYKHUl4W1Rn6zYZFg625BO5dDK5gu1NVKo7sAfiiWrVqERISwpdffgnAwoULPVqvc+fOzJ8/H4CdO3eyY8eO7PeK6wH06dOHOXPmALB48WK6deuW77bG06dPc/Lkyezv165dm339IWfRWrp0aZnubFKqQjS/DeJftZ4P0LkEyp0jewCeXgNwmpkzZ5KQkEBAQABdunTBk1NYY8eOZcSIEbRo0YIWLVrQrl07j/d3//33M2TIEJo2bUqdOnWyi86hQ4cYOXIkq1at4pdffiE+Ph6w7kC655576NmzJwCPPfYYSUlJiAiRkZG89tprpfiplapgbfpbg8et+pM1l0D8axCgn13Lgzj5rpHY2FiT97bFPXv20KJFC5sSFe3UqVPUrFkTgEmTJnH48GEmT55scypVHCf/TqkcNr0IH/8drh1pTS6jD/EVSkS2GmNii2vnyB6Ar1q5ciXPP/886enpNGnShNmzZ9sdSanK44Y/Wj2BzZN1LoFyogWgHA0YMIABAwbYHUOpykkEuj+jcwmUIy0ASinfIQK9XtK5BMqJIwuAr14EVkp5QdZcAud1LoGycuSldF94EEwpZaNAF/SfA02ut+YS+PYDuxP5JEcWAKWUKpar2sW5BN4ZBj9usjvCgR6bAAAOgklEQVSRz9EC4CDXX399uWxn9uzZPPjgg2XaxrJly5g0aVKB72Xd6up0kZGR/Prrr/z++++8/PLLdsdRFSHnXAILBsLBrXYn8ilaALzIGENmZmah73/22WdeTFO49PR0+vTpk29Au5JKTk6ma9euJV6vuONUUloAKrnsuQTquOcS2GN3Ip/hyALg5AlhHn/8caZNm5b9euLEifz73//m1KlTxMXFcc011xAdHc37778PWH8EmzVrxtChQ2ndujXPPvssDz988da1GTNmMGHCBODiJ+sNGzbQtWtX+vbtS/Pmzbn33nuzh3letWoVzZs3p127dowbN67Q8XsOHTpEz549ufLKK3nssceyl8+cOZOrrrqK9u3bk5CQkN1TGD58OGPGjKFDhw489thjuXoRP/74Ix07diQ6Opqnnnoqe1uHDx+mc+fOxMTE0Lp1azZt8qwLnvV0ctu2bWnbti2fffZZvuN04MABatasyaOPPkqrVq3o3r07W7ZsoWvXrlx++eUsW2bNPpq3t9O7d282bNiQ77/Z999/nz1Mt6qEQhvB0PchsIo1guhvP9qdyCc48i4gY8xyYHlsbGxCkQ1XPw4/f1O+O780Gm4t+NQHWPf6P/zwwzzwwAMAvP3226xZs4bg4GCWLl1KaGgov/76K9dddx19+vQBrKGj58yZw3XXXcepU6do27YtL7zwAi6Xi1mzZhU4BMP27dvZtWsXjRo1olOnTmzevJnY2FhGjx7Nxo0biYqKYtCgQYXmTEpKYvv27VStWpVmzZrx0EMPERgYyLPPPsu2bdsICQmhW7dutG3bNnudlJQUPvvsMwIDA3M9xDZ+/HjGjh3L0KFDcxW/+fPnc8stt/Dkk0+SkZHBmTNnPDrE48aNo0uXLixdupSMjAxOnTrF8ePHcx0nsMYu6tatGy+88ALx8fE89dRTfPjhh+zevZthw4ZlH9/iTJo0iZ07d5KUlORRe+Wj6lxuTSgz61ZrQpkRH0BoQ7tTOZojewBOdvXVV3PkyBEOHTrE119/Te3atWncuDHGGJ544gnatGlD9+7dOXjwIL/88gsATZo0yf6jVrNmTbp168aKFSvYu3cvaWlpREdH59tP+/btiYiIICAggJiYGJKTk9m7dy+XX345UVFRAEUWgLi4OMLCwggODqZly5b89NNPbNmyhS5dulCnTh1cLhf9+vXLtU6/fv0IDAzMt63Nmzdn72vIkCHZy6+99lpmzZrFxIkT+eabbwgJCQEgPj6emJgYbrvtNhITE7OnnJw1axYA69atY+zYsYA1kmrW3V45jxNAlSpVssctio6OpkuXLrhcLqKjowud/lL5uQYtrWsCp3+FN+PhzG92J3I0R/YAPFbEJ/WK1K9fPxYvXszPP/+c/eTvvHnzOHr0KFu3bsXlchEZGcm5c+eA/ENBjxw5kueee47mzZszYkTBD7HkHW46PT29RBlLs35RQ1YXNHl6586d2bhxIytXrmT48OHZw1QvXboUsE5/DR8+PN8pGU/373K5svcbEBCQ/TMFBARk/zxBQUG5rhdkHXPlx7LmEnirr3VNYNgy61kBlY/2AEphwIABLFy4kMWLF2d/ik5NTaV+/fq4XC7Wr1/PTz/9VOj6HTp04MCBA8yfP7/IT/F5NWvWjB9++CH70++iRYtKlPvaa6/lk08+4fjx46Snp7NkyRKP1uvUqVP2SKPz5s3LXv7TTz/RoEEDEhISGDlyJNu2bfNoe3FxcbzyyiuANbVmWa71REZGkpSURGZmJgcOHGDLli352oSEhGQPia38RNZcAoe/tuYSSPN8qlV/ogWgFFq1asXJkycJDw+nYUPrHOO9995LYmIi0dHRzJ07l+bNmxe5jf79+9OpU6cSTcFYrVo1Xn75ZXr27Em7du0ICQnxaMjpLOHh4TzxxBO0b9+eTp06ERkZ6dH6kydPZtq0aURHR3Pw4MHs5Rs2bKBt27ZcffXVLFq0iPHjx3uUY/Lkyaxfv57o6GjatWvH7t27Pf4Z8urUqRNRUVG0bNmScePG5ZpSM0vdunXp1KkTrVu31ovA/iTnXALv6FwCBTLGOParXbt2Jq/du3fnW+aLevXqZT766KMSr3fy5EljjDGZmZlm7Nix5qWXXirV+mlpaaZ3797m3XffLXGGyqay/E6pQnw53Zi/hRqz+i92J/EaINF48DdWewBe9vvvv3PVVVdRrVo14uLiSrz+jBkziImJoVWrVqSmpjJ69OgSrT9x4sTs2zajoqK48847S5xBKZ/SPgEu6wiHPDtF6U8cOSFMjsHgEvbt25frPZ28Q5U3/Z3yA3Nut04B3ecfYwZ5OiGMI3sARgeDU0qpCufIAlAcJ/ZalG/S3yXlz3yuAAQHB3Ps2DH9H1eVmTGGY8eOERwcbHcU5Q36NyMfn3sQLCIigpSUFI4ePWp3FFUJBAcHExERYXcMpWzhcwXA5XJlD4WglFKeyf8ku/LBU0BKKaXKhxYApZTyU1oAlFJ+Qi8C5+XIB8GyiMhRoPBR1Yp2CfBrOcapaL6U15eygm/l9aWs4Ft5/SlrE2NMveIaOboAlIWIJHryJJxT+FJeX8oKvpXXl7KCb+XVrPnpKSCllPJTWgCUUspPVeYCMN3uACXkS3l9KSv4Vl5fygq+lVez5lFprwEopZQqWmXuASillCqCTxcAEQkWkS0i8rWI7BKRZwpoM1xEjopIkvtrpB1Zc+QJFJHtIrKigPeqisgiEdkvIl+KSKT3E+bLVFRexxxbEUkWkW/cORILeF9EZIr72O4QkfxzR3qRB3m7ikhqjmP7tB053VlqichiEdkrIntEpGOe9512bIvL64hjKyLNcmRIEpETIvJwnjYVemx9biygPM4D3Ywxp0TEBXwqIquNMV/kabfIGPOgDfkKMh7YA4QW8N79wHFjTFMRGQj8ExjgzXAFKCovOOvY3mSMKeze6VuBK91fHYBX3P/aqai8AJuMMb29lqZwk4EPjDF9RaQKUD3P+047tsXlBQccW2PMt0AMWB+0gIPA0jzNKvTY+nQPwD395Sn3S5f7y7EXNUQkAugFvF5IkzuAOe7vFwNxImLbKFYe5PUldwBz3b8zXwC1RKSh3aGcTkTCgM7ATABjzAVjzO95mjnm2HqY14nigO+NMXkffK3QY+vTBQCyT1EkAUeAD40xXxbQ7G5392mxiDT2csSc/gM8BmQW8n44cADAGJMOpAJ1vROtQMXlBeccWwOsFZGtIjKqgPezj61binuZXYrLC9DRfXpztYi08ma4HKKAo8As96nA10WkRp42Tjq2nuQFZxzbnAYCCwpYXqHH1ucLgDEmwxgTA0QA7UWkdZ4my4FIY0wb4EMufsL2KhHpDRwxxmy1Y/8l5WFeRxxbtxuMMddgdZkfEJHONmbxRHF5t2E9zt8W+C/wnrcDugUB1wCvGGOuBk4Dj9uUxROe5HXKsQXAfZqqD/COt/ft8wUgi7ubtx7omWf5MWPMeffL14F23s7m1gnoIyLJwEKgm4i8lafNQaAxgIgEAWHAMW+GzKHYvA46thhjDrr/PYJ1HrV9nibZx9Ytwr3MFsXlNcacyDq9aYxZBbhE5BKvB7U+cabk6FkvxvoDm5OTjm2xeR10bLPcCmwzxvxSwHsVemx9ugCISD0RqeX+vhrQA9ibp03O82V9sC5oep0x5i/GmAhjTCRWd2+dMWZwnmbLgGHu7/u629hyTcOTvE45tiJSQ0RCsr4HbgZ25mm2DBjqvqviOiDVGHPYy1EBz/KKyKVZ139EpD3W/6te/zBgjPkZOCAizdyL4oDdeZo55th6ktcpxzaHQRR8+gcq+Nj6+l1ADYE57ivoAcDbxpgVIvJ3INEYswwYJyJ9gHTgN2C4bWkLkCfrTOBNEdmPlXWgreEK4NBj2wBY6v5/OgiYb4z5QETGABhjXgVWAbcB+4EzwAibsoJnefsCY0UkHTgLDLTrwwDwEDDPfariB2CEg48tFJ/XMcfW/QGgBzA6xzKvHVt9ElgppfyUT58CUkopVXpaAJRSyk9pAVBKKT+lBUAppfyUFgCllPJTWgCUz3Lfz71QRL53D6mwSkSuqqB9nXL/Gyki95TTNiNFJO/zCkp5jRYA5ZPcD/IsBTYYY64wxrQD/oJ1j33OduX9rEskUC4FQCm7aQFQvuomIM39sAwAxpivjTGbxBrvfZOILAN2uz9p7xWR2SLynYjME5HuIrJZRPa5nwZFRCaKyJ+yticiOyX/nAyTgBvFGr99Qs433L2RXjlezxaRvu79bxKRbe6v6/P+MGLNrTA1x+sVItLV/f3NIvK5e913RKRmGY6bUtm0AChf1RooaqC6a4DxxpisU0JNgReB5u6ve4AbgD8BT5Rgv49jjSUfY4z5vzzvLQL6Q/YAX3HASqyRanu4B38bAEzxdGfuMWqeArq7108E/liCvEoVyteHglCqMFuMMT/meP2jMeYbABHZBXxsjDEi8g3WaZ3ysBqYLCJVsQYl3GiMOSvWGPVTRSQGyABKcp3iOqAlsNk9dEQV4PNyyqv8nBYA5at2YY3pUpjTeV6fz/F9Zo7XmVz8/yCd3L3i4JIEMsacE5ENwC1Yn/QXut+aAPwCtHVv/1wBqxe2b8Ga52JQSbIo5Qk9BaR81TqgquSYTEVE2ojIjWXYZjLuoYPFmns1qoA2J4GQIraxCGvArhuBD9zLwoDDxphMYAgQWMi+Y0QkQKyJdbKGh/4C6CQiTd25alTUnU7K/2gBUD7JPXpjPNDdfRvoLuB54OcybHYJUMe9rQeB7wposwPIEGs2qQkFvL8W6AJ8ZIy54F72MjBMRL7Guv6Qt3cCsBn4EWvo4ilYk5ZgjDmKNcrqAhHZgXX6p3npfjylctPRQJVSyk9pD0AppfyUFgCllPJTWgCUUspPaQFQSik/pQVAKaX8lBYApZTyU1oAlFLKT2kBUEopP/X/ARurh7sNU8FCAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.xlabel('Crmult value')\n",
    "plt.plot(crmult,energies_c-emin,label='hgrid=0.55')\n",
    "plt.plot(crmult,energies_hc-emin,label='varying hgrids+crmult')\n",
    "plt.yscale('log')\n",
    "plt.legend(loc='best')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can the analyze results exactly as in a non-AiiDa notebook, similarly to our previous [tutorial](./CH4.ipynb) on the methane molecule."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}