{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dissociation curve\n",
    "\n",
    "MyBigDFT comes with some classes implementing particular workflows of calculations. These workflows define a queue of jobs, that can easily be run sequentially, without having to worry about the Job context manager. They also generally define a particular post-processing procedure, run after all the BigDFT calculations in order to extract some meaningful imformation.\n",
    "\n",
    "The example provided here shows how to obtain the dissociation curve of the N$_2$ molecule by using the [Dissocation](https://mmoriniere.gitlab.io/MyBigDFT/dissociation.html) class."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initialization\n",
    "\n",
    "You first need to import some useful classes to define a ground state calculation as well as the `Dissociation` class, that is a workflow:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from mybigdft import Posinp, Atom, InputParams\n",
    "from mybigdft.workflows import Dissociation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The [Dissocation](https://mmoriniere.gitlab.io/MyBigDFT/dissociation.html) class\n",
    "\n",
    "This class allows to compute the dissocation curve of two subsystems. To that end, you need to define two subsystems via ``Posinp`` instances (here, both subsystems are the same nitrogen atom) and a set of distances between those two subsystems. From these data, the ``Dissociation`` class defines a queue of jobs, that differ only by the initial positions: for each specified distance between both subsytem, the second fragment is translated by that amount along the $y$ direction. This means you could study the dissociation curve of two surfaces or of a molecule or of a single atom on top of that surface. \n",
    "\n",
    "Of course, you can specify the input parameters you want to use for all these calculations, via the ``inputparams`` argument. You may also want to give a specific name to the runs and a specific folder where to run them, as usual."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "frag1 = Posinp([Atom('N', [0.0, 0.0, 0.0])], units=\"angstroem\", boundary_conditions=\"free\")\n",
    "distances = np.arange(0.95, 1.25, 0.05)\n",
    "inp = InputParams()  # You might want to specify some non-default input parameters\n",
    "dc = Dissociation(frag1, frag1, distances, inputparams=inp,\n",
    "                  name=\"N2\", run_dir=\"N2/dissociation_curve\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Running the workflow is done as usual:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/dissociation_curve/y_0.95\n",
      "Logfile log-N2.yaml already exists!\n",
      "\n",
      "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/dissociation_curve/y_1.0\n",
      "Logfile log-N2.yaml already exists!\n",
      "\n",
      "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/dissociation_curve/y_1.05\n",
      "Logfile log-N2.yaml already exists!\n",
      "\n",
      "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/dissociation_curve/y_1.1\n",
      "Logfile log-N2.yaml already exists!\n",
      "\n",
      "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/dissociation_curve/y_1.1500000000000001\n",
      "Logfile log-N2.yaml already exists!\n",
      "\n",
      "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/dissociation_curve/y_1.2000000000000002\n",
      "Logfile log-N2.yaml already exists!\n",
      "\n",
      "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/dissociation_curve/y_1.2500000000000002\n",
      "Logfile log-N2.yaml already exists!\n",
      "\n"
     ]
    }
   ],
   "source": [
    "dc.run(nmpi=6, nomp=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can check that the initial positions are as expected:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/dissociation_curve/y_0.95 N2\n",
      "2   angstroem\n",
      "free\n",
      "N   0.0   0.0   0.0\n",
      "N   0.0   0.95   0.0\n",
      "\n",
      "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/dissociation_curve/y_1.0 N2\n",
      "2   angstroem\n",
      "free\n",
      "N   0.0   0.0   0.0\n",
      "N   0.0   1.0   0.0\n",
      "\n",
      "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/dissociation_curve/y_1.05 N2\n",
      "2   angstroem\n",
      "free\n",
      "N   0.0   0.0   0.0\n",
      "N   0.0   1.05   0.0\n",
      "\n",
      "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/dissociation_curve/y_1.1 N2\n",
      "2   angstroem\n",
      "free\n",
      "N   0.0   0.0   0.0\n",
      "N   0.0   1.1   0.0\n",
      "\n",
      "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/dissociation_curve/y_1.1500000000000001 N2\n",
      "2   angstroem\n",
      "free\n",
      "N   0.0   0.0   0.0\n",
      "N   0.0   1.15   0.0\n",
      "\n",
      "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/dissociation_curve/y_1.2000000000000002 N2\n",
      "2   angstroem\n",
      "free\n",
      "N   0.0   0.0   0.0\n",
      "N   0.0   1.2   0.0\n",
      "\n",
      "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/dissociation_curve/y_1.2500000000000002 N2\n",
      "2   angstroem\n",
      "free\n",
      "N   0.0   0.0   0.0\n",
      "N   0.0   1.25   0.0\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for job in dc.queue:\n",
    "    print(job.run_dir, job.name)\n",
    "    print(job.posinp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ``Dissociation`` class also has an ``energies`` attribute, storing the energy of each calculation. This is useful to plot the dissociation curve (see below)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[-19.805444659275025,\n",
       " -19.85497382791818,\n",
       " -19.878933352041976,\n",
       " -19.884549270710195,\n",
       " -19.877164837418295,\n",
       " -19.86087438302971,\n",
       " -19.838574516454937]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dc.energies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A more readable summary can be done as follows (thanks to the ``distance`` attribute of each job of the queue):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "d = 0.95 angstroem: -19.8054 Ha\n",
      "d = 1.00 angstroem: -19.8550 Ha\n",
      "d = 1.05 angstroem: -19.8789 Ha\n",
      "d = 1.10 angstroem: -19.8845 Ha\n",
      "d = 1.15 angstroem: -19.8772 Ha\n",
      "d = 1.20 angstroem: -19.8609 Ha\n",
      "d = 1.25 angstroem: -19.8386 Ha\n"
     ]
    }
   ],
   "source": [
    "for job in dc.queue:\n",
    "    print(f\"d = {job.distance:.2f} {job.posinp.units}: {job.logfile.energy:.4f} Ha\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The job with the minimal energy is also stored under the ``minimum`` attribute:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.1, -19.884549270710195)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dc.minimum.distance, dc.minimum.logfile.energy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is then easy to plot the actual dissociation curve:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZoAAAEKCAYAAAArYJMgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XtY1GXex/H3DIgoCHhC5egJlfMEA460HldDHytLszKtFMms1exJJbfd1l13XWtzF7eThVoeUitT0tLVMMXSVBBFTTRPoGgqeAAUVBjmfv5w41kWEFSG3wDf13VxXfI7zHxuBufL7/7dc986pZRCCCGEsBK91gGEEEI0bFJohBBCWJUUGiGEEFYlhUYIIYRVSaERQghhVVJohBBCWJUUGiGEEFYlhUYIIYRVSaERQghhVfZaB7AFbdq0oWPHjlrHEEKIeiUrK4uLFy9We5wUGqBjx47s2bNH6xhCCFGvGI3GGh0nXWdCCCGsSgqNEEIIq5JCI4QQwqqk0AghhLAqKTRCCCGsSgqNEEIIq5JCI4QQwqqk0NyDE7nXmPVVBsVmi9ZRhBDCZkmhuQenLxXx0Y5MkjIuaB1FCCFslhSae9CnW1u8Wjbjk12ntI4ihBA2SwrNPbDT63iqpw87T17ieM41reMIIYRNkkJzjx43etPETsfy3XJVI4QQlZFCc4/aODdlcFAHVqed4XpxqdZxhBDC5kihqQVjevpQcMPMV/t/1jqKEELYHCk0tSCyUyu6tXPmE+k+E0KICqTQ1AKdTsfonr4cOJPPgTN5WscRQgibokmhWbVqFYGBgej1+nILjhUXFzNu3DiCg4MJDQ0lOTm50vPT09MxmUwYDAaMRiMpKSkAHDlyhF69etG0aVPmzp1bF00p82iYJ82a2MlQZyGE+C+aFJqgoCDWrFlDnz59ym1fsGABAAcPHiQpKYmpU6disVT81H1cXBwzZ84kPT2dWbNmERcXB0CrVq14++23mTZtmvUb8V9cHJvwyH0erNv/M/lFJXX+/EIIYas0KTT+/v507969wvaMjAwGDBgAgLu7O25ubpUusazT6SgoKAAgPz8fDw+PsnMiIiJo0qSJFdNXbXRPX26UWFi994wmzy+EELbIpu7RhIaGsm7dOsxmM5mZmaSlpZGdnV3huHnz5jF9+nS8vb2ZNm0ac+bM0SBtRUGerhi83Vi++xRKKa3jCCGETbBaoRk4cCBBQUEVvtauXVvlOTExMXh5eWE0Gnn55ZeJiorCzs6uwnHz588nPj6e7Oxs4uPjGT9+/B3nS0hIwGg0YjQayc3NvePzqzLG5MuJ3EJ2nrxUa48phBD1mb21Hnjz5s13fI69vT3x8fFl30dFRdGtW7cKxy1ZsoR//vOfAIwcOZLY2Ng7fq4JEyYwYcIEAIxG4x2fX5UHQzrw568zWL77NFFd2tTa4wohRH1lU11nRUVFFBYWApCUlIS9vT0BAQEVjvPw8GDbtm0AbNmyBT8/vzrNeTuOTewYGe7Fph/Pk3P1htZxhBBCc1a7ormdxMREJk+eTG5uLkOHDsVgMLBp0yZycnKIjo5Gr9fj6enJsmXLys6JjY1l4sSJGI1GFixYwJQpUzCbzTg6OpKQkADA+fPnMRqNFBQUoNfrmTdvHhkZGbi4uNRp+0abfFm4PZPPU7OZNMB2iqAQQmhBp+SuNUajsdLRbfdizMLdZF4s5Lu4/tjpdbX62EIIYQtq+t5pU11nDcnonj6czbvO1iM5WkcRQghNSaGxkoEB7XBv0VTmPxNCNHpSaKykiZ2eJyN92HY0l+zLRVrHEUIIzUihsaJRkd7odTqW7z6tdRQhhNCMFBor6uDajF/3cOfzPdncNMuiaEKIxkkKjZWNMflyubCYjT+e1zqKEEJoQgqNlf2qaxt8WzeX5QOEEI2WFBor0+t1jO7pQ2rWFY6cL9A6jhBC1DkpNHVgZLg3DvZ6VsigACFEIySFpg60dHLgweAOrNl7lsKbZq3jCCFEnZJCU0dGm3y5dtPM2vSftY4ihBB1SgpNHQnzccO/gwuf7JJF0YQQjYsUmjqi0+kYY/Ih41wB+7LztI4jhBB1RgpNHRpm8MTJwU6GOgshGhUpNHXIuak9j4Z58vWBc1wpLNY6jhBC1AkpNHVsjMmXYrOFL9LOaB1FCCHqhCaFZtWqVQQGBqLX68stmlNcXMy4ceMIDg4mNDSU5OTkSs9PT0/HZDJhMBgwGo2kpKQAsHz5ckJCQggODiYqKor9+/fXRXPuSI/2Lhh9W7J89yksFhkUIIRo+DQpNEFBQaxZs4Y+ffqU275gwQIADh48SFJSElOnTsVisVQ4Py4ujpkzZ5Kens6sWbOIi4sDoFOnTmzbto2DBw/y+uuvM2HCBOs35i6MMfmSdamIHScuah1FCCGsTpNC4+/vT/fu3Stsz8jIYMCAAQC4u7vj5uZW6TKhOp2OgoJb07nk5+fj4eEBQFRUFC1btgTAZDJx5oxtdk8NCW5PKycHGRQghGgUbOoeTWhoKOvWrcNsNpOZmUlaWhrZ2dkVjps3bx7Tp0/H29ubadOmMWfOnArHLFq0iCFDhtRF7DvW1N6OkUYvNh/O4Xz+Da3jCCGEVVmt0AwcOJCgoKAKX2vXrq3ynJiYGLy8vDAajbz88stERUVhZ2dX4bj58+cTHx9PdnY28fHxjB8/vtz+rVu3smjRIt58880qnyshIQGj0YjRaCQ3N/fuG3qXRkf6YlGKlSky/5kQomHTKQ0/pt6vXz/mzp2L0WisdH9UVBQLFy4kICCg3HZXV1fy8vLQ6XQopXB1dS3rSjtw4ACPPvoo//rXv+jWrVuNchiNxkq76Kzt2Y9SOHK+gB2vDsDezqYuLoUQolo1fe+0qXe3oqIiCgsLAUhKSsLe3r5CkQHw8PBg27ZtAGzZsgU/Pz8ATp8+zfDhw1m2bFmNi4yWxph8uVBwk82Hc7SOIoQQVmOvxZMmJiYyefJkcnNzGTp0KAaDgU2bNpGTk0N0dDR6vR5PT0+WLVtWdk5sbCwTJ07EaDSyYMECpkyZgtlsxtHRkYSEBABmzZrFpUuXePHFF281zt5ekyuVmhrQwx0PV0eW7z7F4KD2WscRQgir0LTrzFZo1XUG8M63x/h70lG2TutHpzZOmmQQQoi7US+7zhqjJyK8sdfrWLFbhjoLIRomKTQac3dx5IHAdqxKO8ONklKt4wghRK2TQmMDxvT0Ja+ohPUHzmkdRQghap0UGhvQq0trOrd14hPpPhNCNEBSaGyATqdjdE9f9p3O49DP+VrHEUKIWiWFxkY8FuaFYxM9n+ySmQKEEA2LFBob4dq8CQ+FeLA2/SxXb5RoHUcIIWqNFBobMsbkS1FxKYn7zmodRQghao0UGhsS6u1GsKcry3edRj5HK4RoKKTQ2JgxJh9+unCVPaeuaB1FCCFqhRQaG/NQqActHO1lUTQhRIMhhcbGNHewZ0SYF/86eJ5L125qHUcIIe6ZFBobNMbkQ3Gphc/32OZS1EIIcSek0Nigru4tMHVuxYqUU1gsMihACFG/SaGxUaN7+pJ9+TrbjtX9MtNCCFGbpNDYqOjA9rRxbspyGRQghKjnNCk0q1atIjAwEL1eX27RnOLiYsaNG0dwcDChoaEkJydXen56ejomkwmDwYDRaCQlJQWAtWvXEhISUrZ9+/btddEcq3Cw1/NEhBdbjuRwNu+61nGEEOKuaVJogoKCWLNmDX369Cm3fcGCBQAcPHiQpKQkpk6disViqXB+XFwcM2fOJD09nVmzZhEXFwfAr3/9a/bv3096ejofffQRsbGx1m+MFY2K9EEBK3fL/GdCiPpLk0Lj7+9P9+7dK2zPyMhgwIABALi7u+Pm5lbpMqE6nY6CggIA8vPz8fDwAMDZ2RmdTgdAYWFh2b/rK6+WzRnQ3Z1PU7MpNlcsuEIIUR/Y1D2a0NBQ1q1bh9lsJjMzk7S0NLKzsyscN2/ePKZPn463tzfTpk1jzpw5ZfsSExPp0aMHQ4cO5aOPPqrL+FYxxuTLxWs3+SbjvNZRhBDirlit0AwcOJCgoKAKX2vXrq3ynJiYGLy8vDAajbz88stERUVhZ2dX4bj58+cTHx9PdnY28fHxjB8/vmzfo48+ypEjR/jyyy95/fXXq3yuhIQEjEYjRqOR3FzbHdnVp1tbvFo2k5kChBD1lk5pOHtjv379mDt3LkajsdL9UVFRLFy4kICAgHLbXV1dycvLQ6fToZTC1dW1rCvtP3Xu3JmUlBTatGlz2xxGo7HSLjpb8X7ycf628Sc2v9KXru7OWscRQgig5u+dNtV1VlRURGFhIQBJSUnY29tXKDIAHh4ebNu2DYAtW7bg5+cHwPHjx8tmPd67dy83b96kdevWdZTeeh43etPETsdyWepZCFEP2WvxpImJiUyePJnc3FyGDh2KwWBg06ZN5OTkEB0djV6vx9PTk2XLlpWdExsby8SJEzEajSxYsIApU6ZgNptxdHQkISEBgNWrV7N06VKaNGlCs2bN+Oyzz+r9gACANs5NGRLUgdVpZ4iL7kEzh4rdiUIIYas07TqzFbbedQaQknmZxz/cyd9GhPB4hLfWcYQQon52nYmqRXRsSbd2znwi3WdCiHpGCk09odPpGGPy5cCZfA6cydM6jhBC1JgUmnrk0fs8ae5gJ0OdhRD1ihSaeqSFYxOGGTxYt/9n8otKtI4jhBA1UmWhcXFxue1XixYt6NatW11mFdxaPuBGiYXVe2VRNCFE/VBloenSpQsFBQVVfl29ehUnJ6e6zCqAIE9XDN5uLN99ChkwKISoD6osNKtXr6725JocI2rfGJMvJ3IL2XnyktZRhBCiWlUWms6dO1d7ck2OEbXvwZAOuDZrwvJdsnyAEML2VTsYYNeuXURERODs7IyDgwN2dna4uLjURTZRBccmdowM92LTofPkXL2hdRwhhLitagvNpEmTWLlyJX5+fly/fp2FCxfym9/8pi6yidsYbfLFbFF8nlpxGQUhhLAlNRre3LVrV0pLS7Gzs2PcuHFs3LjR2rlENTq1ceJXXduwMiWbUosMChBC2K5qC03z5s0pLi7GYDAQFxdHfHx8pcsri7o3xuTD2bzrbD2So3UUIYSoUrWFZtmyZZSWlvLuu+/i5OREdna2jDazEQP929HOpanMfyaEsGnVLhPg6+sLQLNmzZg5c6bVA4mas7fT82SED29vOUb25SK8WzXXOpIQQlRQZaEJDg6+7VouBw4csEogcWdGRfrw7tbjLN99mhlDemgdRwhRjyil6mTNrioLzddff10WZOjQoWzYsMHqYcSda+/qyEB/dz7fk83/DvKjqb0siiaEqN5NcykvfLKXR+/z5KFQD6s+V5X3aHx9ffH19aVjx440bdq07Ptfvu7FqlWrCAwMRK/Xl1s0p7i4mHHjxhEcHExoaCjJycmVnp+eno7JZMJgMGA0GklJSSm3PzU1FXt7e7744ot7yllfjO7py+XCYjb+eF7rKEKIesBiUUz9fD9bjuRw02z9wV2azN4cFBTEmjVr6NOnT7ntCxYsAODgwYMkJSUxderUSke4xcXFMXPmTNLT05k1axZxcXFl+0pLS3n11Vd54IEHrNsIG/Krrm3wbd1clg8QQlRLKcWfvjrE1wfO8dshPXgs3Mvqz1ll19nevXvL/n39+nX27dtXbhLHsLCwu35Sf3//SrdnZGQwYMAAANzd3XFzc2PPnj1ERkaWO06n01FQUABAfn4+Hh7/f9n3zjvvMGLECFJTU+86X32j1+sY3dOHv244wpHzBfRoLzM3CCEq9+6W4yzZeYrnenfi+b5d6uQ5qyw0U6dOLft3+/bteeWVV8q+1+l0bNmypdbDhIaGsm7dOkaNGkV2djZpaWlkZ2dXKDTz5s0jOjqaadOmYbFY+OGHHwA4e/YsiYmJbN26tVEVGoCR4d7M/eYoy3ed5s+PBGkdRwhhg1bsPs3fk44y/D5Pfjuk8j/4raHKQrN169Z7euCBAwdy/nzFewazZ89m2LBhlZ4TExPD4cOHMRqN+Pr6EhUVhZ1dxZvb8+fPJz4+nhEjRvD5558zfvx4Nm/ezMsvv8ybb76JXl99j2BCQgIJCQkA5Obm3mHrbE9LJwceDO5A4r6zzBjSA6em1Y5cF0I0Iht/PMfvvzxI/+5tefOxEPR66482K6OqkJaWVtWuOzrmdvr27atSU1Or3N+rVy916NChCttdXFyUxWJRSillsVhUixYtlFJKdezYUfn6+ipfX1/l5OSk2rZtqxITE6vNER4efpctsC17si4r31e/Vst3ndI6ihDChuw8cVH5/W6DeuS97arwZkmtPW5N3zur/NN/3LhxXLlyhcuXL1f5NX78+FotekVFRRQWFgKQlJSEvb09AQEBFY7z8PBg27ZtAGzZsgU/Pz8AMjMzycrKIisri8cee4z333+fRx55pFYz2rIwHzf8O7jwyS5ZFE0Iccuhn/N5bskefFo156NnI2juUPe9HVU+Y35+PuHh4bd9w2rbtu1dPWliYiKTJ08mNzeXoUOHYjAY2LRpEzk5OURHR6PX6/H09GTZsmVl58TGxjJx4kSMRiMLFixgypQpmM1mHB0dy7rAGjudTscYkw+/S/yRfdl5hPm01DqSEEJDpy8V8exHqTg72rM0JpKWTg6a5NAp+dMXo9FY7vM89VnhTTM9//otDwS24x+PG7SOI4TQSO7Vmzz2wQ/kXy/hi4m96Oreotafo6bvnZp8jkZYj1NTex69z5OvD5zjSmGx1nGEEBq4eqOEsR+nkFNwk4/GRlilyNwJKTQN0BiTL8VmC1+kndE6ihCijt00l/L8sjR+On+V98eE2UQXuhSaBqh7+xZEdGzJ8t2nsMiiaEI0GqUWxf9+ls4PJy7x1sgQ+nd31zoSUINCM3z4cNavXy+LndUzY0y+ZF0qYseJi1pHEULUAaUUM9f9yIaD5/n9UH8evc/6U8vUVLWF5sUXX2TFihX4+fkxY8YMfvrpp7rIJe7R4KD2tHJykPnPhGgk/vntMT7ZdZrn+3YmtndnreOUU22hGThwIMuXL2fv3r107NiRgQMHEhUVxccff0xJSUldZBR3oam9HSONXmw+nMP5/BtaxxFCWNGyXaeYt/kYj4V7MWOw7a1LVaN7NJcuXWLx4sUsXLiQ++67jylTprB3714GDRpk7XziHoyO9MWiFCtTTmsdRQhhJesPnOMPa3/k1z3ceWP47Res1Eq1hebRRx+ld+/eFBUV8dVXX7Fu3TqeeOIJ3nnnHa5du1YXGcVd8mndnD5+bfk09TQlpXKPTYiG5ofjF/nfz9IJ92nJu0+FYW9nm+O7qp2L4KWXXqJ///6V7msoH3JsyMaYfHlu6R6+PZzD4KD2WscRQtSSH8/m89zSPXRq48SiZyNo5mC7q+tWW2iuXLnCmjVrym1zdXUlODgYd3fbGDonqjaghzsero4s331KCo0QDUTWxULGfpyCW3MHlsRE4tq8idaRbqvaQrNo0SJ27txZdlWTnJxMeHg4mZmZ/OEPf+Dpp5+2ekhx9+z0OkZF+vD3pKNkXiykUxsnrSMJIe5BTsENnv5oNxYFS8dH0t7VUetI1aq2Q6+kpITDhw+zevVqVq9eTUZGBjqdjt27d/Pmm2/WRUZxj56I9MZer2PFbhnqLER9VnCjhGc/TuXStWI+HhtBl7bOWkeqkWoLzZkzZ2jXrl3Z9+7u7mRnZ9OqVSuaNLHtyzVxi3sLR6ID27Mq7Qw3Skq1jiOEuAs3Skp5bskejudc5YMx4YR6u2kdqcaq7Trr168fDz74ICNHjgRg9erV9OvXj8LCQtzc6k9DG7vRJh/WHzzH+gPnGBFuO58YFkJUr9SimPLpPnZnXuafTxro0+3ulmjRSrWF5r333mPNmjVs374dgGeeeYYRI0ag0+nueblnUXd6dW5Nl7ZOfLL7lBQaIeoRpRS///Igmw5dYOZDAQwzeGod6Y7dttCUlpYycOBAtm7dyogRI+oqk7ACnU7H6J6+zPo6g0M/5xPo4ap1JCFEDfwj6SgrU7L5Tf8ujLu/k9Zx7spt79HY2dmh1+vJz8+v1SddtWoVgYGB6PX6cp/FKS4uZty4cQQHBxMaGkpycnKl56enp2MymTAYDBiNRlJSUoBbI+JcXV0xGAwYDAZmzZpVq7nruxFhXjg20fPJLpkpQIj6YPGOTN7ZcpwnI7yZ9kB3rePctWq7zpydnQkODmbQoEE4Of3/0Ni33377rp80KCiINWvW8Pzzz5fbvmDBAgAOHjxITk4OQ4YMITU1Fb2+fD2Mi4tj5syZDBkyhA0bNhAXF1dWlHr37s3XX39919kaMtfmTXgoxIO16Wd57X960MJRBnMIYau+2v8zf/o6gwcC2vGXR4JscmqZmqq20AwfPpzhw4fX6pP6+/tXuj0jI4MBAwYAt0a3ubm5sWfPHiIjI8sdp9PpKCgoACA/Px8PD49azdeQjTH5sirtDIn7zvJMr45axxFCVOL7Y7m88nk6Eb6teHvUfTY7tUxNVVtonn32Wa5fv87p06fp3t26l26hoaGsW7eOUaNGkZ2dTVpaGtnZ2RUKzbx584iOjmbatGlYLBZ++OGHsn07d+4kNDQUDw8P5s6dS2BgoFUz1zeh3m4Ee7ryya5TPG3yrdd/JQnREB04k8fzy9Lo0taZBc8acWxiu1PL1FS1ZfKrr77CYDAwePBg4Nb9kYcffrjaBx44cCBBQUEVvtauXVvlOTExMXh5eWE0Gnn55ZeJiorCzq7iD3n+/PnEx8eTnZ1NfHw848ePByAsLIxTp06xf/9+Jk+ezCOPPFLlcyUkJGA0GjEajeTm5lbbnoZkjMmHoxeukZp1ResoQoj/cDL3GmM/TqWVkwNLYyJxbdZAurdVNcLCwlReXp4yGAxl2wIDA6s7rUb69u2rUlNTq9zfq1cvdejQoQrbXVxclMViUUopZbFYVIsWLSo939fXV+Xm5labIzw8vIaJG4bCmyUqaOZG9dLKvVpHEUL82/n86ypqzrcqbNY36mTuNa3j1EhN3zurvaJp0qQJrq7lh8L+98352lJUVERhYSEASUlJ2NvbExAQUOE4Dw8Ptm3bBsCWLVvw8/MD4Pz58yilAEhJScFisdC6dWurZK3PmjvYMyLMi38dPM+laze1jiNEo5d/vYRnP0ohr6iYxeMiG9ychNXeowkMDGTFihWUlpZy7Ngx3n77baKiou7pSRMTE5k8eTK5ubkMHToUg8HApk2byMnJITo6Gr1ej6enJ8uWLSs7JzY2lokTJ2I0GlmwYAFTpkzBbDbj6OhIQkICAF988QXz58/H3t6eZs2a8emnn8o9iCqMMfmw+IcsPt9zhhf6ddE6jhCN1o2SUmKXpHIi9xofj40k2KvhfcZNp365BKhCUVERs2fP5ptvvkEpRXR0NK+//jqOjrY/Y2hNGY3GRrm2zpMJOzmbd51t0/qj10tBFqKumUstTPxkL98eucA7o+7jwZD6NYK2pu+d1V7RNG/enNmzZzN79uxaCSZsxxiTL5NW7GPbsVz6d5e1hYSoS0opXks8yObDF5g1LLDeFZk7UW2hOXr0KHPnziUrKwuz2Vy2fcuWLVYNJqzvgYD2tHFuyvJdp6TQCFHH3tr0E5/vOcNLA7o2+M+0VVtoRo4cycSJE4mNja10qLGovxzs9TwZ4c37ycc5m3cdT7dmWkcSolH4aHsm7yefYFSkD/87qJvWcayu2kJjb2/PCy+8UBdZhAZG9fTh/eTjrNx9mmnR9XcuJSHqiy/3nWXW1xkMDmxf76eWqalqxyk/9NBDvP/++5w7d47Lly+XfYmGwdOtGQN6uPNpajbFZovWcYRo0JJ/ymHaqv2YOrdi3pMG7BrJIJxqr2iWLFkCwFtvvVW2TafTcfLkSeulEnVqtMmXzYdT+SbjfIO+ISmElvadvsILn+zFr10LEp5pGFPL1FS1hSYzM7MucggN9fFri1fLZny8I4shQR0azV9ZQtSV4znXiFmcStsWTVkSE4FLI5s5vcqus7/97W9l/161alW5fa+99pr1Eok6Z6fXMXlAV9JOXeFvG49oHUeIBuVc/nWeWbQbO72OZeMjcW/RcD6DWFNVFppPP/207N9z5swpt2/jxo3WSyQ08USED8/08uXD707yaYosjCZEbcgrKuaZRSkU3DCzeFwkvq0b1tQyNVVl19l/Thjw35MHVDOZgKin/vBgAKcuFfH7L3/Ep1Vzorq20TqSEPXW9eJSxi/Zw6lLRSyOiSDIs+FNLVNTVV7R/OeQu/8eftcYhuM1RvZ2et556j46t3Vi4idpnMi9pnUkIeqlklILv1mxl72nrzDvSQNRXRr3H21VFpr9+/fj4uJCixYtOHDgAC4uLmXfHzx4sC4zijrk4tiERc9G0MROT8ziVK4UFmsdSYh6RSnFjNUH2XIkhz8PC+J/gjtoHUlzVRaa0tJSCgoKuHr1KmazmYKCgrLvS0pK6jKjqGPerZqT8IyRc/k3eP6TNG6aS7WOJES98ca/jrB67xleHujHGJOv1nFsQv1eiFpYTbhvS956LISUzMu8tuZHuS8nRA0s+O4kH353kqdNvkz5tZ/WcWxGtZ+jEY3XMIMnmRcLmbf5GJ3bOvGb/l21jiSEzVqddobZGw7zP8Ht+ePDgXIv+z9IoRG3NeXXfmReLOStTT/RqY2T9DcLUYmtR3KIW32AqC6tiX+i8UwtU1OadJ2tWrWKwMBA9Hp9uUVziouLGTduHMHBwYSGhpKcnFzp+enp6ZhMJgwGA0ajkZSUlLJ9ycnJGAwGAgMD6du3r7Wb0uDpdDreHBFCuG9L/vezdPZn52kdSQibknbqCi8sT8O/Qws+fDqcpvaNZ2qZmtKk0AQFBbFmzRr69OlTbvuCBQsAOHjwIElJSUydOhWLpeJEj3FxccycOZP09HRmzZpFXFwcAHl5ebz44ousW7eOQ4cOVZjRQNwdxyZ2fPh0OG1bNCV26R7O5l3XOpIQNuHYhavELE6lnYsjH4+NpEUjm1qmpjQpNP7+/nTvXnFK+oyMDAYMGACAu7s7bm5ulS4TqtPpKCgoACA/Px8Pj1sTQa5YsYJTIKdCAAAgAElEQVThw4fj4+NT9hiidrRxbsrHYyO4UVzK+MWpXLtprv4kIRqwn/Ou88xHKTjY61kW05O2LZpqHclm2dSos9DQUNatW4fZbCYzM5O0tDSys7MrHDdv3jymT5+Ot7c306ZNK5si5+jRo1y5coV+/foRHh7O0qVL67oJDZpfuxa8NzqMYznXmLJyH6UWGYkmGqcrhcU8vWg3126YWTIuEp/WzbWOZNOsNhhg4MCBnD9/vsL22bNnM2zYsErPiYmJ4fDhwxiNRnx9fYmKiqp0Vc/58+cTHx/PiBEj+Pzzzxk/fjybN2/GbDaTlpbGt99+y/Xr1+nVqxcmk4lu3SquYJeQkEBCQgIAubm599jaxqNPt7b88eFAXv/yR/664TCvPxigdSQh6lRRsZlxi1PJvnKdpTGRBHi4aB3J5lmt0GzevPmOz7G3tyc+Pr7s+6ioqEqLxJIlS/jnP/8J3FpqOjY2FgAvLy9at26Nk5MTTk5O9OnTh/3791f6GBMmTGDChAkAGI3GO87amD1t8uVEzjUWbc+kc1snRveUD6WJxuHStZu8sHwvB87k8f7ocEydW2sdqV6wqa6zoqIiCgsLAUhKSsLe3p6AgIp/MXt4eLBt2zYAtmzZgp/frQ9GDRs2jO3bt2M2mykqKmL37t34+/vXXQMakdcfDKB/97b8Ye0hvj8mV4Si4fvxbD4Pv7uD/dl5xD9hYHBQe60j1RuaFJrExES8vLzYuXMnQ4cOJTo6GoCcnBzCwsLw9/fnzTffZNmyZWXnxMbGlg0MWLBgAVOnTiU0NJTXXnutrAvM39+fwYMHExISQmRkJLGxsQQFBdV9AxsBO72Od54Kw8/dmReX7+V4zlWtIwlhNYn7zjBi/g8opVj9QhTDDJ5aR6pXdErmFsFoNFY6uk1U72zedYa9u4NmDnq+fPF+WjvLyBvRcJhLLcz51xEWbc+kZ6dWvDc6jDbyO16mpu+dNtV1JuofT7dmLHzWSE7BTZ5fJhNwiobjcmExz3yUwqLtmYyN6sgnsT2lyNwlKTTinhm83fjH4wb2nLrCjNUHZQJOUe/9eDafh97Zzp5TV5g7MpQ/PhxIEzt5u7xbMteZqBVDQzqQebEbc785Sqc2TrwkM9eKempt+lleXX2Als0dWPV8L0K93bSOVO9JoRG15jf9u3LyYiH/SLpVbB4K9dA6khA1Zi618ObGIyz4PpPIjrfux8in/WuHFBpRa3Q6HXOGB3Pm8nWmrtqPZ8tmhPm01DqWENW6UljM5JX72H78Is/08uX3QwNwsJeustoiP0lRq5ra2/HB0+G0d3FkwtI9nLlSpHUkIW4r4+cCHn5vOymZl/nbiBBmDQuSIlPL5Kcpal0rJwc+GhvBTbOF8Yv3cPWGLP0tbNPXB35mxPwfKDZb+Ox5E49HeGsdqUGSQiOsoqu7M/NHh3M89xqTV+7DXFpxuQchtFJqUbzxryNMWrGPAA8Xvpr8K+6Tbl6rkUIjrOZXfm3487Agkn/K5S/rD2sdRwgA8oqKGbc4lQ+2nWB0Tx9WPmfCvYWj1rEaNBkMIKzqqZ4+nMy9xsJ/T8D5TK+OWkcSjdiR8wVMWJrGufzrzBkezKhIH60jNQpSaITV/fZ//Mm6VMQf1x3Cp1Vz+nWXBelE3dtw8BzTVu3Hqak9n04wEe7bSutIjYZ0nQmrs9Pr+OeTBnq0d2HSin38dF4m4BR1p9SieGvTEV5cvpfu7Vvw9eRfSZGpY1JoRJ1wamrPorFGmjvYEbM4ldyrN7WOJBqB/OsljF+SyntbT/BkhDefTjDRzkXux9Q1KTSiznRwbcaiZyO4VHiTCcv2cKNEJuAU1nPswlUeeW8H249d5C+PBDFneDBN7Suu2CusTwqNqFPBXq7Me8LAvtN5TP/igEzAKaxi44/neeS9HVy9YWblBBNjTL7odDqtYzVaUmhEnRsc1IFXB/fgq/0/M2/zMa3jiAbEYlH845ufmPhJGl3bteCryfcT0VHux2hNk0KzatUqAgMD0ev15RbNKS4uZty4cQQHBxMaGkpycnKl56enp2MymTAYDBiNRlJSUgB46623MBgMGAwGgoKCsLOz4/Lly3XRJHGHJvbtzMhwL/757THWpp/VOo5oAApulPDc0j28veU4I8O9+GyCiQ6uzbSOJQCUBjIyMtSRI0dU3759VWpqatn2d999V40dO1YppdSFCxdUWFiYKi0trXD+oEGD1IYNG5RSSq1fv1717du3wjHr1q1T/fv3r1Ge8PDwu2iFuFc3S0rV4x/8oPxe26D2ZF3SOo6ox45duKr6z92quvx2vVryQ6ayWCxaR2oUavreqckVjb+/P927d6+wPSMjgwEDBgDg7u6Om5tbpcuE6nQ6CgoKAMjPz8fDo+J09CtXrmTUqFG1nFzUJgd7PR+MCcezZTMmLE0j+7JMwCnuXFLGBR55bwf5RSUsj+3JM706yv0YG2NT92hCQ0NZt24dZrOZzMxM0tLSyM7OrnDcvHnzmD59Ot7e3kybNo05c+aU219UVMTGjRsZMWJEXUUXd6mlkwOLnjVitihiFqdSIBNwihqyWBTzNh/luaV76NTGia8m/4qenVtrHUtUwmqFZuDAgQQFBVX4Wrt2bZXnxMTE4OXlhdFo5OWXXyYqKgo7u4rDEefPn098fDzZ2dnEx8czfvz4cvu/+uor7r//flq1qvomYEJCAkajEaPRSG5u7t03VNyzzm2dmT8mjMyLhfxm+V6ZgFNU6+qNEp7/JI15m48xPMyTVRN74eEm92NslU4p7caX9uvXj7lz52I0GivdHxUVxcKFCwkICCi33dXVlby8PHQ6HUopXF1dy7rSAB599FFGjhzJU089VaMcRqOx0i46Ubc+T80mbvUBxph8+POwIOn+EJU6mXuN55buIetSEb/7H3/G3S9dZVqp6XunTXWdFRUVUVhYCEBSUhL29vYVigyAh4cH27ZtA2DLli34+f3/+vT5+fls27aNYcOG1U1oUWsej/Dm+b6d+WTXaT7ekaV1HGGDthy5wLB3d3ClqIRl4yOJ+VUnKTL1gCaTaiYmJjJ58mRyc3MZOnQoBoOBTZs2kZOTQ3R0NHq9Hk9PT5YtW1Z2TmxsLBMnTsRoNLJgwQKmTJmC2WzG0dGRhISEco/9wAMP4OTkpEXTxD16NboHWRcL+cv6DDq2ac6AHu20jiRsgMWieG/rcf6x+SgBHVz48OlwvFo21zqWqCFNu85shXSd2ZaiYjOPf7iTzNxCvnghCv8OLlpHEhq6dtPMtM/3s/HQeR4xeDBneAjNHGQqGVtQL7vOhABo7mDPomcjaOHYhPGLU8kpuKF1JKGRrIuFPPreDr7JOM/vh/oT/4RBikw9JIVG2KR2Lo4sfNbIlaJbn/aWCTgbn+Sfcnj43e3kXrvJ0piexPbuLPdj6ikpNMJmBXm68s8nDRw4m8/Uz/djsTT6Xt5GQSnF+8nHGbc4Fc+Wzflq0q/4lV8brWOJeyCFRti0BwLb89shPVh/8Bz/SDqqdRxhZYU3zUxasY+/bfyJocEdWP1CL7xbyU3/+k6WchY277nenTmZW8i7W4/TqY0TI8K9tI4krOD0pSImLNvD0QtX+e2QHkzoI11lDYUUGmHzdDodf34kiNOXi5ix5gDerZoT2Ummfm9Ivjuay+SV+wBYPC6SPt3aapxI1CbpOhP1QhM7PfNHh+PdqjnPL9tD1sVCrSOJWqCU4sNtJxj7cQodXB1ZN+l+KTINkBQaUW+4Nm/CR89GoICYJankF8kEnPVZUbGZlz5NZ86/jjA4qD2rX4jCt7V80LohkkIj6pWObZz4cEw42ZeLeGF5GiUyAWe9lH25iBHzd/L1gZ+JG9yd954Kw6mp9OQ3VFJoRL3Ts3Nr3hgewg8nLvGHtT8ik1vULzuOX+Shd7dz9koRH4+N4MV+XeWmfwMnf0KIemlEuBcnL17jva0n6NzGmef6dNY6kqiGUopF2zP564bDdHV3JuFpIx3bSFdZYyCFRtRbUwd1J/NiIX/912F8WzfngcD2WkcSVcgvKmHmuh/5Mv1nBge2Z+7joThLV1mjIa+0qLf0eh1/H2ng7JWdTPk0nVUTexHk6ap1LPFvxWYLW3/KIXHvWbYcyaHEYmHaA934TX/pKmtspNCIeq2Zgx0LnjHyyHs7iF2yh7WT7qedi6PWsRotpRT7svNI3HuWrw/8zJWiEto4OzDa5MPIcG8CPGQm7sZICo2o99xdHFn4bAQjP/iB2CV7+Ox5E80d5Fe7Lp2+VETivrN8mX6WzIuFNLXX80Bge4bf50lvvzbY28m4o8ZM/jeKBiHAw4W3R93Hc0v38Mpn+3l/dBh6vXTPWFN+UQlfH/yZxL1n2XPqCgC9OrfmhX5dGBLUnhaOTTROKGyFJn9mrFq1isDAQPR6fblFc4qLixk3bhzBwcGEhoaSnJxc6fnp6emYTCYMBgNGo5GUlBTg1jLODz30EKGhoQQGBvLxxx/XRXOEjfi1fzt+NzSAjYfO82TCLhZ8d5Kfzl+V4c+1qNhs4ZtD53nhkzQiZm/md4k/kne9hOnR3dkxYwArJ5h43OgtRUaUo8kVTVBQEGvWrOH5558vt33BggUAHDx4kJycHIYMGUJqaip6ffl6GBcXx8yZMxkyZAgbNmwgLi6O5ORk3nvvPQICAvjqq6/Izc2le/fujB49GgcHhzprm9BWzP0dKbVY+HzPGWZvOMzsDYdp7+JIb7829OnWll91bUNLJ/l9uBO3u+8y/D4vgjxd5Oa+uC1NCo2/v3+l2zMyMhgwYAAA7u7uuLm5sWfPHiIjI8sdp9PpKCgoAG5dxXh4eJRtv3r11l+w165do1WrVtjbS+9gY6LT6ZjQpwsT+nTh57zrfH8sl++OXuSbjAusSjuDTgchnq706daWPt3aYvB2o4ncP6jU7e67/MqvjfzcRI3Z1LtwaGgo69atY9SoUWRnZ5OWlkZ2dnaFQjNv3jyio6OZNm0aFouFH374AYBJkybx8MMP4+HhwdWrV/nss88qXA2JxsPDrRlPRPjwRIQPpRbF/jN5fHc0l++O5vLe1uO8s+U4LZraE9W1Nb392tK3W9tGv/ZJflEJ6w+eI3HfGVKzbt13MXVuJfddxD2xWqEZOHAg58+fr7B99uzZDBs2rNJzYmJiOHz4MEajEV9fX6KiorCzq7g++Pz584mPj2fEiBF8/vnnjB8/ns2bN7Np0yYMBgNbtmzhxIkTDBo0iN69e+PiUnFIZUJCAgkJCQDk5ubeY2uFrbPT6wjzaUmYT0teHtiN/Osl/HD8It/9+4pn06ELAHRq40Sff3ezmTq3bhTzbxWbLST/lEPivrN8eziH4lILXd2dmR7dnUfu88TTrZnWEUU9p1Ma3int168fc+fOxWg0Vro/KiqKhQsXEhAQUG67q6sreXl56HQ6lFK4urpSUFDA0KFDmTFjBr179wZgwIABvPHGGxWuiP6b0WgsNyhBNC5KKU7kFv67my2XXScvc72klCZ2OsJ9W97qZvNrS0AHlwYzkk0pRXp2Hmv+677LQ6Eect9F1FhN3ztt6s+1oqIilFI4OTmRlJSEvb19hSID4OHhwbZt2+jXrx9btmzBz88PAB8fH7799lt69+7NhQsX+Omnn+jcWebAEren0+no6u5MV3dnxt3fiZvmUvZkXeG7o7lsO5rL3zb+xN82/kQbZwd6+7Wlt18bevu1pW2LplpHv2PZl2/dd0nc9//3XQYFtGN4mCe9/drKfRdhFZpc0SQmJjJ58mRyc3Nxc3PDYDCwadMmsrKyiI6ORq/X4+npyaJFi/D19QUgNjaWiRMnYjQa2b59O1OmTMFsNuPo6Mj7779PeHg4P//8M2PHjuXcuXMopZgxYwZjxoypNo9c0YjbySm4wffHbnWzfX/sIpcLiwEI6ODy70EFbQj3bUlT+4rdvLagqvsuw+/zYnBwe1zkvou4SzV979S068xWSKERNWWxKA79XMB3x25d7ew9dQWzRdHcwQ5T59Zl93c6tXHStOup2Gxh29FcEvedYXPGrfsuXdo6MTzMS+67iFpTL7vOhLB1er2OYC9Xgr1c+U3/rly7aWbniUu3RrMdy2XLkRwAvFo2+/dItjZEdW1TJ1cNv9x3Sdx3lq/237rv0trJgad6+jA8zJNgT1e57yI0IYVGiHvg3NSeQQHtGBTQDoBTlwr57thFvjuay1f7f2Zlymns9Dru83ajT7db93dCvNywq8VBBb/cd/ly31lOyn0XYYOk6wzpOhPWUVJqYe+pK2VDqH/8OR+lwK15E+7v2oa+fm3p3a0NHVzvvBsr/3oJGw6eY81eue8itCNdZ0JorImdnp6dW9Ozc2umR8OlazfZfvwi3x29NbBg/YFzAHRr50xvv1szFfTs1ArHJpUPKih33+VwDsXmW/ddpkd3Z5jBA6+WjfvDpsJ2SaERoo60dm7KMIMnwwyeKKU4cv4q3x29NZJt2c5TLNqeSVN7PZGdWtG3W1t6+7WlWzvnyu+7RMp9F1F/SKERQgM6nQ7/Di74d3Dh+b5duF5cyq7MS2VT5Pxl/WHgME4OdhQWl+Jgr+cBue8i6ikpNELYgGYOdvTv7k7/7u4AnM27zvdHc9l3Oo8wXzeGBHeQ+y6i3pJCI4QN8nRrxpORPjwZ6aN1FCHumVx/CyGEsCopNEIIIaxKCo0QQgirkkIjhBDCqqTQCCGEsCopNEIIIaxKCo0QQgirkkIjhBDCqmT2ZqBNmzZ07Njxrs7Nzc2lbdu2tRtII9IW29RQ2tJQ2gHSll9kZWVx8eLFao+TQnOPGtISA9IW29RQ2tJQ2gHSljslXWdCCCGsSgqNEEIIq7L74x//+EetQ9R34eHhWkeoNdIW29RQ2tJQ2gHSljsh92iEEEJYlXSdCSGEsCopNLexceNGunfvTteuXXnjjTcq7D916hS//vWvCQkJoV+/fpw5c6Zsn52dHQaDAYPBwMMPP1yXsSuIiYnB3d2doKCgSvcrpXjppZfo2rUrISEh7N27t2zfkiVL8PPzw8/PjyVLltRV5CrdS1ts6TWB6tty5MgRevXqRdOmTZk7d265fdX9btale2lHx44dCQ4OxmAwYDQa6yLubVXXluXLlxMSEkJwcDBRUVHs37+/bJ8tvSZwb22p9ddFiUqZzWbVuXNndeLECXXz5k0VEhKiDh06VO6Yxx57TC1evFgppdS3336rxowZU7bPycmpTvPezrZt21RaWpoKDAysdP/69evV4MGDlcViUTt37lSRkZFKKaUuXbqkOnXqpC5duqQuX76sOnXqpC5fvlyX0Su427YoZVuviVLVt+XChQsqJSVFvfbaa+qtt94q216T3826dLftUEopX19flZubWxcxa6S6tuzYsaPs/8CGDRvKfr9s7TVR6u7bolTtvy5yRVOFlJQUunbtSufOnXFwcODJJ59k7dq15Y7JyMhgwIABAPTv37/CflvRp08fWrVqVeX+tWvX8swzz6DT6TCZTOTl5XHu3Dk2bdrEoEGDaNWqFS1btmTQoEFs3LixDpNXdLdtsUXVtcXd3Z2IiAiaNCm/hHNNfjfr0t22wxZV15aoqChatmwJgMlkKuvFsLXXBO6+LdYghaYKZ8+exdvbu+x7Ly8vzp49W+6Y0NBQ1qxZA0BiYiJXr17l0qVLANy4cQOj0YjJZOLLL7+su+B3oaq21uRnYGtul7k+vSa3Ux9fl6rodDoeeOABwsPDSUhI0DrOHVm0aBFDhgwB6v9r8p9tgdp/Xezv+REasblz5zJp0iQWL15Mnz598PT0xM7ODrh1/8bT05OTJ08yYMAAgoOD6dKli8aJGzd5TWzP9u3b8fT0JCcnh0GDBtGjRw/69Omjdaxqbd26lUWLFrF9+3ato9yzytpS26+LXNFUwdPTk+zs7LLvz5w5g6enZ7ljPDw8WLNmDfv27WP27NkAuLm5lZ0P0LlzZ/r168e+ffvqKPmdq6qtNfkZ2JrbZa5Pr8nt1MfXpSq/5HZ3d+fRRx8lJSVF40TVO3DgALGxsaxdu5bWrVsD9fc1qawtUPuvixSaKkRERHDs2DEyMzMpLi7m008/rTBS6eLFi1gsFgDmzJlDTEwMAFeuXOHmzZtlx+zYsYOAgIC6bcAdePjhh1m6dClKKXbt2oWrqysdOnQgOjqab775hitXrnDlyhW++eYboqOjtY57W1W1pb69JrdTk9/N+qCwsJCrV6+W/fubb76pcoSUrTh9+jTDhw9n2bJldOvWrWx7fXxNqmqLVV6XWhtW0ACtX79e+fn5qc6dO6u//OUvSimlXn/9dbV27VqllFKrVq1SXbt2VX5+fmr8+PHqxo0bSqlbozmCgoJUSEiICgoKUgsXLtSsDUop9eSTT6r27dsre3t75enpqRYuXKjmz5+v5s+fr5RSymKxqBdffFF17txZBQUFqdTU1LJzFy1apLp06aK6dOmiPvroI62aUOZu22Jrr4lS1bfl3LlzytPTU7Vo0UK5uroqT09PlZ+fr5Sq/HdTK3fbjhMnTqiQkBAVEhKiAgICNG+HUtW3Zfz48crNzU2Fhoaq0NBQFR4eXnauLb0mSt19W6zxusjMAEIIIaxKus6EEEJYlRQaIYQQViWFRgghhFVJoRFCCGFVUmiEEEJYlRQa0WD8MjtzYGAgoaGh/P3vfy/7nNOePXt46aWXqjw3KyuLFStW1FXUCs/drFkzDAZDnT3nX//61zp7rl+cOHECg8GAs7NznT+30JYMbxYNhrOzM9euXQMgJyeHp556ivvvv58//elP1Z6bnJzM3Llz+frrr60ds4KsrCwefPBBfvzxxzp7zv/8Wf0npRRKKfR66/0NWtVzi4ZLrmhEg+Tu7k5CQgLvvvsuSimSk5N58MEHAdi2bVvZujT33XcfV69eZcaMGXz//fcYDAbi4+PJysqid+/ehIWFERYWxg8//ADcKkj9+vXjscceo0ePHowePZpf/lZLTU0lKiqK0NBQIiMjuXr1KqWlpUyfPp2IiAhCQkL48MMPa5T/kUceITw8nMDAwHKTGjo7O/O73/2O0NBQTCYTFy5cAG5dLZhMJoKDg/n9739fdtVw7tw5+vTpg8FgICgoiO+//54ZM2Zw/fp1DAYDo0ePJisri+7du/PMM88QFBREdnY2K1euJDg4mKCgIF599dWy5//mm2/o1asXYWFhjBw5sqxgdOzYkd/+9rdl65fs3buX6OhounTpwgcffHCPr6ao9+75I59C2IjK1ptxdXVV58+fV1u3blVDhw5VSin14IMPqu3btyullLp69aoqKSkpt18ppQoLC9X169eVUkodPXq07FPTW7duVS4uLio7O1uVlpYqk8mkvv/+e3Xz5k3VqVMnlZKSopRSKj8/X5WUlKgPP/xQ/fnPf1ZKKXXjxg0VHh6uTp48WS5jZmZmhTVDLl26pJRSqqioSAUGBqqLFy8qpZQC1Lp165RSSk2fPr3ssYcOHapWrFihlFJq/vz5ZT+LuXPnln2y22w2q4KCggo/q8zMTKXT6dTOnTuVUkqdPXtWeXt7q5ycHFVSUqL69++vEhMTVW5ururdu7e6du2aUkqpN954Q/3pT39SSt1av+T9999XSin18ssvq+DgYFVQUKBycnKUu7t7ta+TaNhk9mbR6Nx///288sorjB49muHDh+Pl5VXhmJKSEiZNmkR6ejp2dnYcPXq0bF9kZGTZOQaDgaysrLI51SIiIgBwcXEBbl0BHDhwgC+++AKA/Px8jh07RqdOnW6b8e233yYxMRGA7Oxsjh07RuvWrXFwcCi7MgsPDycpKQmAnTt3li198NRTTzFt2jTg1hxcMTExlJSU8Mgjj1R5H8jX1xeTyQTcujLr168fbdu2BWD06NF899132Nvbk5GRwf333w9AcXExvXr1KnuMX+b2Cg4O5tq1a7Ro0YIWLVrQtGlT8vLyyiacFY2PFBrRYJ08eRI7Ozvc3d05fPhw2fYZM2YwdOhQNmzYwP3338+mTZsqnBsfH0+7du3Yv38/FosFR0fHsn1NmzYt+7ednR1ms7nKDEop3nnnnTuajDQ5OZnNmzezc+dOmjdvTr9+/bhx4wYATZo0QafT1ei54dbiV9999x3r169n7NixvPLKKzzzzDMVjnNycqo2l1KKQYMGsXLlykr3//Jz0ev15X5Ger2+2pyiYZN7NKJBys3NZeLEiUyaNKnsjfkXJ06cIDg4mFdffZWIiAiOHDlCixYtymashVtXHh06dECv17Ns2TJKS0tv+3zdu3fn3LlzpKamAnD16lXMZjPR0dHMnz+fkpISAI4ePUphYeFtHys/P5+WLVvSvHlzjhw5wq5du6ptr8lkYvXq1QB8+umnZdtPnTpFu3bteO6554iNjWXv3r3ArYL1S6b/FhkZybZt27h48SKlpaWsXLmSvn37YjKZ2LFjB8ePHwduzez7n1d6QlRFrmhEg/HLDe6SkhLs7e15+umneeWVVyocN2/ePLZu3YperycwMJAhQ4ag1+uxs7MjNDSUsWPH8uKLLzJixAiWLl3K4MGDq/2L38HBgc8++4zJkydz/fp1mjVrxubNm4mNjSUrK4uwsDCUUrRt27ba1T0HDx7MBx98gL+/P927dy/r0rqdefPmMWbMGGbPns3gwYNxdXUFbl0dvfXWWzRp0gRnZ2eWLl0KwIQJEwgJCSEsLKxsLaVfdOjQgTfeeIP+/fujlGLo0KEMGzYMgMWLFzNq1KiyJRf+8pe/lJtiXojKyPBmITRWG8Obi4qKaNasGTqdjk8//ZSVK1dqvmZ9VWR4c+MjVzRCaMzOzo78/HwMBgPp6el39RhpaWlMmjQJpRRubm589NFHtZzy3p04cYIRI0bQrl07raOIOiZXNEIIIaxKBgMIIYSwKik0QgghrEoKjRBCCKuSQiOEEMKqpOf16NgAAAATSURBVNAIIYSwKik0QgghrOr/AJeFmXBHs9rKAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig=plt.figure()\n",
    "fig.patch.set_facecolor('white')  # useful when dark background used on jupyter lab\n",
    "plt.plot(dc.distances, dc.energies)\n",
    "plt.xlabel(\"Distance [{}]\".format(dc.queue[0].posinp.units))\n",
    "plt.ylabel(\"Energy [Ha]\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results found here are as expected: when performing a geometry optimization of the N$_2$ molecule, you actually get a distance of 1.0935 $\\unicode[serif]{xC5}^4$ with default input parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some errors are raised:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* ### If periodic boundary conditions are used...\n",
    "\n",
    "... because it is not possible to define a distance between two fragments under such boundary conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "Cannot compute a dissociation curve with periodic boundary conditions:\n1   angstroem\nperiodic\nN   0.0   0.0   0.0\n",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-9-4e208da571b1>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mfrag2\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_boundary_conditions\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m\"periodic\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mdistances\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0.95\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1.25\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0.05\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 6\u001b[0;31m \u001b[0mdc3\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mDissociation\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfrag1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfrag2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdistances\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/Documents/Python/MyBigDFT/mybigdft/workflows/dissociation.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, fragment1, fragment2, distances, inputparams, name, run_dir)\u001b[0m\n\u001b[1;32m     69\u001b[0m                 raise ValueError(\n\u001b[1;32m     70\u001b[0m                     \u001b[0;34m\"Cannot compute a dissociation curve with periodic \"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 71\u001b[0;31m                     \"boundary conditions:\\n{}\".format(frag))\n\u001b[0m\u001b[1;32m     72\u001b[0m         \u001b[0;31m# Make sure both fragments use the same units (could actually be\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     73\u001b[0m         \u001b[0;31m# implemented properly in the __add__ method of posinp)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: Cannot compute a dissociation curve with periodic boundary conditions:\n1   angstroem\nperiodic\nN   0.0   0.0   0.0\n"
     ]
    }
   ],
   "source": [
    "from copy import deepcopy\n",
    "frag1 = Posinp([Atom('N', [0.0, 0.0, 0.0])], units=\"angstroem\", boundary_conditions=\"free\")\n",
    "frag2 = deepcopy(frag1)\n",
    "frag2._boundary_conditions = \"periodic\"\n",
    "distances = np.arange(0.95, 1.25, 0.05)\n",
    "dc3 = Dissociation(frag1, frag2, distances)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* ### If units differ...\n",
    "\n",
    "... because the unit conversion is not yet implemented!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "ename": "NotImplementedError",
     "evalue": "Unit conversion of positions needed",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNotImplementedError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-10-072a4efcd7d5>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mfrag2\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_units\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m\"atomic\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mdistances\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0.95\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1.25\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0.05\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 6\u001b[0;31m \u001b[0mdc2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mDissociation\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfrag1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfrag2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdistances\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/Documents/Python/MyBigDFT/mybigdft/workflows/dissociation.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, fragment1, fragment2, distances, inputparams, name, run_dir)\u001b[0m\n\u001b[1;32m     74\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mfragment1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0munits\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0mfragment2\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0munits\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     75\u001b[0m             raise NotImplementedError(\n\u001b[0;32m---> 76\u001b[0;31m                 \"Unit conversion of positions needed\")  # pragma: no cover\n\u001b[0m\u001b[1;32m     77\u001b[0m         \u001b[0;31m# Set the base attributes that are specific to this workflow\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     78\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfragment1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfragment1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNotImplementedError\u001b[0m: Unit conversion of positions needed"
     ]
    }
   ],
   "source": [
    "from copy import deepcopy\n",
    "frag1 = Posinp([Atom('N', [0.0, 0.0, 0.0])], units=\"angstroem\", boundary_conditions=\"free\")\n",
    "frag2 = deepcopy(frag1)\n",
    "frag2._units = \"atomic\"\n",
    "distances = np.arange(0.95, 1.25, 0.05)\n",
    "dc2 = Dissociation(frag1, frag2, distances)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}