{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "02171dee",
   "metadata": {},
   "source": [
    "# M1.D5: Electronic energy of molecules\n",
    "\n",
    "## Learning Objectives\n",
    "\n",
    "* Use Cartesian coordinates to specify the geometry of a diatomic molecule.\n",
    "* Run \"single point\" energy calculations of diatomic molecules using the pyscf package. \n",
    "* Run \"restricted/closed shell\" and \"unrestricted/open shell\" calculations. \n",
    "* Compare the energy of the singlet and triplet electronic state for the diatomic molecules of oxygen and nitrogen \n",
    "* Given bond distances and bond angles, use internal coordinates to specify the geometry of water.\n",
    "* Compare the homolytic and heterolytic bond breaking of H2.\n",
    "\n",
    "# Electronic energy of molecules\n",
    "\n",
    "\n",
    "So far we have solved the following quantum systems:\n",
    "\n",
    "* Particle in box: 1 particle's kinetic energy, no potential. The solution was analytical (a function).\n",
    "* Hydrogen atom: 2 particle's kinetic energy, one attractive potential between the two. The solution albeit more complicated was still analytical.\n",
    "* Polyelectronic atoms: many particle's kinetic energy (independent of each other), atttactive potential and an extra electron-electron repulsion. The solution was not analytical and had to be solved computationally.\n",
    "\n",
    "Now we want to know how to solve molecules or groups of atoms interacting with each other. As usual the hamiltonian is the sum of kinetic and potential energy\n",
    "\n",
    "$$\n",
    "\\hat{H} = \\text{Kinetic Energy} + \\text{Potential Energy}\n",
    "$$\n",
    "\n",
    "For molecules in atomic units:\n",
    "\n",
    "$$\n",
    "\\hat{H} = \\sum_{k=nuclei}-\\frac{1}{2m_k}\\nabla_{k} -\\frac{1}{2}\\sum_{i=elec} \\nabla_i + \\sum_{k=nuclei}\\sum_{i=elec} \\frac{q_{i}q_{k}}{r_{ik}} + \\sum_{i,j=elec} \\frac{q_{i}q_{j}}{r_{ij}} + \\sum_{k,l=nuc}\\frac{q_{k}{q_{l}}}{r_{kl}}\n",
    "$$\n",
    "\n",
    "Translated into English we have:\n",
    "\n",
    "$$\n",
    "\\hat{H} = \\sum\\text{Kinetic of nucleus} + \\sum{\\text{Kinetic of e-}} + \\sum{\\text{attraction of elec-nucl}} + \\sum{\\text{Repulsion elec-elect}} + \\sum{\\text{repulsion between nuclei}}\n",
    "$$\n",
    "\n",
    "## The Born-Oppenheimer approximation and potential energy surface\n",
    "\n",
    "\"Let us define a molecule by a geometry structure determined by the nuclear positions. If the nuclei have fixed positions, the nuclear kinetic term vanishes and the nuclear repulsion term becomes a constant. The Hamiltonian expression of equation has a shortened form that we label as electronic Hamiltonian\"\n",
    "\n",
    "$$\n",
    "\\hat{H}_{elec} = -\\frac{1}{2}\\sum_{i=elec} \\nabla_i + \\sum_{k=nuclei}\\sum_{i=elec} \\frac{q_{i}q_{k}}{r_{ik}} + \\sum_{i,j=elec} \\frac{q_{i}q_{j}}{r_{ij}}\n",
    "$$\n",
    "\n",
    "So we can define the electronic energy $E_{elec}$ of a molecule as the eigenvalues of the electronic Hamiltonian\n",
    "\n",
    "$$\n",
    "\\hat{H}_{elec}\\Psi_{elec} = E_{elec}\\Psi_{elec} \n",
    "$$\n",
    "\n",
    "Without going into much detail it can be shown that the motion of the nuclei is mostly decoupled from the motion of the electrons, so we can imagine that the nuclei are moving on a single electronic state surface, the most common one being the ground electronic state. This surface is called the potential energy surface and the approximation that led us here is called the Born-Oppenheimer approximation.\n",
    "\n",
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/8/81/Potential_Energy_Surface_and_Corresponding_Reaction_Coordinate_Diagram.png\" width=\"500\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8c56a0d",
   "metadata": {},
   "source": [
    "# Single point energy of molecules\n",
    "\n",
    "The calculation below is the \"single point\" energy calculation of H2. \n",
    "The output is saved in the file 'h2_sp.txt'. Open it and try to rationalize the information in it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce724177",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "overwrite output file: m1d5_h2_sp.out\n",
      "-1.0790270257473757\n"
     ]
    }
   ],
   "source": [
    "from pyscf import gto\n",
    "\n",
    "molh2 = gto.M(\n",
    "    atom = 'H 0 0 0; H 0 0 0.53',  # in Angstrom\n",
    "    basis = '6-31g',\n",
    "    symmetry = True,\n",
    "    charge=0, \n",
    "    spin=0,\n",
    "    verbose=4,\n",
    "    output='m1d5_h2_sp.out')  \n",
    "\n",
    "h2 = molh2.RHF().run() #RHF Restricted Hartree-Fock. All electrons are paired\n",
    "\n",
    "print(h2.e_tot)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f1ced51",
   "metadata": {},
   "source": [
    "* Why did we use rhf and not uhf?\n",
    "* The geometry is introduced using cartesian coordinates (xyz). What is the bond distance in this case?\n",
    "\n",
    "Let's calculate the energy of the molecule of water. Notice the alternative way to give the geometry. These are called internal coordinates and they are convenient to introduce bond distances and angles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "56967165",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "overwrite output file: m1d5_h2o_sp.out\n",
      "-76.00900934217299\n"
     ]
    }
   ],
   "source": [
    "#entering the geometry with internal coordinates\n",
    "molh2o = gto.M(\n",
    "    atom = '''\n",
    "    O\n",
    "    H    1  0.96\n",
    "    H    1  0.96 2 104.5\n",
    "    ''',\n",
    "    basis = '6-31g*',\n",
    "    symmetry = True,\n",
    "    charge=0, \n",
    "    spin=0,\n",
    "    verbose=4,\n",
    "    output='m1d5_h2o_sp.out'\n",
    ")\n",
    "h2o = molh2o.RHF().run()\n",
    "print(h2o.e_tot)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c72f451",
   "metadata": {},
   "source": [
    "## Assessing the energy of different spin states \n",
    "\n",
    "We cannot always foresee the most stable spin. If we want to calculate the energy of the molecule of oxygen, by drawing the Lewis structure we would see that all electrons are paired, therefore we would expect that the ground state is a singlet.\n",
    "\n",
    "**Let's compare the energy for the singlet and triplet of N2 and O2**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "4f4793c9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "overwrite output file: m1d5_o2_1.out\n",
      "overwrite output file: m1d5_o2_3.out\n",
      "-149.4533513245776\n",
      "-149.54579333299034\n"
     ]
    }
   ],
   "source": [
    "molo2_singlet = gto.M(\n",
    "    atom = 'O 0 0 0; O 0 0 1.2',  # in Angstrom\n",
    "    basis = '6-31g',\n",
    "    symmetry = True,\n",
    "    charge=0, \n",
    "    spin=0,\n",
    "    verbose=3,\n",
    "    output='m1d5_o2_1.out')  \n",
    "\n",
    "o2_singlet = molo2_singlet.RHF().run() #RHF Restricted Hartree-Fock. All electrons are paired\n",
    "o2_singlet.analyze()\n",
    "\n",
    "molo2_triplet = gto.M(\n",
    "    atom = 'O 0 0 0; O 0 0 1.2',  # in Angstrom\n",
    "    basis = '6-31g',\n",
    "    symmetry = True,\n",
    "    charge=0, \n",
    "    spin=2, #two unpaired electrons\n",
    "    verbose=3,\n",
    "    output='m1d5_o2_3.out')  \n",
    "\n",
    "o2_triplet = molo2_triplet.UHF().run() #UHF Unrestricted Hartree-Fock. Two unpaired electrons\n",
    "o2_triplet.analyze()\n",
    "\n",
    "print(o2_singlet.e_tot)\n",
    "print(o2_triplet.e_tot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "3a5e646b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "overwrite output file: m1d5_n2_1.out\n",
      "overwrite output file: m1d5_n2_3.out\n",
      "-108.86761837305848\n",
      "-108.59455017417186\n"
     ]
    }
   ],
   "source": [
    "moln2_singlet = gto.M(\n",
    "    atom = 'N 0 0 0; N 0 0 1.1',  # in Angstrom\n",
    "    basis = '6-31g',\n",
    "    symmetry = True,\n",
    "    charge=0, \n",
    "    spin=0,\n",
    "    verbose=3,\n",
    "    output='m1d5_n2_1.out')  \n",
    "\n",
    "n2_singlet = moln2_singlet.RHF().run() #RHF Restricted Hartree-Fock. All electrons are paired\n",
    "n2_singlet.analyze()\n",
    "\n",
    "moln2_triplet = gto.M(\n",
    "    atom = 'N 0 0 0; N 0 0 1.1',  # in Angstrom\n",
    "    basis = '6-31g',\n",
    "    symmetry = True,\n",
    "    charge=0, \n",
    "    spin=2, #two unpaired electrons\n",
    "    verbose=3,\n",
    "    output='m1d5_n2_3.out')  \n",
    "\n",
    "n2_triplet = moln2_triplet.UHF().run() #UHF Unrestricted Hartree-Fock. Two unpaired electrons\n",
    "n2_triplet.analyze()\n",
    "\n",
    "print(n2_singlet.e_tot)\n",
    "print(n2_triplet.e_tot)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a884bb6c",
   "metadata": {},
   "source": [
    "* What spin state is most stable for N2 and O2? Compare it with their Lewis structures and assess if the Lewis structures can be used to predict the electronic stability.\n",
    "* Look at the orbital energies of the output to explain this phenomenon. In the case of O2, the molecules has 16 electrons. Check the \"Orbital Energies\" and fill them up with 16 electrons."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b6d38f7",
   "metadata": {},
   "source": [
    "## Population analysis\n",
    "\n",
    "Partial charges are not really a quantum property because the electrons are not constrained to one atom. Rather, the square of the wavefunction extents over all space. That being said, there are methods to decide on what atom the electrons are most probable to be found.\n",
    "\n",
    "Below we can calculate the charges according using the Mulliken method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d719c5b4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "overwrite output file: m1d5_co_sp.out\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "((array([1.99999186e+00, 1.73393973e+00, 1.17103529e-02, 4.46250338e-04,\n",
       "         4.13382584e-01, 4.13382584e-01, 1.01917143e+00, 2.29978258e-04,\n",
       "         2.29978258e-04, 1.11210209e-02, 2.19064334e-07, 2.19064334e-07,\n",
       "         4.58726884e-05, 9.82468610e-35, 4.68933508e-03, 7.18710213e-03,\n",
       "         4.68933508e-03, 2.25139348e-33, 1.99999328e+00, 1.79866648e+00,\n",
       "         4.66845337e-03, 9.99958974e-05, 1.57936995e+00, 1.57936995e+00,\n",
       "         1.40799362e+00, 4.82202347e-04, 4.82202347e-04, 1.44029039e-03,\n",
       "         4.04889972e-04, 4.04889972e-04, 9.93620981e-05, 3.78733935e-34,\n",
       "         1.44084349e-03, 3.42491238e-03, 1.44084349e-03, 9.88197247e-33]),\n",
       "  array([ 0.37978216, -0.37978216])),\n",
       " array([ 0.        ,  0.        , -1.67449812]))"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "molCO = gto.M(\n",
    "    atom = 'C 0 0 0; O 0 0 1.4',  # in Angstrom\n",
    "    basis = '6-31g',\n",
    "    symmetry = True,\n",
    "    charge=0, \n",
    "    spin=0,\n",
    "    verbose=4,\n",
    "    output='m1d5_co_sp.out')  \n",
    "\n",
    "co = molCO.RHF().run() #RHF Restricted Hartree-Fock. All electrons are paired\n",
    "co.analyze()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f494852",
   "metadata": {},
   "source": [
    "* Open the output file m1d5_co_sp.out. At the bottom of the file we added the dipole and Mulliken charges. What is the charge on C and O in the molecule of CO? Does it match with its Lewis structure?\n",
    "* Charges depend a lot on the size of the basis set and the level of theory. You could try to redo the above calculation with b3lyp/6-311G** and see if the charges are the same.\n",
    "\n",
    "Notice though that we do not know whether or not the geometry we introduced so far is the actual bond distance. At different geometries we may obtain different results."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b71bac46",
   "metadata": {},
   "source": [
    "# Scanning a distance to find bond distance\n",
    "\n",
    "We want to produce a binding energy curve for the $H_2$ molecule using different levels of theory. The binding energy is given by:\n",
    "\n",
    " $$E_{bind} = E(H_2) - 2E(H)  \\tag{1} $$\n",
    "\n",
    "Or if we build the binding curve, the binding energy would be the energy difference between the bottom of the well and the energy at which the H-H bond is infinitely apart."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "2f927696",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import constants\n",
    "\n",
    "hartree_value = 1  # Value in Hartree\n",
    "electronvolt_value = hartree_value * constants.physical_constants[\"Hartree energy in eV\"][0]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0975a63d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.3\n",
      "overwrite output file: m1d5_h2_scan0.3UHF.out\n",
      "0.4\n",
      "overwrite output file: m1d5_h2_scan0.4UHF.out\n",
      "0.5\n",
      "overwrite output file: m1d5_h2_scan0.5UHF.out\n",
      "0.6\n",
      "overwrite output file: m1d5_h2_scan0.6UHF.out\n",
      "0.7\n",
      "overwrite output file: m1d5_h2_scan0.7UHF.out\n",
      "0.8\n",
      "overwrite output file: m1d5_h2_scan0.8UHF.out\n",
      "0.9\n",
      "overwrite output file: m1d5_h2_scan0.9UHF.out\n",
      "1.0\n",
      "overwrite output file: m1d5_h2_scan1.0UHF.out\n",
      "1.1\n",
      "overwrite output file: m1d5_h2_scan1.1UHF.out\n",
      "1.2\n",
      "overwrite output file: m1d5_h2_scan1.2UHF.out\n",
      "1.3\n",
      "overwrite output file: m1d5_h2_scan1.3UHF.out\n",
      "1.4\n",
      "overwrite output file: m1d5_h2_scan1.4UHF.out\n",
      "1.5\n",
      "overwrite output file: m1d5_h2_scan1.5UHF.out\n",
      "1.6\n",
      "overwrite output file: m1d5_h2_scan1.6UHF.out\n",
      "1.7\n",
      "overwrite output file: m1d5_h2_scan1.7UHF.out\n",
      "1.8\n",
      "overwrite output file: m1d5_h2_scan1.8UHF.out\n",
      "1.9\n",
      "overwrite output file: m1d5_h2_scan1.9UHF.out\n",
      "2.0\n",
      "overwrite output file: m1d5_h2_scan2.0UHF.out\n"
     ]
    }
   ],
   "source": [
    "myGeom = 'H 0 0 0; H 0 0 mydist'  #using mydist as a string to replace in the loop by the actual distance\n",
    "myoutput = 'm1d5_h2_scanmydist.out'\n",
    "\n",
    "def runThisDist(mydist,method):\n",
    "    molh2 = gto.M(\n",
    "        atom = myGeom.replace('mydist',str(mydist)), #str: need to convert float to string\n",
    "        basis = '6-31g',\n",
    "        symmetry = True,\n",
    "        charge=0, \n",
    "        spin=0,\n",
    "        verbose=2,\n",
    "        output=myoutput.replace('mydist',str(mydist)+method))  \n",
    "    if method == \"UHF\": \n",
    "        h2 = molh2.UHF().run()\n",
    "    elif method == \"RHF\":\n",
    "        h2 = molh2.RHF().run()\n",
    "    else:\n",
    "        print(\"Wrong method\")\n",
    "    return(h2.e_tot*electronvolt_value)\n",
    "\n",
    "distances = [0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. , 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2. ]\n",
    "#the distance array can be done with np.linspace(0.3,2.0,18)\n",
    "listEne = []\n",
    "for r in distances:\n",
    "    print(r)\n",
    "    listEne.append( runThisDist(r,\"UHF\"))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "id": "be5a0b41",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj4AAAHFCAYAAADyj/PrAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAUoVJREFUeJzt3XlYVGX/BvB7GGDABQRBQdkE951wwyUhU8i9Ul81zQVNTTPLNu19BS2z0rLyLZdfhJqZ5laZueWSmmjqq5kbioaAQAoiICDr+f3xNAMjiwwMc2a5P9c1FzNnzgzfOY4zN895FoUkSRKIiIiILICV3AUQERERGQqDDxEREVkMBh8iIiKyGAw+REREZDEYfIiIiMhiMPgQERGRxWDwISIiIovB4ENEREQWg8GHiIiILAaDD1EtWrt2LRQKheZiZ2cHNzc3BAcHY8mSJbh9+3aZx0REREChUMhQrW4UCgUiIiJ0flxOTg4iIiJw+PDhMvepj1dcXFyN66Oqyc/Px/Tp0+Hu7g6lUonOnTtXuO/EiRNRr169Cu+vV68eJk6cqLl99epVvPbaawgICECDBg3g7OyMXr16YevWrXp8BUS6sZa7ACJLEBUVhdatW6OgoAC3b9/GsWPH8MEHH2DZsmXYvHkznnzySc2+U6ZMQWhoqIzVVk10dDQ8PDx0flxOTg4WLlwIAAgKCtK6b9CgQYiOjoa7u7s+SqQqWLlyJVavXo0VK1YgICCg0mCjq3379mHXrl0YP348unbtisLCQmzevBkjR47EwoULsWDBAr39LqKqYvAhMoD27dujS5cumtvPPvssXnnlFfTu3RvPPPMMrl27hsaNGwMAPDw8qhUoDK1Hjx56f05XV1e4urrq/XmNQUFBARQKBaytjetj98KFC7C3t8esWbP0/tyjR4/GzJkztVown3rqKaSmpuKDDz7Am2++CZVKpfffS1QZnuoikomXlxc++ugjZGVlYfXq1Zrt5Z3qOnjwIIKCgtCwYUPY29vDy8sLzz77LHJycjT7rFy5Ep06dUK9evVQv359tG7dGvPnz9d6ngsXLmDYsGFwcnKCnZ0dOnfujHXr1pWp7d69e5g7dy58fX2hUqnQqFEjDBw4EFeuXNHs8/Cprjt37uDFF19E27ZtUa9ePTRq1AhPPPEEjh49qtknLi5OE2wWLlyoOQWoPj1S0amur776Cp06dYKdnR2cnZ3x9NNP4/Lly1r7qE/DxMbGYuDAgahXrx48PT0xd+5c5OXlVfIvUWLjxo0IDAxEvXr1UK9ePXTu3BmRkZGa+318fLRO5agFBQVptV4dPnwYCoUCX3/9NebOnYumTZtCpVLh4sWLUCgUWs+ptnv3bigUCvz444+abdeuXcPYsWPRqFEjqFQqtGnTBp9//nmVXsuDBw8wb948NGvWDLa2tmjatClmzpyJe/fuafZRKBT48ssvkZubq/m3WLt2bZWevypcXFzKPW3brVs35OTk4O7du3r7XURVZVx/ehBZmIEDB0KpVOLIkSMV7hMXF4dBgwahT58++Oqrr9CgQQPcunULe/bsQX5+PurUqYNNmzbhxRdfxEsvvYRly5bBysoKsbGxuHTpkuZ5YmJi0LNnTzRq1AifffYZGjZsiA0bNmDixIn4+++/8cYbbwAAsrKy0Lt3b8TFxeHNN99E9+7dcf/+fRw5cgTJyclo3bp1uXWqv8TCw8Ph5uaG+/fvY8eOHQgKCsKBAwcQFBQEd3d37NmzB6GhoQgLC8OUKVMAoNJWniVLlmD+/PkYM2YMlixZgrS0NERERCAwMBCnTp1CixYtNPsWFBRg6NChCAsLw9y5c3HkyBG88847cHR0fORplQULFuCdd97BM888g7lz58LR0REXLlzAzZs3K31cZebNm4fAwECsWrUKVlZW8PT0hL+/P6KiohAWFqa179q1azUBEwAuXbqEnj17agKym5sb9u7di9mzZyM1NRXh4eEV/l5JkjB8+HAcOHAA8+bNQ58+fXD+/HmEh4cjOjoa0dHRUKlUiI6OxjvvvINDhw7h4MGDAAA/P79Hvq7CwsJqHxMAOHToEFxdXdGoUaMaPQ9RtUhEVGuioqIkANKpU6cq3Kdx48ZSmzZtNLfDw8Ol0v81t27dKgGQzp07V+FzzJo1S2rQoEGltYwePVpSqVRSfHy81vannnpKqlOnjnTv3j1JkiRp0aJFEgBp//79lT4fACk8PLzC+wsLC6WCggKpX79+0tNPP63ZfufOnQofqz5ef/31lyRJkpSeni7Z29tLAwcO1NovPj5eUqlU0tixYzXbJkyYIAGQvvvuO619Bw4cKLVq1arS13Ljxg1JqVRKzz33XKX7eXt7SxMmTCizvW/fvlLfvn01tw8dOiQBkB5//PEy+3722WcSACkmJkaz7e7du5JKpZLmzp2r2RYSEiJ5eHhIGRkZWo+fNWuWZGdnJ929e7fCOvfs2SMBkD788EOt7Zs3b5YASGvWrNFsmzBhglS3bt0Kn6s09TGu7FLe8Snt//7v/yQA0qefflql30mkbzzVRSQzSZIqvb9z586wtbXFCy+8gHXr1uHGjRtl9unWrRvu3buHMWPG4IcffkBqamqZfQ4ePIh+/frB09NTa/vEiRORk5OD6OhoAOKUS8uWLbU6XFfVqlWr8Nhjj8HOzg7W1tawsbHBgQMHypyWqqro6Gjk5uaWOb3k6emJJ554AgcOHNDarlAoMGTIEK1tHTt2fGSrzf79+1FUVISZM2dWq86KPPvss2W2Pffcc1CpVFqnlL799lvk5eVh0qRJAMRpqgMHDuDpp59GnTp1UFhYqLkMHDgQDx48wIkTJyr8verWm4eP28iRI1G3bt0yx00X9vb2OHXqVLkXe3v7Sh+7e/duzJw5EyNGjMBLL71U7RqIaoLBh0hG2dnZSEtLQ5MmTSrcx8/PD7/88gsaNWqEmTNnws/PD35+fvj00081+4wfPx5fffUVbt68iWeffRaNGjVC9+7dsX//fs0+aWlp5Y6WUv/utLQ0AKKvTnU6V3/88ceYMWMGunfvjm3btuHEiRM4deoUQkNDkZubq/Pzla6porrV96vVqVMHdnZ2WttUKhUePHhQ6e+5c+cOAOi9U3l5dTs7O2Po0KFYv349ioqKAIjTXN26dUO7du0AiNddWFiIFStWwMbGRuuiPhVWXrhVS0tLg7W1dZlTiAqFAm5ubmWOmy6srKzQpUuXci9WVhV/pezduxfPPPMM+vfvj2+++cYkpmwg88Q+PkQy2rVrF4qKisoM635Ynz590KdPHxQVFeH06dNYsWIF5syZg8aNG2P06NEAgEmTJmHSpEnIzs7GkSNHEB4ejsGDB+Pq1avw9vZGw4YNkZycXOa5k5KSAIiOqIDob5OYmKjza9mwYQOCgoKwcuVKre1ZWVk6P5daw4YNAaDCutU115Q6ICQmJpZpESvNzs6u3I7Sqamp5dZS0Zf7pEmTsGXLFuzfvx9eXl44deqU1nFzcnKCUqnE+PHjK2yFatasWYV1NmzYEIWFhbhz545W+JEkCSkpKejatWuFj60Ne/fuxfDhw9G3b19s27YNtra2Bv39RKWxxYdIJvHx8Xjttdfg6OiIadOmVekxSqUS3bt314zs+d///ldmn7p16+Kpp57C22+/jfz8fFy8eBEA0K9fPxw8eFATdNTWr1+POnXqaIanP/XUU7h69armdElVKRSKMkOTz58/rzmFpqbepyqtQIGBgbC3t8eGDRu0ticmJmpO3enDgAEDoFQqy4S2h/n4+OD8+fNa265evYqYmBidf1/Tpk0RFRWFqKgo2NnZYcyYMZr769Spg+DgYJw9exYdO3Yst3VFHQrLoz4uDx+3bdu2ITs7W2/HrSr27duH4cOHo3fv3vj+++85fJ1kxxYfIgO4cOGCpo/G7du3cfToUURFRUGpVGLHjh2VjmpatWoVDh48iEGDBsHLywsPHjzAV199BQCafjhTp06Fvb09evXqBXd3d6SkpGDJkiVwdHTU/HUfHh6On376CcHBwViwYAGcnZ3xzTffYNeuXfjwww/h6OgIAJgzZw42b96MYcOG4a233kK3bt2Qm5uLX3/9FYMHD0ZwcHC5dQ4ePBjvvPMOwsPD0bdvX8TExGDRokVo1qyZ1iig+vXrw9vbGz/88AP69esHZ2dnuLi4wMfHp8xzNmjQAP/5z38wf/58PP/88xgzZgzS0tKwcOFC2NnZVTqySRc+Pj6YP38+3nnnHeTm5mLMmDFwdHTEpUuXkJqaqplwcfz48Rg3bhxefPFFPPvss7h58yY+/PBDneceUiqVeP755/Hxxx/DwcEBzzzzjOb4q3366afo3bs3+vTpgxkzZsDHxwdZWVmIjY3Fzp07Kw2m/fv3R0hICN58801kZmaiV69emlFd/v7+GD9+vO4HqRqOHTuG4cOHw83NDfPnz8e5c+e07m/bti0cHBwMUguRhty9q4nMmXqUkvpia2srNWrUSOrbt6/03nvvSbdv3y7zmIdHdUVHR0tPP/205O3tLalUKqlhw4ZS3759pR9//FGzz7p166Tg4GCpcePGkq2trdSkSRNp1KhR0vnz57We+88//5SGDBkiOTo6Sra2tlKnTp2kqKioMjWkp6dLL7/8suTl5SXZ2NhIjRo1kgYNGiRduXJFsw8eGpmVl5cnvfbaa1LTpk0lOzs76bHHHpO+//57acKECZK3t7fW8//yyy+Sv7+/pFKptEYCPTyqS+3LL7+UOnbsKNna2kqOjo7SsGHDpIsXL2rtU9HopIePZ2XWr18vde3aVbKzs5Pq1asn+fv7ax2f4uJi6cMPP5R8fX0lOzs7qUuXLtLBgwcrHNW1ZcuWCn/X1atXNe+LikbQ/fXXX9LkyZOlpk2bSjY2NpKrq6vUs2dP6d13333ka8nNzZXefPNNydvbW7KxsZHc3d2lGTNmSOnp6Vr76Tqqq7J969atqzWqS33sK7ocOnSoSr+XSJ8UkvSIISVEREREZoJ9fIiIiMhiMPgQERGRxWDwISIiIovB4ENEREQWg8GHiIiILAaDDxEREVkMTmD4kOLiYiQlJaF+/fpcS4aIiMhESJKErKwsNGnSpNJ14xh8HpKUlFTpWj1ERERkvBISEipdcJjB5yH169cHIA4cp1InIiIyDZmZmfD09NR8j1eEwech6tNbDg4ODD5EREQm5lHdVNi5mYiIiCwGgw8RERFZDAYfIiIishgMPkRERGQxGHyIiIjIYjD4EBERkcVg8CEiIiKLweBDREREFoPBh4iIiCwGgw8RERFZDAYfIiIishgMPoaSmAgcOiR+EhERkSwYfAwhMhLw9gaeeEL8jIyUuyIiIiKLxOBT2xITgRdeAIqLxe3iYmDaNLb8EBERyYDBp7Zdu1YSetSKioDYWHnqISIismAMPrWtRQvA6qHDrFQCzZvLUw8REZEFY/CpbR4ewJo1gEIhbisUwOrVYjsREREZFIOPIYSFAR98IK4//ri4TURERAbH4GMoXbqIn0lJ8tZBRERkwRh8DMXXV/y8eVN0biYiIiKDY/AxFA8PwNoayM9nqw8REZFMTCb4LF68GD179kSdOnXQoEGDcvc5deoU+vXrhwYNGsDJyQkDBgzAuXPnDFpnhZRKwMdHXL9xQ9ZSiIiILJXJBJ/8/HyMHDkSM2bMKPf+rKwshISEwMvLCydPnsSxY8fg4OCAkJAQFBQUGLjaCjRrJn4y+BAREcnCWu4CqmrhwoUAgLVr15Z7f0xMDNLT07Fo0SJ4enoCAMLDw9GxY0fEx8fDz8/PUKVWTN3P56+/5K2DiIjIQplMi8+jtGrVCi4uLoiMjER+fj5yc3MRGRmJdu3awdvbu8LH5eXlITMzU+tSa9TBhy0+REREsjCb4FO/fn0cPnwYGzZsgL29PerVq4e9e/fi559/hrV1xQ1bS5YsgaOjo+aibi2qFQw+REREspI1+EREREChUFR6OX36dJWeKzc3F5MnT0avXr1w4sQJ/Pbbb2jXrh0GDhyI3NzcCh83b948ZGRkaC4JCQn6enllsY8PERGRrGTt4zNr1iyMHj260n181COhHmHjxo2Ii4tDdHQ0rP5ZG2vjxo1wcnLCDz/8UOHvUalUUKlUOtVdbeoWn7//BrKzgbp1DfN7iYiICIDMwcfFxQUuLi56ea6cnBxYWVlBoV4TC9DcLn54dXS5ODkBDRoA9+4BcXFAu3YyF0RERGRZTKaPT3x8PM6dO4f4+HgUFRXh3LlzOHfuHO7fvw8A6N+/P9LT0zFz5kxcvnwZFy9exKRJk2BtbY3g4GCZqy+F/XyIiIhkYzLBZ8GCBfD390d4eDju378Pf39/+Pv7a/oAtW7dGjt37sT58+cRGBiIPn36ICkpCXv27IG7u7vM1ZfCfj5ERESyMZl5fNauXVvhHD5q/fv3R//+/Q1TUHWxxYeIiEg2JtPiYzY4iSEREZFsGHwMjS0+REREsmHwMbTSwUeS5K2FiIjIwjD4GJqXF6BQALm5Yj4fIiIiMhgGH0OztQXUy2Kwnw8REZFBMfjIgf18iIiIZMHgIwcGHyIiIlkw+MiBkxgSERHJgsFHDmzxISIikgWDjxw4iSEREZEsGHzkoA4+iYlAXp68tRAREVkQBh85uLoCdeuKCQxv3pS7GiIiIovB4CMHhYIdnImIiGTA4CMX9vMhIiIyOAYfuXBkFxERkcEx+MiFwYeIiMjgGHzkwj4+REREBsfgI5fSLT6SJG8tREREFoLBRy4+PuJnZiaQni5rKURERJaCwUcudeoA7u7iOk93ERERGQSDj5zYz4eIiMigGHzkxLl8iIiIDIrBR04c0k5ERGRQDD5yYvAhIiIyKAYfOTH4EBERGRSDj5zUnZtv3gQKC+WthYiIyAIw+MipSRPA1hYoKgISE+WuhoiIyOwx+MjJyopD2omIiAyIwUdu7OdDRERkMAw+cmOLDxERkcEw+MiNkxgSEREZDIOP3Hiqi4iIyGAYfOTG4ENERGQwDD5yU/fxSU0FsrLkrYWIiMjMMfjIzcEBaNhQXGc/HyIiolrF4GMMeLqLiIjIIBh8jAGDDxERkUEw+BgDBh8iIiKDYPAxBuoOzuzjQ0REVKsYfIwBW3yIiIgMwiSCT1xcHMLCwtCsWTPY29vDz88P4eHhyM/P19ovPj4eQ4YMQd26deHi4oLZs2eX2ccolZ69ubhY3lqIiIjMmLXcBVTFlStXUFxcjNWrV6N58+a4cOECpk6diuzsbCxbtgwAUFRUhEGDBsHV1RXHjh1DWloaJkyYAEmSsGLFCplfwSN4egJKJZCXByQnA02byl0RERGRWVJIkiTJXUR1LF26FCtXrsSNf04P7d69G4MHD0ZCQgKaNGkCANi0aRMmTpyI27dvw8HBoUrPm5mZCUdHR2RkZFT5MXrh6ytafI4cAfr0MdzvJSIiMgNV/f42iVNd5cnIyICzs7PmdnR0NNq3b68JPQAQEhKCvLw8nDlzRo4SdcPFSomIiGqdSQaf69evY8WKFZg+fbpmW0pKCho3bqy1n5OTE2xtbZGSklLhc+Xl5SEzM1PrIgt2cCYiIqp1sgafiIgIKBSKSi+nT5/WekxSUhJCQ0MxcuRITJkyRes+hUJR5ndIklTudrUlS5bA0dFRc/H09NTPi9MVgw8REVGtk7Vz86xZszB69OhK9/Hx8dFcT0pKQnBwMAIDA7FmzRqt/dzc3HDy5Emtbenp6SgoKCjTElTavHnz8Oqrr2puZ2ZmyhN+GHyIiIhqnazBx8XFBS4uLlXa99atWwgODkZAQACioqJgZaXdWBUYGIjFixcjOTkZ7u7uAIB9+/ZBpVIhICCgwudVqVRQqVTVfxH6wkkMiYiIap1JjOpKSkpC37594eXlhfXr10OpVGruc3NzAyCGs3fu3BmNGzfG0qVLcffuXUycOBHDhw/XaTi7bKO60tIAdQjMyQHs7Q33u4mIiExcVb+/TWIen3379iE2NhaxsbHw8PDQuk+d25RKJXbt2oUXX3wRvXr1gr29PcaOHauZ58foOTsDDg5AZiYQFwe0aSN3RURERGbHJFp8DEm2Fh8A8PcHzp0DfvoJGDTIsL+biIjIhJn9PD5mif18iIiIahWDjzHhyC4iIqJaxeBjTBh8iIiIahWDjzFh8CEiIqpVDD7GRN3H58YNgH3OiYiI9I7Bx5h4ewMKBZCdDaSmyl0NERGR2WHwMSZ2dkDTpuI6T3cRERHpHYOPsWE/HyIiolrD4GNsGHyIiIhqDYOPseEkhkRERLWGwcfYsMWHiIio1jD4GBsGHyIiolrD4GNs1MEnIQEoKJC3FiIiIjPD4GNsGjcG7O2B4mIgPl7uaoiIiMwKg4+xUSi0Z3AmIiIivWHwMUbs50NERFQrGHyMEYMPERFRrWDwMUacy4eIiKhWMPgYI7b4EBER1QoGH2PE4ENERFQrGHyMkfpUV3q6uBAREZFeMPgYo7p1xXw+APv5EBER6RGDj7FiB2ciIiK9Y/AxVuznQ0REpHcMPsaKwYeIiEjvGHyMFYMPERGR3jH4GCv28SEiItI7Bh9jpW7xiYsDiopkLYWIiMhcMPgYq6ZNARsboKAAuHVL7mqIiIjMAoOPsVIqAR8fcZ39fIiIiPSCwceYqU93sZ8PERGRXjD4GDN1B2e2+BAREekFg48x45B2IiIivWLwMWYMPkRERHrF4GPMGHyIiIj0isHHmKn7+Ny+DWRny1sLERGRGWDwMWYNGgBOTuI6R3YRERHVGIOPsePpLiIiIr1h8DF2DD5ERER6w+Bj7LhYKRERkd4w+Bg7tvgQERHpjUkEn7i4OISFhaFZs2awt7eHn58fwsPDkZ+fr9nnjz/+wJgxY+Dp6Ql7e3u0adMGn376qYxV6wmDDxERkd5Yy11AVVy5cgXFxcVYvXo1mjdvjgsXLmDq1KnIzs7GsmXLAABnzpyBq6srNmzYAE9PTxw/fhwvvPAClEolZs2aJfMrqIHS63VJEqBQyFsPERGRCVNIkiTJXUR1LF26FCtXrsSNSlpCZs6cicuXL+PgwYNVft7MzEw4OjoiIyMDDg4O+ii1ZgoKADs7oLgYSE4G3NzkroiIiMjoVPX72yROdZUnIyMDzs7ONd4nLy8PmZmZWhejYmMDeHqK6zzdRUREVCMmGXyuX7+OFStWYPr06RXuEx0dje+++w7Tpk2r9LmWLFkCR0dHzcVTHTKMCfv5EBER6YWswSciIgIKhaLSy+nTp7Uek5SUhNDQUIwcORJTpkwp93kvXryIYcOGYcGCBejfv3+lNcybNw8ZGRmaS0JCgt5en94w+BAREemFrJ2bZ82ahdGjR1e6j4+Pj+Z6UlISgoODERgYiDVr1pS7/6VLl/DEE09g6tSp+Pe///3IGlQqFVQqlU51G1zpDs5ERERUbbIGHxcXF7i4uFRp31u3biE4OBgBAQGIioqClVXZxqqLFy/iiSeewIQJE7B48WJ9lysf9SSGbPEhIiKqEZMYzp6UlISgoCB4eXlh2bJluHPnjuY+t39GOV28eBHBwcEYMGAAXn31VaSkpAAAlEolXF1dZalbb3iqi4iISC9MIvjs27cPsbGxiI2NhYeHh9Z96tH4W7ZswZ07d/DNN9/gm2++0dzv7e2NuLg4Q5arf+rgc+sW8OCBGN5OREREOjPZeXxqi9HN4wOIiQsdHID794ErV4BWreSuiIiIyKiY/Tw+FkWhYAdnIiIiPWDwMRXs4ExERFRjDD6mgh2ciYiIaozBx1Qw+BAREdUYg4+pYB8fIiKiGmPwMRWl+/hwIB4REVG1MPiYCvXSHZmZwN27spZCRERkqhh8TIW9PdCkibjOfj5ERETVwuBjStjPh4iIqEYYfEwJ5/IhIiKqEQYfU8Ih7URERDXC4GNKGHyIiIhqhMHHlDD4EBER1QiDjylRB5/4eKCwUN5aiIiITBCDjylxcwNUKqCoCEhIkLsaIiIik8PgY0qsrDiyi4iIqAYYfEwN+/kQERFVG4OPqeEkhkRERNXG4GNqeKqLiIio2hh8TA1PdREREVUbg4+pYfAhIiKqNgYfU6M+1ZWWBmRmylsLERGRiWHwMTX16wMuLuI6OzgTERHphMHHFPF0FxERUbUw+JgiBh8iIqJqYfAxRZzLh4iIqFoYfEwRW3yIiIiqhcHHFHESQyIiomph8DFFpU91FRfLWwsREZEJYfAxRR4egLU1kJ8PJCXJXQ0REZHJYPAxRdbWgLe3uM4OzkRERFXG4GOq2M+HiIhIZ9bVfWBeXh5+//13xMXFIScnB66urvD390cz9Rcy1S6O7CIiItKZzsHn+PHjWLFiBb7//nvk5+ejQYMGsLe3x927d5GXlwdfX1+88MILmD59OurXr18bNRPA4ENERFQNOp3qGjZsGEaMGIGmTZti7969yMrKQlpaGhITE5GTk4Nr167h3//+Nw4cOICWLVti//79tVU3cRJDIiIinenU4jNgwABs2bIFtra25d7v6+sLX19fTJgwARcvXkQSRxzVHrb4EBER6UwhSZKkywMKCwthbV3trkFGLzMzE46OjsjIyICDg4Pc5VTs7l2gYUNxPTsbqFNH3nqIiIhkVNXvb51Hdbm7u+O1117DpUuXalQg1ZCTE+DoKK7HxclaChERkanQOfi8+uqr2LlzJzp06IDAwEBERkbi/v37tVEbVUahYD8fIiIiHekcfObNm4eYmBgcPnwYrVu3xpw5c+Du7o5Jkybht99+q40aqSLs50NERKSTak9g2KdPH0RFRSElJQWffPIJYmNj0adPH7Rq1QoffvihPmtEXFwcwsLC0KxZM9jb28PPzw/h4eHIz88vd/+0tDR4eHhAoVDg3r17eq3FqHASQyIiIp3UeObmunXrIiwsDEePHsXOnTuRmpqKefPm6aM2jStXrqC4uBirV6/GxYsXsXz5cqxatQrz588vd/+wsDB07NhRrzUYJbb4EBER6aTGw7NycnKwefNmREVF4bfffoOfnx9ef/11fdSmERoaitDQUM1tX19fxMTEYOXKlVi2bJnWvitXrsS9e/ewYMEC7N69W691GB328SEiItJJtYPP0aNHERUVha1bt6KoqAgjRozAu+++i8cff1yf9VUoIyMDzs7OWtsuXbqERYsW4eTJk7hRxVaQvLw85OXlaW5nZmbqtc5aVbrFR5JEh2ciIiKqkM6nut577z20bNkSQUFBuHjxIpYuXYrk5GSsW7fOYKHn+vXrWLFiBaZPn67ZlpeXhzFjxmDp0qXw8vKq8nMtWbIEjo6Omounp2dtlFw7vLxE2MnOBu7ckbsaIiIio6dz8Fm+fDkGDRqEP/74AydPnsS0adOqPdFfREQEFApFpZfTp09rPSYpKQmhoaEYOXIkpkyZotk+b948tGnTBuPGjdOphnnz5iEjI0NzSUhIqNZrkYVKBXh4iOvs50NERPRIOs/cXFBQABsbG61tDx48gJ2dnc6/PDU1FampqZXu4+Pjo3nupKQkBAcHo3v37li7di2srEpyW+fOnfHnn39C8c/pHkmSUFxcDKVSibfffhsLFy6sUk0mM3OzWlAQ8OuvwDffAGPHyl0NERGRLKr6/a1zHx916CkuLsbixYuxatUq/P3337h69Sp8fX3xn//8Bz4+PggLC3vkc7m4uMDFxaVKv/fWrVsIDg5GQEAAoqKitEIPAGzbtg25ubma26dOncLkyZNx9OhR+Pn56fAKTYyvrwg+7OBMRET0SNUezv7uu+9i7dq1+PDDD7UWLe3QoQO+/PJLvRSnlpSUhKCgIHh6emLZsmW4c+cOUlJSkJKSotnHz88P7du311ya/TPHTZs2bdCoUSO91mNUOKSdiIioyqo9qmv9+vVYs2YN+vXrp9XJuGPHjrhy5YpeilPbt28fYmNjERsbCw91n5Z/6HimzvxwEkMiIqIqq3aLz61bt9C8efMy24uLi1FQUFCjoh42ceJESJJU7qUiQUFBkCQJDRo00GstRoctPkRERFVW7eDTrl07HD16tMz2LVu2wN/fv0ZFkQ7UwScxEahgCQ8iIiISqn2qKzw8HOPHj8etW7dQXFyM7du3IyYmBuvXr8dPP/2kzxqpMo0aAXXqADk5QHw8UE4rHBEREQnVbvEZMmQINm/ejJ9//hkKhQILFizA5cuXsXPnTvTv31+fNVJlFAr28yEiIqqiGq3VFRISgpCQEH3VQtXl6wtcvMjgQ0RE9Ag1Xp29MhY/4spQuFgpERFRlegUfNq0aYONGzci/xGdaK9du4YZM2bggw8+qFFxVEUc2UVERFQlOp3q+vzzz/Hmm29i5syZGDBgALp06YImTZrAzs4O6enpuHTpEo4dO4ZLly5h1qxZePHFF2urbiqNwYeIiKhKdF6rCwCOHz+OzZs348iRI4iLi0Nubi5cXFzg7++PkJAQjBs3zmTnzzG5tboA0b+nfXugQQMgPV3uaoiIiAyuqt/f1Qo+5swkg09ODlC3rrh+9y7g5CRvPURERAZW1e9vnTs3Dxw4EBkZGZrbixcvxr179zS309LS0LZtW12flmqiTh3AzU1cZwdnIiKiCukcfPbu3Yu8vDzN7Q8++AB3797V3C4sLERMTIx+qqOqYz8fIiKiR9I5+Dx8ZoxnyowEJzEkIiJ6pFqdx4cMiC0+REREj6Rz8FEoFFAoFGW2kcw4iSEREdEj6bxkhSRJmDhxIlQqFQDgwYMHmD59Our+M6qodP8fMiC2+BARkbFLTASuXQNatAA8PGQpQefgM2HCBK3b48aNK7PP888/X/2KqHrUfXzi4oCiIkCplLUcIiIiLZGRwAsvAMXFgJUVsGYNEBZm8DI4j89DTHIeH0CEnTp1gPx8EX68veWuiIiISEhMFN9LxcUl25RK8X2lp5afWpvHh4yUUgn4+Ijr7OdDRETGIC0NWLECCA7WDj2A+IM9NtbgJTH4mBP28yEiIrkVFQG7dwOjRgFNmgCzZ5cfcJRKoHlzg5encx8fMmIMPkREJJerV4G1a4F164CkpJLtnTsDkyeLFp+5c0v6oa5eLUsHZwYfc8JJDImIyJCysoAtW4CvvgJ++61ku7MzMG4cMGmSCD5qzz4rWn+aNzedUV1kxDiXDxER1TZJAo4eBaKiROjJzhbbrayA0FARdoYMAf6Z9kaLh4dsgUeNwcec8FQXERHVloQEYP16EXiuXy/Z3rKlCDvjxwNNm8pXXxUx+JgT9amu27eB+/eBevXkrYeIiEzbgwfADz+IU1n794vWHkB8v/zrX6LvTmAgYEIrODD4mBNHR3Fe9e5dcbqrQwe5KyIiIlNQekblpk2B//1PhJ2NG4F790r269tXhJ1nnwX+WbHB1DD4mBtfXxF8btxg8CEiokcrPaOyQiGGoN+6VXK/pycwcSIwYQLg5ydbmfrC4GNufH2B06eBffuAgADZO5EREZERi48Hpk4tOYUlSSL02NqKVp1Jk4AnnjCrZZA4gaG5ycoSP7/4QkwPHhkpbz1ERGR8/v4bWLIE6NatJPSUtnWrOM3Vv79ZhR6ALT7mJTER2LOn5HZxMTBtGhASwpYfIiJLV1wMHDokJg7csQMoLCx/P6US8Pc3bG0GxOBjTq5dK5vc1WuhMPgQEVmm1FQxo/Lq1dpLRwQGij+Oc3KAl16SfUZlQ2HwMSctWogJpB5e/VaGtVCIiEhG6kkGV68Wp63y88X2+vXFfDvTpgEdO5bsP2SI7DMqGwqDjznx8ADWrCnpqKZQmH1yJyKiUtLTxSSDq1cDly+XbA8IAKZPB0aPLn+ONyOYUdlQGHzMTViY+DllingTT54sbz1ERFS7JAk4cUKEnc2bxaSDAFCnDjB2rGjd6dJF3hqNiEKSyuvObbkyMzPh6OiIjIwMODg4yF1O9WRnAy4u4s3/559A+/ZyV0RERPqWmQl88w2wahVw/nzJ9g4dROvOc8+JiW0tRFW/v9niY47q1gWefBL46Sfgxx8ZfIiIzMmZMyLsfPttyQKhdnZiCYlp04AePUxqCQlD4zw+5mrIEPFz50556yAioupJTBTDzxMTxfqLX34pTll16SKuZ2cDbdoAn3wCJCWJkVsmtm6WHNjiY64GDxY/T54EUlIANzd56yEioqp7eBkJW1sgL0/cp55Vefp0oE8fBh0dscXHXDVpAnTtKjq97doldzVERFRVf/0lRueqpyaRJBF6fHyADz8ULUAbNwKPP87QUw0MPuZs6FDx88cf5a2DiIge7fZt4N13S/5ofVhkJPD664Crq+FrMyMMPuZM3c9n/34gN1feWoiIqHynT4uVzz09gf/8B0hLK7uPUgm0bGn42syQSQSfuLg4hIWFoVmzZrC3t4efnx/Cw8ORr56JspS1a9eiY8eOsLOzg5ubG2bNmiVDxUaiY0fAy0uEngMH5K6GiIjUCgqATZuAnj1FC8/69WJ25e7dgQ0bxKgt9eKgFrCMhCGZROfmK1euoLi4GKtXr0bz5s1x4cIFTJ06FdnZ2Vi2bJlmv48//hgfffQRli5diu7du+PBgwe4ceOGjJXLTKEQp7v++19xukvd4ZmIiOTx999ihv2VK4HkZLHNxgYYNUqsl9W9e8m+gwZZzDIShmSyExguXboUK1eu1ASb9PR0NG3aFDt37kS/fv2q/bxmMYFhafv3AwMGiFFdt26JtbyIiMiwTp8GPvtMzKysPlvh5iZGZk2bxpG3elDV72+T/RbMyMiAs7Oz5vb+/ftRXFyMW7duoU2bNvDw8MCoUaOQkJAgY5VGoG9fsShdSoqY9IqIiAwjP19MMqg+nfX11yWns775Brh5EwgPZ+gxMJMMPtevX8eKFSswffp0zbYbN26guLgY7733Hj755BNs3boVd+/eRf/+/cvtC6SWl5eHzMxMrYtZsbUFQkPFdY7uIiKqfX//DSxaJIafjx0LREeL01njxom51U6cENttbeWu1CLJGnwiIiKgUCgqvZw+fVrrMUlJSQgNDcXIkSMxZcoUzfbi4mIUFBTgs88+Q0hICHr06IFvv/0W165dw6FDhyqsYcmSJXB0dNRcPD09a+31yobD2omIat+pU8Dzz4tBJeHhog+PmxuwcCEQHy9afLp1k7tKiydr5+ZZs2Zh9OjRle7j4+OjuZ6UlITg4GAEBgZizZo1Wvu5u7sDANq2bavZ5urqChcXF8THx1f4/PPmzcOrr76quZ2ZmWl+4WfgQDEq4Px50bTq7S13RUREpisxEbh2DWjRAmjUCNi2TfTfOXGiZJ8ePYDZs8UMy2zZMSqyBh8XFxe4uLhUad9bt24hODgYAQEBiIqKgtVDnXR79eoFAIiJiYHHP73f7969i9TUVHhX8kWvUqmgUqmq+QpMhLMz0KsXcOSIWLvLkof4ExHVxMNLSdSvL1ZJB8TprNGjxeisrl3lrZMqZBJ9fJKSkhAUFARPT08sW7YMd+7cQUpKClJSUjT7tGzZEsOGDcPLL7+M48eP48KFC5gwYQJat26N4OBgGas3EjzdRURUM4mJJaEHELMrZ2aKVp9Fi4CEBDEfD0OPUTOJ4LNv3z7Exsbi4MGD8PDwgLu7u+ZS2vr169G9e3cMGjQIffv2hY2NDfbs2QMbGxuZKjci6uBz+DCQkSFrKUREJqW4GPj5ZzEbvjr0lLZhg5hxuXFjw9dGOjPZeXxqi9nN41NamzbAlStiHolRo+SuhojIuD14IIadf/wxcOlS+fsolUBcHCcYNAJmP48PVYN67a6dO+Wtg4jImKWlAYsXi+HoU6aI0FO/PvDaa2J1dC4lYdJMYskK0pOhQ4GlS4Fdu4DCQsCa//xERBrXrwPLlwNRUUBOjtjm6Qm8/LIIQI6OYtuYMVxKwoTxm8+SBAYCDRuKv2Z++03M6kxEZOmio4GPPgK2bxcdlgHA31+08IwcKUZrlebhwcBjwniqy5IolSULlXJ0FxFZsqIiEXR69RJLSmzbJkLPwIHAgQNiiZ+xY8uGHjJ5DD6WRt3P58cfS/6yISKyFDk5wBdfAK1aickFjx8XEwxOngxcuCC6AjzxhJijh8wST3VZmgEDxH/y2FggJgZo3VruioiIat/ffwOffy5CT1qa2ObkBLz4opjUlQuFWgwGH0tTv774a2bPHtHqw+BDRObs8mUxHP3rr4G8PLHN1xd45RVg0iSgbl156yOD46kuS1T6dBcRkTlITAQOHRI/JUlM1jp4MNC2LfDllyL09OgBbN0KXL0qWnkYeiwSJzB8iFlPYKiWkCBWD7ayAlJSAFdXuSsiIqq+h9fP8vISCzID4vbw4WKEVs+espZJtYsTGFLFPD3FUE31NOxERKaqvPWzbt4EVCpgxgzRl3H7doYe0mDwsVRctJSITN3du2KNrPLWz9q8WXRkbtHC8HWRUWPwsVTqfj5794r1aIiITEVyMvD664C3N7B2bdn7lUogIMDgZZFpYPCxVI89BjRpAmRni06ARETG7sYNYPp0sYbWsmXA/ftAx47iVBfXz6IqYvCxVAoFT3cRkWn480/guefEaavVq4H8fDHj8k8/AefOiW1xcWJUV1wcEBYmc8FkzBh8LFnp4MPBfURkbE6cAIYNE606GzeKvjyhocCRI8CxY8CgQSUzLHt4AEFBbOmhR2LwsWTBwWIei1u3gLNn5a6GiEj8EbZ/v5hoNTBQ/GGmUIjFQs+cAXbvBvr0kbtKMmEMPpbMzk4sYQEAO3fKWwsRWbbiYjHsvGtX8bl06BBgbS3W0Lp8GfjuO9E3kaiGGHwsHfv5EJGcCgqAdeuA9u3FoqFnzgD29sDs2cD162Jywlat5K6SzAjX6rJ06nPk//ufmAiM58eJyBByc0WoWboUiI8X2xwdxVISL7/MGeWp1rDFx9K5uorz6ABPdxFR7cvIAJYsEUPSX3pJhJ7GjYEPPhDX332XoYdqFYMPlZzuYvAhIn0qvXDo7dvA/PliHa3588VtHx/g88+Bv/4C3ngDMNf1Ecmo8FQXieDz1lvAgQNiQrB69eSuiIhM3cMLhyqVQGGhuK9tW/GZM3o0YGMjb51kcdjiQ0Dr1kDz5mJSsH375K6GiExdeQuHFhYCnToBO3aICQnHj2foIVkw+JD4a0y9dhdHdxFRTcTGAlOmlL9w6PLlwPDhgBW/ekg+fPeRoO7ns2sXUFQkby1EZHquXQMmTBAtyHv3lr1fqeRK6WQUGHxI6NULcHICUlPFNPFERFUREyNOW7VuDaxfL/5weuop0YGZC4eSEWLwIcHGBhg4UFzn6S4iepTLl8XCoW3bAhs2iFNbgwYBJ08CP/8MLF7MhUPJKDH4UAn28yGiR7l0CRgzBmjXrmTh0KFDgdOnxWrp3bqV7MuFQ8kIMfhQidBQsTbOlSvifD0RkdqFC8C//iWWlti0SYzUGj5czPr+ww9AQIDcFRJVCYMPlXB0FH+dAZzMkIiE8+eBESOADh3EQqGSBDzzDHD2rBia7u8vd4VEOmHwIW1ctJSIAODcORFwOnUCtm0T20aMAP74Q9zu3FnO6oiqjcGHtKn7+Rw7Bty9K28tRGR4//ufOIXl7y9adBQKYNQoMengli1Ax45yV0hUIww+pM3HRzRpFxUBu3fLXQ0RGcqZM6LFNyBA9NlRKMSSEhcuAJs3i749RGaAwYfK4ukuIvNWevHQU6eAwYOBLl1E3z4rK2DsWODiReDbb8VwdSIzwkVKqayhQ8UcHLt3i/W7bG3lroiI9KX04qGlqQPPv/8NtGolT21EBsAWHyqrSxfAzQ3IygKOHJG7GiLSl4cXD1UbMUJMSPj11ww9ZPYYfKgsKyvR9A3wdBeRufjjDzEPT3mLh86cCbRsafiaiGTA4EPlK93PR5LkrYWIqu/SJTEqq3Nn4PjxsvcrlUDz5gYvi0guDD5Uvn79ADs74OZNMYyViExLbKxYPLR9ezEMHRCjtN55h4uHkkVj52YqX506QP/+YpTHzp2cu4PIVNy8KcLN2rViWgoAePppYOFCMVUFAEycKIJR8+YMPWRx2OJDFeOwdiLTceuW6KvTooUYuVVUBAwcKBYP3b69JPQAXDyULJpJBJ+4uDiEhYWhWbNmsLe3h5+fH8LDw5Gfn6+136lTp9CvXz80aNAATk5OGDBgAM6dOydP0eZA3cH599+B5GR5ayGi8t2+Dbz6KuDnB3zxBVBQIE5VHz8O7NrFxUOJHmISwefKlSsoLi7G6tWrcfHiRSxfvhyrVq3C/PnzNftkZWUhJCQEXl5eOHnyJI4dOwYHBweEhISgoKBAxupNmJsb0K2buP7TT/LWQkTa0tKAefOAZs2A5cuBvDygd28xMeEvvwCBgXJXSGSUFJJkmkN2li5dipUrV+LGjRsAgNOnT6Nr166Ij4+Hp6cnAODPP/9Ex44dERsbCz8/vyo9b2ZmJhwdHZGRkQEHB4daq99kLF4sJjQbMoSnvIiMQUaGCDoffyzm2gKArl1Fv54BA8RSE0QWqKrf3ybR4lOejIwMODs7a263atUKLi4uiIyMRH5+PnJzcxEZGYl27drB29u7wufJy8tDZmam1oVKUffz2b8fyMmRtxYiS3b/PrBkiWjhWbhQhJ5OncQfJCdPAiEhDD1EVWCSwef69etYsWIFpk+frtlWv359HD58GBs2bIC9vT3q1auHvXv34ueff4a1dcWD15YsWQJHR0fNRd1aRP9o314sXPrggWg+JyLDys0VrTu+vsD8+UB6OtCmjRii/r//idZYBh6iKpM1+EREREChUFR6OX36tNZjkpKSEBoaipEjR2LKlCma7bm5uZg8eTJ69eqFEydO4LfffkO7du0wcOBA5ObmVljDvHnzkJGRobkkJCTU2us1SQqF+GAFeKqLyJDy8oD//ld0Wp47F7hzRww/37BBzK01YoSYZZ2IdCJrH5/U1FSkpqZWuo+Pjw/s7OwAiNATHByM7t27Y+3atbAq9Z8+MjIS8+fPR3JysmZ7fn4+nJycEBkZidGjR1epJvbxKccvv4g5fRo3BpKS+GFLVBsSE4Fr10QL64EDwKJFgPoPMW9vYMEC4PnngUpasIksWVW/v2X9H+Ti4gIXF5cq7Xvr1i0EBwcjICAAUVFRWqEHAHJycmBlZQVFqSZf9e3i8tamoap7/HHAwQH4+2/g1Cmge3e5KyIyLxWtmN6kiRhcEBYG2NrKUxuRmTGJP92TkpIQFBQET09PLFu2DHfu3EFKSgpSUlI0+/Tv3x/p6emYOXMmLl++jIsXL2LSpEmwtrZGcHCwjNWbAVtb4KmnxHWe7iLSr4QEYOrUsqEnPFzMrjxjBkMPkR6ZRPDZt28fYmNjcfDgQXh4eMDd3V1zUWvdujV27tyJ8+fPIzAwEH369EFSUhL27NmjtR9VE/v5EOmXJInRkk8+Wf5CwEFBgL29wcsiMncmO49PbWEfnwrcvQs0aiSmwb9xQwypJaLqiY4WI7QOHy7/fqUSiIvjkhJEOjD7eXzIwJydgT59xPWdO+WthchUnT8v5sbq2VOEHltbYM4cMVydK6YTGQSDD1UdFy0lqp7YWOC554DOncUfDkolMGWKGMW1fDnwyiuihefQIfEzLEzmgonMF091PYSnuioRGytWfra2BlJTAUdHuSsiMm63bolh6erV0gFg1CixrVUreWsjMjM81UX617y5mDG2sBDYs0fuaoiMV2oq8Npr4v/MmjUi9AwcKGZa3ryZoYdIRgw+pBue7iKqWGamWEfL1xf46COx1EufPsDRo8CuXYC/v9wVElk8Bh/SjTr4/PwzUFAgby1ExqL0eloREWIBUX9/YPdu4Ndfgd695a6QiP7B4EO66d4dcHEB7t0Djh2TuxoieRUUAP/3f6Lv29y5QFqaOI313XfA6dNAaCgXECUyMgw+pBulEhg8WFznsHayVMXFwLffAm3biqUmbt0CPD1FJ+YLF4CRI7mmHZGR4v9M0l3pfj4cFEiWRJKAn34Sp7HGjhUjHV1dgU8+EUPTJ0/mIqJERo7/Q0l3/fsDKhVw/Tpw+bL4q5fIXKlXTU9NFQHn+HGx3cEBeOMN4OWXgXr1ZC2RiKqOLT6ku3r1gCeeENeXLxdfDETmKDIS8PYW7/dRo0TosbcH3nwT+Osv4O23GXqITAyDD1VPw4bi55dfii+GyEh56yHSt0OHxOzKpVdNVyjE0PT33xfLuBCRyWHwId0lJgIbN5bcLi4Gpk1jyw+Zh7g4YNIkoF+/svdJkhiqTkQmi8GHdHftmvZfwYCYmTY2Vp56iPQhJQWYPRto2RJYu7b8jvtKpZiNmYhMFoMP6a5Fi7JDda2s+IVApik9HZg/H/DzA1asEHPzPPkkcPKkOJXLVdOJzApHdZHuPDzE+kPTppUsvGhtDWRny1sXkS6ys4HPPgM+/FBMyAmICTrfe6+k8363bkBIiGjNbN6coYfIDLDFh6onLEz0hfjlF6BnTyA/Hxg9GsjLk7syosrl5QH//a9o4Zk/X4Se9u2BH34AoqNLQo+ahwcQFMTQQ2QmGHyo+jw8RAfQLVvEMhbnzol5TYiMUVGR6LvTqhXw0kvA33+LtbU2bBDv3aFDubwEkQVg8KGaa9IEWLdOXP/sMy5lQcZFkoBt24AOHcRorZs3AXd3YOVK4MoV4LnnSvrxEJHZY/Ah/Rg4EHjlFXF90iSxdhGRnCQJ2LcP6NoVGDFCzDLu7Cz69MTGAtOnAzY2cldJRAbG4EP6s2QJ8NhjYoXqceNKOj4TGdrx40BwsOiYfOaMmF15wQLgxg3g9deBOnXkrpCIZMLgQ/qjUokVq+vWBQ4fFkGIyJD++AMYMgTo1Qv49VfxnnzlFRF4Fi4EHB3lrpCIZMbgQ/rVsiXwxRfiekQE8NtvspZDFuLaNbFaeufOYvV0pVIsN3HtGvDxx2IFdSIiMPhQbXj++ZJTXWPHigniiPQpMVGspXXqlJhPqk0b0doIAP/6F3DpEvB//wd4espbJxEZHU5gSLXjiy+AEydEJ9IpU4CtWzlUmPQjMhJ44YWyy6YMGgS8+65o9SEiqgBbfKh21K8PbNokRs1s3y6m+ieqqfPngalTy4aebdvEKS6GHiJ6BAYfqj0BAcD774vrr7wCXLggbz1kutLSxCzL3buXv3ios7PhayIik8TgQ7VrzhzgqaeABw9E34ucHLkrIlNy9y7w738DzZqJUYIPHpTdhyumE5EOGHyodllZiWUC3NxEh1P1JIdElUlPF/PuNGsGLF4MZGWJ01jffy86LXPFdCKqJnZuptrXqJFYD6l/f7Gq+5NPAiNHyl0VGaN794BPPgGWLwcyM8W2jh3F1AjDh5d0kA8N5YrpRFQtbPEhw+jXD3jrLXF96lSxsjuRWkYGsGgR4OMjJhrMzBQrpm/dCpw9Czz9tPaoQK6YTkTVxOBDhrNwIdCjh/iSGzsWKCyUuyKSW2amGILu4wOEh4v3Rrt2wHffiVmYn31WnC4lItITfqKQ4djYABs3Ag4OQHS0OH1BlikrS/Td8fEB/vMfcYqrbVtg82YxZH3kSAYeIqoV/GQhw2rWTHROBYD33gMOHpS3HjKsrCwxOsvHR4zWSk8HWrcWsy6fPw+MGsXAQ0S1ip8wZHijRonZnCVJLG1x547cFVFtu38f+OADEXznzxfD1Fu1Ar75RszvNHp0yUgtIqJaxOBD8vj0U7G+UnIyMGlS+ZPSkenLzgaWLhWB5623xESELVoAX38NXLwo+nox8BCRATH4kDzq1BH9OVQqYNcuEYTI9KkXD712DfjoI8DXF3jjDSA1VQw9X7dOzOc0bhwDDxHJgvP4kHw6dAA+/hiYOVN8OT7+OPDYY3JXRdVV0eKhvr6iA/O4cYA1P3KISF5s8SF5zZgh5mgpKBD9PLKy5K6IquOPP8pfPHTpUuDKFWDiRIYeIjIKDD4kL4UC+PJLwNNTnB6ZNUvuikgXly8D06YBXbuW30+rSxcxjQERkZEwmeAzdOhQeHl5wc7ODu7u7hg/fjySkpK09omPj8eQIUNQt25duLi4YPbs2cjPz5epYqoyZ2cxusfKCli/XixvQcZLkoBffgEGDhRz76xZI1rsHsbFQ4nICJlM8AkODsZ3332HmJgYbNu2DdevX8eIESM09xcVFWHQoEHIzs7GsWPHsGnTJmzbtg1z586VsWqqsj59xMy9gDj9de2avPVQWXl5YsHZTp3Eumu7d4sWu6efBo4c4eKhRGQSFJJkmuOIf/zxRwwfPhx5eXmwsbHB7t27MXjwYCQkJKBJkyYAgE2bNmHixIm4ffs2HBwcqvS8mZmZcHR0REZGRpUfQ3pSVAQ88YT4Eg0IAI4fB2xt5a6K7twBVq0CPv8c+Ptvsa1uXWDyZODllwE/v5J9ExO5eCgRyaKq398m2dvw7t27+Oabb9CzZ0/Y/NN/IDo6Gu3bt9eEHgAICQlBXl4ezpw5g+Dg4HKfKy8vD3l5eZrbmeoVocnwlEpxyqtTJ+DMGTHR3bJlcldluS5fFiulr18PPHggtnl4ALNniwkonZzKPsbDg4GHiIyayZzqAoA333wTdevWRcOGDREfH48ffvhBc19KSgoaN26stb+TkxNsbW2RkpJS4XMuWbIEjo6Omounp2et1U9V4OEBfPWVuP7RR+J0ChmOuv/OoEEl/XcePBCdlDduBG7cAF5/vfzQQ0RkAmQNPhEREVAoFJVeTp8+rdn/9ddfx9mzZ7Fv3z4olUo8//zzKH2mTqFQlPkdkiSVu11t3rx5yMjI0FwSEhL0+yJJd8OGlYzumjBBzO5Mtevh/js//6zdf+f334ExYzhCi4hMnqynumbNmoXRo0dXuo+Pj4/muouLC1xcXNCyZUu0adMGnp6eOHHiBAIDA+Hm5oaTJ09qPTY9PR0FBQVlWoJKU6lUUKlUNXodVAuWLgWOHhXzwzz/PLB3LxevrA269N8hIjIDsgYfdZCpDnVLj7p/TmBgIBYvXozk5GS4u7sDAPbt2weVSoWAgAD9FEyGY2cHbNokOjn/8gvw9tvAgAFinSf2Iam56vTfISIyAyYxquv333/H77//jt69e8PJyQk3btzAggULkJycjIsXL0KlUqGoqAidO3dG48aNsXTpUty9excTJ07E8OHDsWLFiir/Lo7qMjJffQWEhZXctrIS/U5Kb6PKJSaK6QGaNweuXhXLhPz8c8n9AQHA3LnAiBE8lUVEJsusRnXZ29tj+/btCA8PR3Z2Ntzd3REaGopNmzZpTlMplUrs2rULL774Inr16gV7e3uMHTsWyzgqyLT17699u7hYrAf15JOAt7c8NZmSitbPUihEX6pXXwV69xa3iYgsgEm0+BgSW3yMzKFDYm6fh7m7i1aKyZN5WqY89+6J01gvv1z2vkmTxFQBnFWZiMxIVb+/2VuUjFuLFuV3ak5OBl57DWjaVLRonD9v+NqMTXKy6KgcEgK4upYfegDRWZyhh4gsFIMPGTcPD9Gnp/RSCF98IbZ16ADk5oqlEjp1Avr2BbZuLX/dKHMVGytGwPXsKULgjBnAvn1AYSHQsmXZU1hcP4uILBxPdT2Ep7qMVHlLIUiSGPK+YgWwY4dY8gIoCQBTpwKNGslXc22QJODcOfF6d+wALlzQvr9HDzH3ztNPi9ayyEixenpRUcn6WewYTkRmqKrf3ww+D2HwMVGJieI0z5o1Ym4aQKzz9a9/ickQu3WTt76aKCoCfvtNBJ3vvwfi4krus7YGgoNF0Bk2DCi1ZIsG188iIgvA4FNNDD4mLi8P+O474L//FbMNq3XrJgLQqFGAKUxYmZcHHDgAbN8O/PhjSZgDAHt7IDQUeOYZsbQEO3cTETH4VBeDjxn5/XcRgDZvBvLzxTZXV9EZevp0+Vs/1PPrqCdlzMwUa5Pt2CHm2cnKKtnXyQkYMkS07AwYANSpI1/dRERGiMGnmhh8zNDt26ID9MqVwK1bYptSKULESy8BffoYfh6b0vPrKBRA+/ZATExJQANEX6Xhw0Wdjz/OyQWJiCrB4FNNDD5mrKAA+OEH0Rn6yJGS7R07itNgY8cC6enarTC6yMsT8+eoL+np5d9OSgJ++qn852jZUgSdZ54RK6JzfTIioiph8KkmBh8Lcf68OA22YYMYEg+IvjMPHoiRU1ZWwJw5QGDgo4OM+qJ+nuqKihKr0XMWZSIinTH4VBODj4VJTxfrgX36KZCQUPPnUygAR0egQYOSi5OT9m0AiIgQAUtNqRSjteTud0REZKLMaq0uolrj5CSWvujUqey6YIDoe9OsWeVBpvRtB4eqnZ5q2rTs/DoMPUREtY7BhwgAWrcWgaX0Yp5KpRhlVRuBJCxMLC3B+XWIiAyKPSeJgPKXxqjtVhgPDyAoiKGHiMiA2OJDpMZWGCIis8fgQ1SahwcDDxGRGeOpLiIiIrIYDD5ERERkMRh8iIiIyGIw+BAREZHFYPAhIiIii8HgQ0RERBaDwYeIiIgsBoMPERERWQwGHyIiIrIYDD5ERERkMRh8iIiIyGJwra6HSJIEAMjMzJS5EiIiIqoq9fe2+nu8Igw+D8nKygIAeHp6ylwJERER6SorKwuOjo4V3q+QHhWNLExxcTGSkpJQv359KBSKcvfJzMyEp6cnEhIS4ODgYOAKjQePQwkeC4HHQeBxEHgcSvBYCLV5HCRJQlZWFpo0aQIrq4p78rDF5yFWVlbw8PCo0r4ODg4W/QZW43EowWMh8DgIPA4Cj0MJHguhto5DZS09auzcTERERBaDwYeIiIgsBoNPNahUKoSHh0OlUsldiqx4HErwWAg8DgKPg8DjUILHQjCG48DOzURERGQx2OJDREREFoPBh4iIiCwGgw8RERFZDAYfIiIishgMPgC++OILNGvWDHZ2dggICMDRo0cr3Hf79u3o378/XF1d4eDggMDAQOzdu1drn7Vr10KhUJS5PHjwoLZfSo3pciwOHz5c7uu8cuWK1n7btm1D27ZtoVKp0LZtW+zYsaO2X0aN6XIcJk6cWO5xaNeunWYfU3xPHDlyBEOGDEGTJk2gUCjw/fffP/Ixv/76KwICAmBnZwdfX1+sWrWqzD6m9n7Q9TiY82eErsfCXD8jdD0O5voZsWTJEnTt2hX169dHo0aNMHz4cMTExDzycXJ/Tlh88Nm8eTPmzJmDt99+G2fPnkWfPn3w1FNPIT4+vtz9jxw5gv79++Pnn3/GmTNnEBwcjCFDhuDs2bNa+zk4OCA5OVnrYmdnZ4iXVG26Hgu1mJgYrdfZokULzX3R0dH417/+hfHjx+OPP/7A+PHjMWrUKJw8ebK2X0616XocPv30U63Xn5CQAGdnZ4wcOVJrP1N7T2RnZ6NTp07473//W6X9//rrLwwcOBB9+vTB2bNnMX/+fMyePRvbtm3T7GOK7wddj4M5f0boeizUzO0zQtfjYK6fEb/++itmzpyJEydOYP/+/SgsLMSAAQOQnZ1d4WOM4nNCsnDdunWTpk+frrWtdevW0ltvvVXl52jbtq20cOFCze2oqCjJ0dFRXyUajK7H4tChQxIAKT09vcLnHDVqlBQaGqq1LSQkRBo9enSN660tNX1P7NixQ1IoFFJcXJxmm6m+J9QASDt27Kh0nzfeeENq3bq11rZp06ZJPXr00Nw2xfdDaVU5DuUxl8+I0qpyLMz1M6K06rwnzPEzQpIk6fbt2xIA6ddff61wH2P4nLDoFp/8/HycOXMGAwYM0No+YMAAHD9+vErPUVxcjKysLDg7O2ttv3//Pry9veHh4YHBgweX+WvP2NTkWPj7+8Pd3R39+vXDoUOHtO6Ljo4u85whISFVPr6Gpo/3RGRkJJ588kl4e3trbTe194SuKvq3Pn36NAoKCirdx1jfD/pgLp8RNWFOnxH6YK6fERkZGQBQ5r1emjF8Tlh08ElNTUVRUREaN26stb1x48ZISUmp0nN89NFHyM7OxqhRozTbWrdujbVr1+LHH3/Et99+Czs7O/Tq1QvXrl3Ta/36VJ1j4e7ujjVr1mDbtm3Yvn07WrVqhX79+uHIkSOafVJSUmp0fA2tpu+J5ORk7N69G1OmTNHaborvCV1V9G9dWFiI1NTUSvcx1veDPpjLZ0R1mONnRE2Z62eEJEl49dVX0bt3b7Rv377C/Yzhc4KrswNQKBRatyVJKrOtPN9++y0iIiLwww8/oFGjRprtPXr0QI8ePTS3e/XqhcceewwrVqzAZ599pr/Ca4Eux6JVq1Zo1aqV5nZgYCASEhKwbNkyPP7449V6TmNR3ZrXrl2LBg0aYPjw4VrbTfk9oYvyjtvD203x/VBd5vgZoQtz/oyoLnP9jJg1axbOnz+PY8eOPXJfuT8nLLrFx8XFBUqlskyKvH37dpm0+bDNmzcjLCwM3333HZ588slK97WyskLXrl2NOrnX5FiU1qNHD63X6ebmVuPnNKSaHAdJkvDVV19h/PjxsLW1rXRfU3hP6Kqif2tra2s0bNiw0n2M9f1QE+b2GaEvpv4ZURPm+hnx0ksv4ccff8ShQ4fg4eFR6b7G8Dlh0cHH1tYWAQEB2L9/v9b2/fv3o2fPnhU+7ttvv8XEiROxceNGDBo06JG/R5IknDt3Du7u7jWuubZU91g87OzZs1qvMzAwsMxz7tu3T6fnNKSaHIdff/0VsbGxCAsLe+TvMYX3hK4q+rfu0qULbGxsKt3HWN8P1WWOnxH6YuqfETVhbp8RkiRh1qxZ2L59Ow4ePIhmzZo98jFG8Tmhly7SJmzTpk2SjY2NFBkZKV26dEmaM2eOVLduXU1v+7feeksaP368Zv+NGzdK1tbW0ueffy4lJydrLvfu3dPsExERIe3Zs0e6fv26dPbsWWnSpEmStbW1dPLkSYO/Pl3oeiyWL18u7dixQ7p69ap04cIF6a233pIASNu2bdPs89tvv0lKpVJ6//33pcuXL0vvv/++ZG1tLZ04ccLgr6+qdD0OauPGjZO6d+9e7nOa4nsiKytLOnv2rHT27FkJgPTxxx9LZ8+elW7evClJUtnjcOPGDalOnTrSK6+8Il26dEmKjIyUbGxspK1bt2r2McX3g67HwZw/I3Q9Fub6GaHrcVAzt8+IGTNmSI6OjtLhw4e13us5OTmafYzxc8Lig48kSdLnn38ueXt7S7a2ttJjjz2mNRRvwoQJUt++fTW3+/btKwEoc5kwYYJmnzlz5kheXl6Sra2t5OrqKg0YMEA6fvy4AV9R9elyLD744APJz89PsrOzk5ycnKTevXtLu3btKvOcW7ZskVq1aiXZ2NhIrVu31vrQM1a6HAdJkqR79+5J9vb20po1a8p9PlN8T6iHIlf0Xi/vOBw+fFjy9/eXbG1tJR8fH2nlypVlntfU3g+6Hgdz/ozQ9ViY62dEdf5vmONnRHnHAIAUFRWl2ccYPycU/xRPREREZPYsuo8PERERWRYGHyIiIrIYDD5ERERkMRh8iIiIyGIw+BAREZHFYPAhIiIii8HgQ0RERBaDwYeIzFJQUBDmzJkDAPDx8cEnn3wiaz1EZBwYfIjI7J06dQovvPBClfZlSCIyb9ZyF0BEVNtcXV3lLoGIjARbfIjI5GVnZ+P5559HvXr14O7ujo8++kjr/odbcSIiIuDl5QWVSoUmTZpg9uzZAMTpsZs3b+KVV16BQqGAQqEw5MsgIgNg8CEik/f666/j0KFD2LFjB/bt24fDhw/jzJkz5e67detWLF++HKtXr8a1a9fw/fffo0OHDgCA7du3w8PDA4sWLUJycjKSk5MN+TKIyAB4qouITNr9+/cRGRmJ9evXo3///gCAdevWwcPDo9z94+Pj4ebmhieffBI2Njbw8vJCt27dAADOzs5QKpWoX78+3NzcDPYaiMhw2OJDRCbt+vXryM/PR2BgoGabs7MzWrVqVe7+I0eORG5uLnx9fTF16lTs2LEDhYWFhiqXiGTG4ENEJk2SJJ329/T0RExMDD7//HPY29vjxRdfxOOPP46CgoJaqpCIjAmDDxGZtObNm8PGxgYnTpzQbEtPT8fVq1crfIy9vT2GDh2Kzz77DIcPH0Z0dDT+/PNPAICtrS2KiopqvW4ikgf7+BCRSatXrx7CwsLw+uuvo2HDhmjcuDHefvttWFmV/3fd2rVrUVRUhO7du6NOnTr4+uuvYW9vD29vbwBiBNiRI0cwevRoqFQquLi4GPLlEFEtY/AhIpO3dOlS3L9/H0OHDkX9+vUxd+5cZGRklLtvgwYN8P777+PVV19FUVEROnTogJ07d6Jhw4YAgEWLFmHatGnw8/NDXl6ezqfSiMi4KST+ryYiIiILwT4+REREZDEYfIiIiMhiMPgQERGRxWDwISIiIovB4ENEREQWg8GHiIiILAaDDxEREVkMBh8iIiKyGAw+REREZDEYfIiIiMhiMPgQERGRxWDwISIiIovx/xgmughzgVu+AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "#the x axis is the atomic number\n",
    "\n",
    "ax.plot(distances,listEne,'.r-')\n",
    "ax.set_xlabel(\"dist\")\n",
    "ax.set_ylabel(\"E(eV)\")\n",
    "ax.set_title(\"Dissociation curve of H2\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c68214f",
   "metadata": {},
   "source": [
    "* Your turn: Build the bonding curve for H2(+) and H2(-).\n",
    "* Plot these curves."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "598d2782",
   "metadata": {},
   "source": [
    "## Homolytic and heterolytic breaking of a bond\n",
    "\n",
    "One can break a bond homolitically, that is, each atom keeps one electron so that at the end you have two unpaired electrons (triplet) and both atoms remain neutral. Or you can break it heterolitically, so that one atom keeps both electrons and no electrons are unpaired. In gas phase (without solvent) the homolytical bond breaking takes less energy.\n",
    "\n",
    "The plot shown above was an homolytical breaking because we allowed the electrons to be separated. If we run the same H2 curve just using RHF we will obtain the heterolytical breaking."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "id": "27914a81",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.3\n",
      "overwrite output file: m1d5_h2_scan0.3RHF.out\n",
      "0.4\n",
      "overwrite output file: m1d5_h2_scan0.4RHF.out\n",
      "0.5\n",
      "overwrite output file: m1d5_h2_scan0.5RHF.out\n",
      "0.6\n",
      "overwrite output file: m1d5_h2_scan0.6RHF.out\n",
      "0.7\n",
      "overwrite output file: m1d5_h2_scan0.7RHF.out\n",
      "0.8\n",
      "overwrite output file: m1d5_h2_scan0.8RHF.out\n",
      "0.9\n",
      "overwrite output file: m1d5_h2_scan0.9RHF.out\n",
      "1.0\n",
      "overwrite output file: m1d5_h2_scan1.0RHF.out\n",
      "1.1\n",
      "overwrite output file: m1d5_h2_scan1.1RHF.out\n",
      "1.2\n",
      "overwrite output file: m1d5_h2_scan1.2RHF.out\n",
      "1.3\n",
      "overwrite output file: m1d5_h2_scan1.3RHF.out\n",
      "1.4\n",
      "overwrite output file: m1d5_h2_scan1.4RHF.out\n",
      "1.5\n",
      "overwrite output file: m1d5_h2_scan1.5RHF.out\n",
      "1.6\n",
      "overwrite output file: m1d5_h2_scan1.6RHF.out\n",
      "1.7\n",
      "overwrite output file: m1d5_h2_scan1.7RHF.out\n",
      "1.8\n",
      "overwrite output file: m1d5_h2_scan1.8RHF.out\n",
      "1.9\n",
      "overwrite output file: m1d5_h2_scan1.9RHF.out\n",
      "2.0\n",
      "overwrite output file: m1d5_h2_scan2.0RHF.out\n"
     ]
    }
   ],
   "source": [
    "listEneRHF = []\n",
    "#using the same distances array defined above\n",
    "for r in distances:\n",
    "    print(r)\n",
    "    listEneRHF.append(runThisDist(r,\"RHF\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "id": "ab5d1980",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj4AAAHFCAYAAADyj/PrAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAbgFJREFUeJzt3XlcVPX+x/HXsIMIiKigIDDu+567QuaSptUty7qZmtqittpq5dLmLc26Wb+0MrXbvlnaqqXkvqa5ryNugIrIoijr+f1BTCKggMhh4P18POZRc+Y7Z95zZjjz8Zzv+X4thmEYiIiIiFQCTmYHEBERESkrKnxERESk0lDhIyIiIpWGCh8RERGpNFT4iIiISKWhwkdEREQqDRU+IiIiUmmo8BEREZFKQ4WPiIiIVBoqfMq5efPmYbFY7DcPDw8CAwOJjIxk6tSpnDhxIt9zJk+ejMViMSFt8VgsFiZPnlzs56WmpjJ58mSioqLyPZa7vaKjo684X2WXkJDAkCFDqFmzJhaLhZtuuumqvl5ERATNmzcv8LH4+Ph835dNmzYxduxYWrRoQdWqValVqxbXXXcdS5cuLZU8ixYtYuDAgdSqVQs3Nzf8/f3p1asXn3zyCRkZGaXyGuVZWFgYw4cPN+W1LRYL48aNK/Cxr7/+GovFkufvf/jw4Xh7exe6Pm9v7zzvJTo6Os9+9cJb+/btS+ttXLGIiAgiIiIu2y4sLCzPe6hSpQpt27bl7bff5uLJGaKiorBYLHz99dcFrmvcuHH5fj8iIiIK3V7bt28v8fszi4vZAaRo5s6dS+PGjcnIyODEiROsXLmSV199lenTp/PFF19w3XXX2duOGjWKfv36mZi2aNasWUNwcHCxn5eamsqUKVMA8u0UBgwYwJo1awgKCiqNiJXaiy++yIIFC/jwww+pV68e/v7+ZkfK47PPPmP9+vXcc889tGrVirNnzzJr1ix69erF/Pnzufvuu0u0XsMwuOeee5g3bx79+/dnxowZhISEkJSUxLJlyxgzZgzx8fE8/PDDpfyOypcFCxbg4+Njdoyr6sEHH+TOO+/Ms+xSBVR51rVrV6ZPnw5ATEwMM2bM4MEHHyQ5OZkJEyZc8fqtViuffPJJvuX16tW74nWXNRU+DqJ58+Z5/iVyyy238Oijj9KtWzf+9a9/sW/fPmrVqgVAcHBwiQqKstapU6dSX2eNGjWoUaNGqa+3PMjIyMBiseDiUjZ/ttu3b6devXr8+9//LpX1GYbB+fPn8fT0LJX1Pfnkk/Ydfa7+/fvTtm1bXnjhhRIXPtOmTWPevHlMmTKFiRMn5nls4MCBPPnkk+zfv7/Eucu7c+fO4enpSZs2bcyOctXVrVv3quyHzODn55fnvVx33XXUrVuX2bNnl0rh4+npWWG2lU51ObC6devy+uuvk5KSwuzZs+3LCzrVtXTpUiIiIqhevTqenp7UrVuXW265hdTUVHubd999l1atWuHt7U3VqlVp3Lhxvj+Y7du3c+ONN1KtWjU8PDxo3bo18+fPz5ctMTGR8ePHY7VacXd3p2bNmvTv35/du3fb21x86uLkyZOMGTOGpk2b4u3tTc2aNbn22mtZsWKFvU10dLS9sJkyZYr9cGvuYezCTnV9+OGHtGrVCg8PD/z9/bn55pvZtWtXnja5h8v3799P//798fb2JiQkhPHjx5OWlnaJT+Ifn376KZ07d8bb2xtvb29at27NnDlz7I8Xdvrg4kPauYej//e//zF+/Hjq1KmDu7s7O3bswGKx5Flnrp9//hmLxcLChQvty/bt28edd95JzZo1cXd3p0mTJrzzzjuXfA+5pwF+++03du3aZd/GuacWEhISGDNmDHXq1MHNzQ2r1cqzzz6bbxvlnq6YNWsWTZo0wd3dvcDvSknVrFkz3zJnZ2fatWvHkSNHSrTOjIwMXn31VRo3bszzzz9fYJvAwEC6detmv1/c7TF37lwaNWqEp6cn7du3Z+3atRiGwbRp0wgPD8fb25trr702X3GVeypwxYoVdOrUCU9PT+rUqcPzzz9PVlZWnrZTpkyhY8eO+Pv74+PjQ9u2bZkzZ06+0x5hYWHccMMNfPvtt7Rp0wYPDw/70dSLv6u538nPPvuMZ599ltq1a+Pj48N1113Hnj178qzXMAxeeeUVQkND8fDwoH379ixZsqTIp27Ks/379zNixAgaNGiAl5cXderUYeDAgWzbti1Pu+Jur9dee82+vdq2bcvPP/98RTl9fHxo2LAhx48fv6L1VEQ64uPg+vfvj7OzM8uXLy+0TXR0NAMGDKB79+58+OGH+Pn5cezYMX755RfS09Px8vLi888/Z8yYMTz44INMnz4dJycn9u/fz86dO+3r2bNnD126dKFmzZq89dZbVK9enY8//pjhw4dz/PhxnnzySQBSUlLo1q0b0dHRPPXUU3Ts2JEzZ86wfPlyYmNjady4cYE5ExISAJg0aRKBgYGcOXOGBQsWEBERwe+//05ERARBQUH88ssv9OvXj5EjRzJq1CiASx7lmTp1KhMmTOCOO+5g6tSpnDp1ismTJ9O5c2c2bNhAgwYN7G0zMjIYNGgQI0eOZPz48SxfvpwXX3wRX1/ffP/6v9jEiRN58cUX+de//sX48ePx9fVl+/btHDp06JLPu5RnnnmGzp07M2vWLJycnAgJCaFNmzbMnTuXkSNH5mk7b948e4EJsHPnTrp06WIvkAMDA/n111956KGHiI+PZ9KkSQW+ZlBQEGvWrGHMmDEkJSXZD283bdqU8+fPExkZyYEDB5gyZQotW7ZkxYoVTJ06lS1btvDjjz/mWdd3333HihUrmDhxIoGBgQUWKxfLzMzMt+ziH/ZLPXfFihU0a9Ysz/KoqCgiIyOZNGnSJfuVbdy4kYSEBEaPHl2kfnLF3R4//PADmzdv5j//+Q8Wi4WnnnqKAQMGMGzYMGw2G2+//TZJSUk89thj3HLLLWzZsiVPjri4OIYMGcLTTz/NCy+8wI8//shLL73E6dOnefvtt+3toqOjue+++6hbty4Aa9eu5cEHH+TYsWP5vsd//vknu3bt4rnnniM8PJwqVapc8j1PmDCBrl278sEHH5CcnMxTTz3FwIED2bVrF87OzgA8++yzTJ06lXvvvZd//etfHDlyhFGjRpGRkUHDhg0vu11LqqDvzqVkZ2fne46zs/MlP/uYmBiqV6/Of/7zH2rUqEFCQgLz58+nY8eObN68mUaNGuVpX5TtNWXKFKZMmcLIkSO59dZbOXLkCKNHjyYrKyvf+ooqMzOTI0eOFLq9C3rvQL7i+OJ1XsjJyQknJwc8fmJIuTZ37lwDMDZs2FBom1q1ahlNmjSx3580aZJx4Uf79ddfG4CxZcuWQtcxbtw4w8/P75JZhgwZYri7uxuHDx/Os/z66683vLy8jMTERMMwDOOFF14wAGPJkiWXXB9gTJo0qdDHMzMzjYyMDKNXr17GzTffbF9+8uTJQp+bu70OHjxoGIZhnD592vD09DT69++fp93hw4cNd3d3484777QvGzZsmAEYX375ZZ62/fv3Nxo1anTJ92Kz2QxnZ2fj3//+9yXbhYaGGsOGDcu3vGfPnkbPnj3t95ctW2YARo8ePfK1feuttwzA2LNnj31ZQkKC4e7ubowfP96+rG/fvkZwcLCRlJSU5/njxo0zPDw8jISEhEtm7dmzp9GsWbM8y2bNmlXgNnr11VcNwFi8eLF9GWD4+vpe9nUufD3gkrdLfV8MwzCeffZZAzC+++67PMujoqIMZ2dnY8qUKZd8/ueff24AxqxZs4qUubjbIzAw0Dhz5ox92XfffWcARuvWrY3s7Gz78jfffNMAjK1bt9qX5W6f77//Ps9rjR492nBycjIOHTpUYMasrCwjIyPDeOGFF4zq1avneZ3Q0FDD2dk5z3fpwscu/K7mficv/lv68ssvDcBYs2aNYRj/fBdvv/32PO3WrFljAHm+54UBjLFjxxb42FdffWUAxrJly+zLcv92L3W78L0cPHiw0HaX229dLDMz00hPTzcaNGhgPProo/blRd1ep0+fNjw8PPLs4wzDMFatWlXk7RUaGmr079/fyMjIMDIyMoxDhw4Zo0ePNlxdXY0ffvghT9vcXJe7Xaiwv83L7e/KKwcs1eRixiUqdIDWrVvj5ubGvffey/z587HZbPnaXHPNNSQmJnLHHXfw/fffEx8fn6/N0qVL6dWrFyEhIXmWDx8+nNTUVNasWQPknHJp2LBhng7XRTVr1izatm2Lh4cHLi4uuLq68vvvv+c7LVVUa9as4dy5c/lOL4WEhHDttdfy+++/51lusVgYOHBgnmUtW7a87FGbJUuWkJWVxdixY0uUszC33HJLvmX//ve/cXd3Z968efZln332GWlpaYwYMQLIORLx+++/c/PNN+Pl5UVmZqb91r9/f86fP8/atWuLnWfp0qVUqVKFW2+9Nc/y3O178fa89tprqVatWpHXX69ePTZs2JDv9ttvv132uR988AEvv/wy48eP58Ybb8zzWM+ePcnMzLzsUbviKu72iIyMzHNEpUmTJgBcf/31eY4y5C6/+HtXtWpVBg0alGfZnXfeSXZ2dp6jvkuXLuW6667D19cXZ2dnXF1dmThxIqdOncp3JWjLli2LdRTm4tdv2bJlnqxr164lLS2N2267LU+7Tp06ERYWVuTXKS5PT88CvzsbNmwotF/Zww8/nK9tx44dL/k6mZmZvPLKKzRt2hQ3NzdcXFxwc3Nj3759Be6nLre91qxZw/nz5/P1pevSpQuhoaFFfv8//fQTrq6uuLq6Ehoayvvvv8/MmTMZMGBAge1fffXVArfVxZ9broL+Nl988cUi5ytPdKrLwZ09e5ZTp07RokWLQtvUq1eP3377jddee42xY8dy9uxZrFYrDz30kP3KlKFDh5KZmcn777/PLbfcQnZ2Nh06dOCll16id+/eAJw6darAq6Vq165tfxxy+urkHmIvjhkzZjB+/Hjuv/9+XnzxRQICAnB2dub5558vceGTm6mw3EuWLMmzzMvLCw8PjzzL3N3dOX/+/CVf5+TJkwCl3qm8oNz+/v4MGjSIjz76iBdffBFnZ2fmzZvHNddcYz/Fc+rUKTIzM5k5cyYzZ84scN0FFbeXc+rUKQIDA/OdCqhZsyYuLi727X2p/JeS2x+kuFnnzp3Lfffdx7333su0adOK9ZoXyv3eHjx4sEjti7s9Lr4yzs3N7ZLLL/7e5V7AcKHAwEB7FoD169fTp08fIiIieP/99wkODsbNzY3vvvuOl19+mXPnzuV5fnE/o+rVq+e57+7uDmBfb26OgrIWtKwgzs7OhZ7ezD3d4urqmme5k5NToZeiF3Y6Jjg4uNiXrz/22GO88847PPXUU/Ts2ZNq1arh5OTEqFGj8m1bKPr2yv0cL1TQssJ069aNN954g6ysLPbt28fzzz/PuHHjaNasWZ4+abmsVmuB772wbgOF/W06IhU+Du7HH38kKyvrsh0Gu3fvTvfu3cnKymLjxo3MnDmTRx55hFq1ajFkyBAARowYwYgRIzh79izLly9n0qRJ3HDDDezdu5fQ0FCqV69ObGxsvnXHxMQAEBAQAOT84Rw9erTY7+Xjjz8mIiKCd999N8/ylJSUYq8rV+5Op7DcuZmvVO7O4ujRo/mOiF3Iw8OjwI7S8fHxBWYprK/BiBEj+Oqrr1iyZAl169Zlw4YNebZbtWrVcHZ2ZujQoYUehQoPD7/keypI9erVWbduHYZh5Ml24sQJMjMz872HshhPau7cuYwaNYphw4Yxa9asK3rN9u3b4+/vz/fff8/UqVMvu67ibo8rVVBH1bi4OHsWgM8//xxXV1d++OGHPEX8d999V+A6S/szys1RWNaiHPWpVasWx44dK/Cx3OVFLaJK28cff8zdd9/NK6+8kmd5fHw8fn5+xV5f7vbK/RwvVNTtBeDr62svTDp27EjHjh1p1aoVY8aMYcuWLY7ZF+cq0ZZwYIcPH+bxxx/H19eX++67r0jPcXZ2pmPHjvYre/788898bapUqcL111/Ps88+S3p6Ojt27ACgV69eLF261F7o5Proo4/w8vKyX+p4/fXXs3fv3mIPJGexWOz/Gsq1detW+ym0XBf/i+lSOnfujKenJx9//HGe5UePHrWfuisNffr0wdnZOV/RdrGwsDC2bt2aZ9nevXvzXeVRlNerU6cOc+fOZe7cuXh4eHDHHXfYH/fy8iIyMpLNmzfTsmVL2rdvn+928b9Ei6JXr16cOXMm34/oRx99ZH+8LM2bN49Ro0Zx11138cEHH1zxj7irqytPPfUUu3fvLvQw/okTJ1i1ahVQ9tsjJSUlz1V7kHMloZOTEz169ACwD3mQ23EWcv5W/ve//5VqlsJ07NgRd3d3vvjiizzL165dW+SO/tdddx3Lli2zH0nNZRgGX331FWFhYdSvX7/UMhdHQfupH3/8sdBC7XI6deqEh4dHvjFyVq9efUUXRjRo0IAnn3ySbdu25fssKjsd8XEQ27dvt/fROHHiBCtWrGDu3Lk4OzuzYMGCS17VNGvWLJYuXcqAAQOoW7cu58+f58MPPwSw98MZPXo0np6edO3alaCgIOLi4pg6dSq+vr506NAByLna6ocffiAyMpKJEyfi7+/PJ598wo8//shrr72Gr68vAI888ghffPEFN954I08//TTXXHMN586d448//uCGG24gMjKywJw33HADL774IpMmTaJnz57s2bOHF154gfDw8DxXE1StWpXQ0FC+//57evXqhb+/PwEBAQX+y8jPz4/nn3+eCRMmcPfdd3PHHXdw6tQppkyZgoeHR6FXNhVXWFgYEyZM4MUXX+TcuXPccccd+Pr6snPnTuLj4+2XCA8dOpS77rqLMWPGcMstt3Do0CFee+21Yo895OzszN13382MGTPw8fHhX//6l3375/rvf/9Lt27d6N69Ow888ABhYWGkpKSwf/9+Fi1aVKIRju+++27eeecdhg0bRnR0NC1atGDlypW88sor9O/fv0T9ukrqq6++YuTIkbRu3Zr77ruP9evX53m8TZs29h+oP/74g169ejFx4sTL9vN54okn2LVrF5MmTWL9+vXceeed9gEMly9fznvvvceUKVPo2rVrmW+P6tWr88ADD3D48GEaNmzITz/9xPvvv88DDzxgP003YMAAZsyYwZ133sm9997LqVOnmD59er4f66vF39+fxx57jKlTp1KtWjVuvvlmjh49ypQpUwgKCirSkYeJEyeyaNEiOnbsyNNPP02DBg2Ii4vj/fffZ8OGDXz55Zdl8E4KdsMNNzBv3jwaN25My5Yt2bRpE9OmTSvxae5q1arx+OOP89JLLzFq1CgGDx7MkSNHmDx5crFOdRXk8ccfZ9asWUyZMoXbbrstTzFcqZnbt1ouJ/cqpdybm5ubUbNmTaNnz57GK6+8Ypw4cSLfcy6+qmvNmjXGzTffbISGhhru7u5G9erVjZ49exoLFy60t5k/f74RGRlp1KpVy3BzczNq165t3HbbbXmuKjEMw9i2bZsxcOBAw9fX13BzczNatWplzJ07N1+G06dPGw8//LBRt25dw9XV1ahZs6YxYMAAY/fu3fY2XHSVTlpamvH4448bderUMTw8PIy2bdsa3333nTFs2DAjNDQ0z/p/++03o02bNoa7u3ueKzYuvqor1wcffGC0bNnScHNzM3x9fY0bb7zR2LFjR542w4YNM6pUqXLZ7XkpH330kdGhQwfDw8PD8Pb2Ntq0aZNn+2RnZxuvvfaaYbVaDQ8PD6N9+/bG0qVLC72q66uvvir0tfbu3XvZK1EOHjxo3HPPPUadOnUMV1dXo0aNGkaXLl2Ml1566bLvpaCrugzDME6dOmXcf//9RlBQkOHi4mKEhoYazzzzjHH+/Pk87bjElTnFeT3DKPhKvstdyXPhdyB3e17uqrALff/998aAAQOMGjVqGC4uLka1atWMyMhIY9asWUZaWpq93ZVsj9yri6ZNm5ZneUGff+72iYqKMtq3b2+4u7sbQUFBxoQJE4yMjIw8z//www+NRo0aGe7u7obVajWmTp1qzJkzJ992CQ0NNQYMGFDg+y/sqq6Lv5O57+Hi7/lLL71kBAcHG25ubkbLli2NH374wWjVqlW+q5cKs2/fPuOuu+6yb1c/Pz+jT58+xu+//56vbWF/u7mqVKlS4FVdF2/3ojh9+rQxcuRIo2bNmoaXl5fRrVs3Y8WKFUX+Gy5se02dOtUICQmxb69FixblW2dhLvU5vvPOOwZgzJ8//5K5co0dO7bAq7oK+9t0RBbDuMwlQSIiYrqIiAji4+Mdcm4kyOkw3rhxYyZNmlQqIwmLlJROdYmISKn666+/+Oyzz+jSpQs+Pj7s2bOH1157DR8fn3wDb4qUNRU+IiJSqqpUqcLGjRuZM2cOiYmJ+Pr6EhERwcsvv2za1VgiuXSqS0RERCoNXc4uIiIilYYKHxEREak0VPiIiIhIpaHOzRfJzs4mJiaGqlWrlslw+yIiInLlDMMgJSWF2rVrX3KgTBU+F4mJibnkXEsiIiJSfh05cuSSI2mr8LlI1apVgZwN5+PjY3IaERERKYrk5GRCQkLsv+OFUeFzkdzTWz4+Pip8REREHMzluqmoc7OIiIhUGip8REREpNJQ4SMiIiKVhvr4iIiImCQrK4uMjAyzYzgEV1dXnJ2dr3g9KnxERETKmGEYxMXFkZiYaHYUh+Ln50dgYOAVjbOnwkdERKSM5RY9NWvWxMvLSwPmXoZhGKSmpnLixAkAgoKCSrwuFT4iIiJlKCsry170VK9e3ew4DsPT0xOAEydOULNmzRKf9lLnZhERkTKU26fHy8vL5CSOJ3ebXUm/KBU+IiIiJtDpreIrjW2mwkdEREQqDRU+IiIiUmmo8CkjRzfEsmzGZo5uiDU7ioiISIlERETwyCOP5Fv+3Xff2U9DzZs3Dz8/vwKfb7FY+O677/Lcv/jWrVu3q5D8H7qqqwzMGb6Ce+d3IZsgnMjivWErGDmvu9mxRERETDd37lz69etnv+/m5nZVX09HfK6yoxti/y56ci67y8aZ++Z31pEfEREpHUePwrJlOf91QLmDEube/P39r+rr6YjPVbZvRRzZ5B1oKQsX9q86TnCHkg/AJCIiFYhhQGpq8Z83fz48+CBkZ4OTE8ycCcOGFf35Xl5Qya4uU+FzlTXoHogTWfYjPgDOZFK/ay0TU4mISLmSmgre3le2juxsGDs251ZUZ85AlSpX9roFSEpKwruI7+eOO+7IMxjhxx9/zE033VTqmXKp8LnKgjsE8d6wFYye3xUDJyxkM3vYGoI7qI+PiIhUTFWrVuXPP//Mt7xBgwb5lr3xxhtcd9119vtXMh1FUajwKQMj53Un4eQynvwpkh6+f6ljs4iI5OXllXP0pTiOHYMmTXKO9ORydoadO6FOnaK/bjH4+PiQlJSUb3liYiI+Pj72+05OTtSvX79I6wwMDCxy29KgwqeMtO/lBz9BTKqf2VFERKS8sViKf8qpYUN47z247z7IysopembPzll+lTRu3Jiff/453/INGzbQqFGjq/a6pUlXdZURa8caABzKqE1WepbJaUREpEIYORKio3Ou6oqOzrl/FY0ZM4YDBw4wduxY/vrrL/bu3cs777zDnDlzeOKJJ67qa5cWHfEpI8EdgnAhg3TcidkcQ0jH2mZHEhGRiiA4OOdWBsLCwlixYgXPPvssffr04fz58zRs2JB58+YxePDgMslwpRzmiM/LL79Mly5d8PLyKnREyA0bNtCrVy/8/PyoVq0affr0YcuWLWWaszDObs6EuR4DwLb2hMlpRERESqZdu3b88ssvHD9+nKSkJDZs2MCQIUPsjw8fPpzExMQCn2sYRp4rti6+XxYcpvBJT09n8ODBPPDAAwU+npKSQt++falbty7r1q1j5cqV+Pj40Ldv3yuavr40hVc9BYDtrxSTk4iIiFRODnOqa8qUKUDOHCAF2bNnD6dPn+aFF14gJCQEgEmTJtGyZUsOHz5MvXr1yipqoay1zkICHNyvPj4iIiJmcJgjPpfTqFEjAgICmDNnDunp6Zw7d445c+bQrFkzQkNDC31eWloaycnJeW5XizUs55JD2xHXq/YaIiIiUrgKU/hUrVqVqKgoPv74Yzw9PfH29ubXX3/lp59+wsWl8ANbU6dOxdfX137LPVp0NVibuANgi/e5TEsRERG5GkwtfCZPnlzglPQX3jZu3FikdZ07d4577rmHrl27snbtWlatWkWzZs3o378/586dK/R5zzzzDElJSfbbkSNHSuvt5RPexg8AW6qmqxARETGDqX18xo0bl6cneEHCwsKKtK5PP/2U6Oho1qxZg5OTk31ZtWrV+P777wt9HXd3d9zd3YuVu6SsXXOG4T6eXZOzJ85SpWbpz48iIiIihTO18AkICCAgIKBU1pWamoqTkxOWC2aZzb2ffeFw3iaqFu6HnyWRRMOP6DWxNLux7IboFhEREQfq43P48GG2bNnC4cOHycrKYsuWLWzZsoUzf89t0rt3b06fPs3YsWPZtWsXO3bsYMSIEbi4uBAZGWly+n9YPWIBsG1MMDmJiIhI5eMwhc/EiRNp06YNkyZN4syZM7Rp04Y2bdrY+wA1btyYRYsWsXXrVjp37kz37t2JiYnhl19+ueozvRZHuH8iALYdhfc7EhERkavDYcbxmTdvXqFj+OTq3bs3vXv3LptAJWStnQbHwGYzO4mIiEjl4zBHfCoKa/2cTX4wzsPkJCIiIsUzfPhw+1XXLi4u1K1blwceeIDTp0/b24SFhfHmm2/me+7kyZNp3bp1nvsFXc3922+/XdX34DBHfCoKa4sq8BnYEv3NjiIiIlJs/fr1Y+7cuWRmZrJz507uueceEhMT+eyzz4q9rmbNmuUrdPz9r+7vowqfMma9JucqNltaHYxsA4uT5TLPEBERKdzRo7BvHzRoUDaTtLu7uxMYGAhAcHAwt99++2W7ohTGxcXFvq6yosKnjNXtGISFbM7hxfHtJwhsWdPsSCIiYjLDgNTU4j9v/nx48EHIzgYnJ5g5E4YNK/rzvbzAcgX//rbZbPzyyy+4ujrOVEwqfMqYm7cbIc5HOZwVzMG1x1X4iIgIqang7X1l68jOhrFjc25FdeYMVCnmWLo//PAD3t7eZGVlcf78eQBmzJiRp81TTz3Fc889l2dZeno6TZs2zbNs27ZteF/wxps2bcr69euLF6iYVPiYwFr1JIcTg7FtSaaz2WFERESKITIyknfffZfU1FQ++OAD9u7dy4MPPpinzRNPPMHw4cPzLHvrrbdYvnx5nmWNGjVi4cKF9vtlMZOCCh8TWGucISoRbHsyzI4iIiLlgJdXztGX4jh2DJo0yTnSk8vZGXbuhDp1iv66xVWlShXq18+ZeeCtt94iMjKSKVOm8OKLL9rbBAQE2NvkKqjTspubW752V5sKHxOE182CfWA77Gx2FBERKQcsluKfcmrYEN57D+67D7Kycoqe2bNzlpelSZMmcf311/PAAw9Qu3btsn3xEtA4PiawNsrpBGY7WdXkJCIi4shGjoToaFi2LOe/I0eWfYaIiAiaNWvGK6+8UvYvXgIqfExgbeMLwMEzNUxOIiIiji44GCIiyuZS9sI89thjvP/++xw5csS8EEVkMQzDMDtEeZKcnIyvry9JSUn4+Phcldc4seMktZrXyLmsPSkDd5+r35lLRETKh/Pnz3Pw4EHCw8Px8NAo/sVxqW1X1N9vHfExQY0mAVThDAZOHFoTY3YcERGRSkOFjwksThbC3WMBsG04ZXIaERGRykOFj0ms1RIAOLijBEN1ioiISImo8DGJNShntEvb/uzLtBQREZHSosLHJFZrzn9tMerYJiJSGenaouIrjW2mwsck4U09AbCd9jM3iIiIlKncCT1TSzIraSWXu82uZFJUjdxsEmv7nKG7beeCMLINLE5XMD2uiIg4DGdnZ/z8/Dhx4gQAXl5eWK5kivRKwDAMUlNTOXHiBH5+fjg7l3zmAxU+JgnrkjOsdzK+nD54Gv961UxOJCIiZSUwMBDAXvxI0fj5+dm3XUmp8DGJV4AXQU5xxGYHYlsVq8JHRKQSsVgsBAUFUbNmTTIyNGF1Ubi6ul7RkZ5cKnxMFF7lBLEpgdg2J9L+brPTiIhIWXN2di6VH3MpOnVuNpE1IBmAg7vTTU4iIiJSOajwMZE1JBMAW7Q+BhERkbKgX1wTWRvmnGm0nahichIREZHKQYWPiaytqgJgSw4wOYmIiEjloMLHROEdawJwKLMOmeczTU4jIiJS8anwMVHtNrVwI40sXDi6Mc7sOCIiIhWeCh8TObk4Ee52DADbupMmpxEREan4VPiYzOp7CgDb1jMmJxEREan4VPiYLLxWzoRrtn1ZJicRERGp+FT4mMwaZgBw8JibyUlEREQqPhU+JrM29QDAdsrH5CQiIiIVnwofk1nb5UxOaku9stlmRURE5PJU+JgsvEsQAPFGACkxKSanERERqdhU+JjMJ9iH6pacK7sOro41OY2IiEjFpsKnHLB65QxeaNt02uQkIiIiFZsKn3LA6p8EgG3nOZOTiIiIVGwqfMoBa3A6ALaDFpOTiIiIVGwqfMqB8Po5H8PBOC+Tk4iIiFRsKnzKAWsLbwBsSf4mJxEREanYHKLwiY6OZuTIkYSHh+Pp6Um9evWYNGkS6enpedodPnyYgQMHUqVKFQICAnjooYfytSmPrB1rAHAwvQ7ZmdkmpxEREam4XMwOUBS7d+8mOzub2bNnU79+fbZv387o0aM5e/Ys06dPByArK4sBAwZQo0YNVq5cyalTpxg2bBiGYTBz5kyT38GlhVwThDOZpOFB7JZY6rQPMjuSiIhIhWQxDMMwO0RJTJs2jXfffRebzQbAzz//zA033MCRI0eoXbs2AJ9//jnDhw/nxIkT+PgUbUqI5ORkfH19SUpKKvJzSoPV9TAHM+uyfOZfdB/XqsxeV0REpCIo6u+3Q5zqKkhSUhL+/v/0iVmzZg3Nmze3Fz0Affv2JS0tjU2bNpkRsVisVeMBOLhVozeLiIhcLQ5Z+Bw4cICZM2dy//3325fFxcVRq1atPO2qVauGm5sbcXFxha4rLS2N5OTkPDczWGudAcC2N9OU1xcREakMTC18Jk+ejMViueRt48aNeZ4TExNDv379GDx4MKNGjcrzmMWSfxwcwzAKXJ5r6tSp+Pr62m8hISGl8+aKyRqa06nZdsQhul2JiIg4JFN/ZceNG8eQIUMu2SYsLMz+/zExMURGRtK5c2fee++9PO0CAwNZt25dnmWnT58mIyMj35GgCz3zzDM89thj9vvJycmmFD/WJm7wK9hOll2/IhERkcrG1MInICCAgICAIrU9duwYkZGRtGvXjrlz5+LklPdgVefOnXn55ZeJjY0lKCjnqqjFixfj7u5Ou3btCl2vu7s77u7uJX8TpSS8tR8AB8/WMDeIiIhIBeYQfXxiYmKIiIggJCSE6dOnc/LkSeLi4vL03enTpw9NmzZl6NChbN68md9//53HH3+c0aNHl+nVWSVl7RIIQEx2EOcSNGeXiIjI1eAQHUoWL17M/v372b9/P8HBwXkey70a39nZmR9//JExY8bQtWtXPD09ufPOO+3j/JR3/vWq4UMSyfgSvTqGJjfUMzuSiIhIheOw4/hcLWaN4wPQxms3W8415odJGxgwuUOZvraIiIgjq/Dj+FRE4X6JABzckWpuEBERkQpKhU85Yq19HgCbTQfhRERErgYVPuWItX7OeEO2GE+Tk4iIiFRMKnzKEWvzKgDYEquZnERERKRiUuFTjoS3y5l7zHa+Nka2TneJiIiUNhU+5Uho59pYyOYs3sTvOWV2HBERkQpHhU854uHnQR3nnEEZbasLn1hVRERESkaFTzlj9T4BgG2LObPEi4iIVGQqfMoZa40UAGy7001OIiIiUvGo8ClnwkOyADh42NnkJCIiIhWPCp9yxtrIFQDbCW+Tk4iIiFQ8KnzKGWvrnPlFbCk1TE4iIiJS8ajwKWesnWsBcCQriIzUDJPTiIiIVCwqfMqZWs1r4Ekq2ThzeG2M2XFEREQqFBU+5YzFyUK4e07BY1sfb3IaERGRikWFTzlk9U0AwLbtrMlJREREKhYVPuWQNegcALb92SYnERERqVhU+JRD4eE5/z0Y425uEBERkQpGhU85ZG3qAYAtwdfkJCIiIhWLCp9yyNreHwDbuSCTk4iIiFQsKnzKofCutQE4bVTj9MFEc8OIiIhUICp8yqEqNatQyylnlvaDq2NNTiMiIlJxqPApp8K9jgNwcEuiuUFEREQqEBU+5ZS1ejIAtp1pJicRERGpOFT4lFPWkJx5umzR+ohERERKi35VyylrA2cAbMe9TE4iIiJScajwKafCW1YF4GBydZOTiIiIVBwqfMopa8caAERn1CErPcvkNCIiIhWDCp9yqk67QFxJJwM3jm2KMzuOiIhIhaDCp5xydnMmzPUYALZ1J01OIyIiUjGo8CnHrD6nADi4NcXkJCIiIhWDCp9yLLxWKgC2ferjIyIiUhpU+JRj1rBsAGxHXE1OIiIiUjGo8CnHrE3dAbCd8jE5iYiISMWgwqccs7atBoDtbKDJSURERCoGFT7lWHiXIABOGDU4e+KsyWlEREQcnwqfcswv1JdqltMAHFwVY3IaERERx6fCp5yzesYCYNuYYHISERERx6fCp5yz+icBYNt53uQkIiIijk+FTzkXXjsNgIMHTQ4iIiJSAajwKees9XM+Ilusp8lJREREHJ9DFD7R0dGMHDmS8PBwPD09qVevHpMmTSI9Pd3e5q+//uKOO+4gJCQET09PmjRpwn//+18TU5cOa4sqANiS/E1OIiIi4vhczA5QFLt37yY7O5vZs2dTv359tm/fzujRozl79izTp08HYNOmTdSoUYOPP/6YkJAQVq9ezb333ouzszPjxo0z+R2UnPWaAAAOptXGyDawOFlMTiQiIuK4LIZhGGaHKIlp06bx7rvvYrPZCm0zduxYdu3axdKlS4u83uTkZHx9fUlKSsLHx/wRkzNSM/Co4kQ2zsT+dYLAljXNjiQiIlLuFPX32yFOdRUkKSkJf/9Ln/4pSpu0tDSSk5Pz3MoTVy9XQpz/vqR9zXGT04iIiDg2hyx8Dhw4wMyZM7n//vsLbbNmzRq+/PJL7rvvvkuua+rUqfj6+tpvISEhpR33ilmrngTAtqV8FWUiIiKOxtTCZ/LkyVgslkveNm7cmOc5MTEx9OvXj8GDBzNq1KgC17tjxw5uvPFGJk6cSO/evS+Z4ZlnniEpKcl+O3LkSKm9v9JirXkGANueDJOTiIiIODZTOzePGzeOIUOGXLJNWFiY/f9jYmKIjIykc+fOvPfeewW237lzJ9deey2jR4/mueeeu2wGd3d33N3di5W7rFlDs2AvHDzibHYUERERh2Zq4RMQEEBAQECR2h47dozIyEjatWvH3LlzcXLKf7Bqx44dXHvttQwbNoyXX365tOOaJryRGywB28mqZkcRERFxaA5xOXtMTAwRERHUrVuX6dOnc/LkSftjgYGBQE7RExkZSZ8+fXjssceIi4sDwNnZmRo1apiSu7RYW+f0Tred0RVdIiIiV8IhCp/Fixezf/9+9u/fT3BwcJ7Hcq/G/+qrrzh58iSffPIJn3zyif3x0NBQoqOjyzJuqbN2ySnujmUFcj7xPB5+HiYnEhERcUwOO47P1VLexvEBMLINfJzPcIaq7P7JRqPrrWZHEhERKVcq/Dg+lYnFyYLVIwaAg5sSTE4jIiLiuFT4OIhwv0QAbNvPmhtERETEganwcRDW2ucAsO3XmUkREZGSUuHjIKzWnMlJbTHq2CwiIlJSKnwchLW5FwAHE/3MDSIiIuLAVPg4iPB2OZOt2s4FYWTrdJeIiEhJqPBxEGFdagOQjC8JB06bnEZERMQxqfBxEJ7+ntR2igXAtjrO5DQiIiKOSYWPA7F6nwDg4JZEc4OIiIg4KBU+DiQ8IAUA2650k5OIiIg4JhU+DsQakgmA7ZA+NhERkZLQL6gDsTbMmVPWdtzb5CQiIiKOSYWPA7G2zpl0zZYSYHISERERx6TCx4FYO9UE4HBmbTLPZ5qcRkRExPGo8HEggS1r4s55snDhyPpYs+OIiIg4HBU+DsTJxYlwt2MA2NadNDmNiIiI41Hh42CsvgkA2LadMTmJiIiI41Hh42CsgakAHNyfbXISERERx6PCx8GEh+dMUGo76mZyEhEREcejwsfBWJt6AmBL8DU5iYiIiONR4eNgrO2qAWBLDTQ5iYiIiONR4eNgwrsEAXDKqE7y0WST04iIiDgWFT4OpmrtqgRY4gE4uFpj+YiIiBSHCh8HZPWKA8C26bTJSURERByLCh8HZK2ec4rLtvO8yUlEREQciwofB2QNTgfg4CGLyUlEREQciwofB2Rt4AyALc7L5CQiIiKORYWPAwpv4Q2ALam6yUlEREQciwofB2TtWAOAg+l1yM7U1BUiIiJFpcLHAQW3D8SFDNJxJ2bzcbPjiIiIOAwVPg7IxcOFUJcYAA6uO2FyGhEREcehwsdBhfvkDGJo+yvF5CQiIiKOw6WkT0xLS2P9+vVER0eTmppKjRo1aNOmDeHh4aWZTwphrXkWEsC2N9PsKCIiIg6j2IXP6tWrmTlzJt999x3p6en4+fnh6elJQkICaWlpWK1W7r33Xu6//36qVq16NTILYA3Lht1gO1Li2lVERKTSKdaprhtvvJFbb72VOnXq8Ouvv5KSksKpU6c4evQoqamp7Nu3j+eee47ff/+dhg0bsmTJkquVu9KzNnEH4GC8j8lJREREHEexDhf06dOHr776Cjc3twIft1qtWK1Whg0bxo4dO4iJiSmVkJKfta0fALazNc0NIiIi4kAshmEYxXlCZmYmLi4V9/RKcnIyvr6+JCUl4eNTfo+mJBw4TfX61QA4ezIVrwCN4iwiIpVXUX+/i31VV1BQEI8//jg7d+68ooByZaqF++FLEgDRq3VkTUREpCiKXfg89thjLFq0iBYtWtC5c2fmzJnDmTNnrkY2uQSLkwWr599j+fx52uQ0IiIijqHYhc8zzzzDnj17iIqKonHjxjzyyCMEBQUxYsQIVq1adTUySiGs/okA2LanmhtERETEQZR4AMPu3bszd+5c4uLiePPNN9m/fz/du3enUaNGvPbaa6WZkejoaEaOHEl4eDienp7Uq1ePSZMmkZ6eXmD7U6dOERwcjMViITExsVSzlCfhQWkA2GwmBxEREXEQVzxyc5UqVRg5ciQrVqxg0aJFxMfH88wzz5RGNrvdu3eTnZ3N7Nmz2bFjB2+88QazZs1iwoQJBbYfOXIkLVu2LNUM5ZG1fs7HZ4v1MDmJiIiIY7jiy7NSU1P54osvmDt3LqtWraJevXo88cQTpZHNrl+/fvTr189+32q1smfPHt59912mT5+ep+27775LYmIiEydO5Oeffy7VHOWNtXnOlVwHE6uZnERERMQxlLjwWbFiBXPnzuXrr78mKyuLW2+9lZdeeokePXqUZr5CJSUl4e/vn2fZzp07eeGFF1i3bh22Ip7/SUtLIy0tzX4/OTm5VHNeTdZrAgCwna+NkW1gcbKYnEhERKR8K/aprldeeYWGDRsSERHBjh07mDZtGrGxscyfP7/Mip4DBw4wc+ZM7r//fvuytLQ07rjjDqZNm0bdunWLvK6pU6fi6+trv4WEhFyNyFdF3Y5BWMjmLN6c3BVvdhwREZFyr9iFzxtvvMGAAQP466+/WLduHffdd1+JB/qbPHkyFovlkreNGzfmeU5MTAz9+vVj8ODBjBo1yr78mWeeoUmTJtx1113FyvDMM8+QlJRkvx05cqRE78UM7j7uBDvHAmBbc9zkNCIiIuVfsUduzsjIwNXVNc+y8+fP4+FR/A628fHxxMdf+khFWFiYfd0xMTFERkbSsWNH5s2bh5PTP3Vb69at2bZtGxZLzukewzDIzs7G2dmZZ599lilTphQpk6OM3Jwrwm8LfyS15pMxq7jzna5mxxERETFFUX+/i93HJ7foyc7O5uWXX2bWrFkcP36cvXv3YrVaef755wkLC2PkyJGXXVdAQAABAQFFet1jx44RGRlJu3btmDt3bp6iB+Cbb77h3Llz9vsbNmzgnnvuYcWKFdSrV68Y79CxWGum8EcSHNybYXYUERGRcq/El7O/9NJLzJs3j9deey3PpKUtWrTggw8+KJVwuWJiYoiIiCAkJITp06dz8uRJ4uLiiIuLs7epV68ezZs3t9/Cw8MBaNKkCTVrVtyJPK11swCwHXI2OYmIiEj5V+Kruj766CPee+89evXqlaeTccuWLdm9e3ephMu1ePFi9u/fz/79+wkODs7zWDHP1FU44Q1d4XewnfQ2O4qIiEi5V+IjPseOHaN+/fr5lmdnZ5ORUbqnXYYPH45hGAXeChMREYFhGPj5+ZVqlvLG2jrnPKYtpYbJSURERMq/Ehc+zZo1Y8WKFfmWf/XVV7Rp0+aKQknRWbsEAnA0K4j0MwVP4SEiIiI5Snyqa9KkSQwdOpRjx46RnZ3Nt99+y549e/joo4/44YcfSjOjXELNpgF4cZZUqnB43VHq9wo1O5KIiEi5VeIjPgMHDuSLL77gp59+wmKxMHHiRHbt2sWiRYvo3bt3aWaUS7A4WQh3jwHAtl6DGIqIiFzKFc3V1bdvX/r27VtaWaSErH4J7DgOtm1nzY4iIiJSrl3x7OyXUtmvuCor1qDzABw8kG1yEhERkfKtWIVPkyZN+PTTT0lPv3Qn2n379vHAAw/w6quvXlE4KRqrNee/tmPu5gYREREp54p1quudd97hqaeeYuzYsfTp04f27dtTu3ZtPDw8OH36NDt37mTlypXs3LmTcePGMWbMmKuVWy5gbeYJ34Itwc/sKCIiIuVasQqfa6+9lg0bNrB69Wq++OILPv30U6Kjozl37hwBAQG0adOGu+++m7vuuqvCj59TnoS38wfAdj7I5CQiIiLlW7EnKa3oHG2SUoDU+FSq1PACIMGWSLVwP3MDiYiIlLGi/n4Xu3Nz//79SUpKst9/+eWXSUxMtN8/deoUTZs2Le5q5Qp4BXgR6HQcgIOrY01OIyIiUn4Vu/D59ddfSUtLs99/9dVXSUhIsN/PzMxkz549pZNOisxaJafwsf2ZaG4QERGRcqzYhc/FZ8Z0pqx8CK+eAoBtV9plWoqIiFReV3UcHyk71pCciWFt0fpIRUREClPsX0mLxYLFYsm3TMxlbZhzgd7BE14mJxERESm/ij1lhWEYDB8+HHf3nMHyzp8/z/3330+VKlUA8vT/kbJjbVUVAFtygMlJRERECnZ0Qyz7VsTRoHsgwR3MGYKl2IXPsGHD8ty/66678rW5++67S55ISiT8mhoARGfUISs9C2c3Z5MTiYiI/GPO8BXcO78L2QThRBbvDVvByHndyzyHxvG5iCOO4wOQlZ6Fl3sm6bgTvfIooV2DzY4kIiIC5BzpCb2mJtn8849yZzKJXn+y1I78XLVxfKR8cnZzJsw1BoCDG+JNTiMiIgKn9iUw89Y/iOx6Pk/RA5CFC/tXHS/zTMU+1SXll9U3nr3x4dj+SiHC7DAiIlIpZaVnsfg/fzL3/Uy+P9qWdHr+/YgB/HMxlDOZ1O9aq8zzqfCpQKyBqRAPtn1ZZkcREZFKZu+vB5k35RDz1zUiJruDfXlrz93c0/842dkG4xd0IwsXnMlk9rA1BHco+z4+KnwqkPBQA7aD7air2VFERKQSSIlJ4atnt/DhN76sSmkJhAPgb0ngrpbbGPFMIK1vbww0BuCWDbHsX3Wc+l1rmVL0gAqfCsXa1AN+hIOnHKdTtoiIOBYj22DFO1uZ+1YyX+1vw1lyChgnsuhX409GDM1k4KS2uPv0zPfc4A5Bpl3GnkuFTwVibVcNANvZQJOTiIhIRXNkXQwfPbeXuVHhHMhsZV/e0PUgI649xNCXGlGnfYdLrKF8UOFTgYR3zil4Thg1OBN3Bu9Ab5MTiYiIIzufeJ7vJ/3Jh5+4seRUWwxqA+BNCrc33MI9j/nReXRzLE7hJictOhU+FYhvXV/8LQkkGP4cXBVDi1samh1JREQcwIUjKtdpF8ifn+7mw1dP8OmOViQaXeztevpu4Z7bznDLS22oUtOcPjpXSoVPBWP1jCMh1R/bptO0uMXsNCIiUt5dOKKyhWxqO8VyLLsJ0ASAEOdjDO+yj2GTw6l3bWtTs5YGFT4VjLV6IhtTYfHC87S7Odb0TmQiIlJ+HV5zjNHzu2L8PZ6xgRPHsmvjxnluCd3EiPs9uPax1ji71TE5aenRyM0VTMq5nEvZ/29HT0Kvqcmc4StMTiQiIuXN8e0nmdo3imu6utiLngt9/dxWPo3uSu+n21W4uR9V+FQgRzfE8kt8O/v9bJy5b35njm6INTGViIiUB9mZ2fw+7U9uC1lDcAs/JiyO4LhRi5wRlf/hTCZtBoWYE7IM6FRXBbJvRRwGeU9t5c6FolNeIiKVU/yeU8wbv43Zi8PYn9HWvryz9zbuG5JE6pksHvy8q+kjKpcVFT4VSIPugTiRlW/2WzPmQhEREfPkDjI4+/UzfH2oPel/z+BYlWSGNt/MfZOCaHlrC3v7gY+ZP6JyWVHhU4EEdwjivWEr7B3VLGRX+MpdRET+cfpgIh+N/4vZPwazK/2fQQbbee3k/ltPMeTVNngHls8RlcuKCp8KZuS87sByRs3vQbBTLPd82M3sSCIichUZ2QZrP9jO7NeS+OJAO87/PRu6F2e5s9Em7psQQPu7m5qcsvywGIZhXL5Z5ZGcnIyvry9JSUn4+DjmnFdnT5wloJYT5/Fk27f7aH5zA7MjiYhIKUs+mswnT2xh1ne12Hq+kX15C4+93H9jLP9+rTW+dX1NTFi2ivr7rau6KqAqNatwXc1tACx895jJaUREpDRt+ngXoxsvp3aIE2M+78HW843w4BzD6q1k9ext/HW2AWM+71mpip7i0KmuCmpg7/P88AksWlOdCWaHERGRYrtwGgm/kKp8/tRmZn1dnU2pTckdVbmJ2wHuG3CUu19vRbVwdW0oChU+FdQNjzWET2DdmWbEbT1BYMuaZkcSEZEiungaCTfSSCPnQhU30rgldCP3j/em+9iWWJzqmZzWsehUVwVVu20gHarswMCJH2fsMTuOiIgU0cHlhxk9v6t9aBIDJ9LwJMz5MK/1j+LozhQ+je5KjwdbYXGymJzW8ajwqcAGdToJwMJf3UxOIiIil3Nix0leui6KDj2rFDiNxJxX43nixwhqNAkwIV3FocKnAht4b86YDEviWnAu4ZzJaUREpCAbP9rJsHorCWnuw/O/R3CK6hQ0jUTDHpVjnJ2rzSEKn+joaEaOHEl4eDienp7Uq1ePSZMmkZ6enq/tvHnzaNmyJR4eHgQGBjJu3DgTEpcPLW9tSF3no5zDi9/f3GZ2HBER+VtGagafP7SaLlW30WFYUz6ydSMddzpW2c7HD6xi1p0rcCYT4IJpJFT4lAaH6Ny8e/dusrOzmT17NvXr12f79u2MHj2as2fPMn36dHu7GTNm8PrrrzNt2jQ6duzI+fPnsdlsJiY3l8XJwqCmB3h7WzALvzrPDS+YnUhEpHI7vv0k7z28g3ejGhOb3QUAV9K5LWwDDz7nS8eRze1tBzxSeaaRKEsOO4DhtGnTePfdd+2FzenTp6lTpw6LFi2iV69eJV5vRRjA8EJL/rOJPs+0I9DpOMfSauDk4hAH+UREKpSNH+3krSkJfGHrQDruAAQ6Hef+Hru4779NdeVtKajwAxgmJSXh7+9vv79kyRKys7M5duwYTZo0ITg4mNtuu40jR46YmNJ8Pce1oCrJxGXXYtMnu82OIyJSaaSfSeezB/85nfW/C05nfTJmFYeSqjFpWYSKnjLmkIXPgQMHmDlzJvfff799mc1mIzs7m1deeYU333yTr7/+moSEBHr37l1gX6BcaWlpJCcn57lVJG7ebvQL3gHAwg9OmJxGRKTiO779JC9cG0WYbwJ3vt2FNWda4Eo6d4WvZN2HO1h7pjl3vtMVN29dcWsGUwufyZMnY7FYLnnbuHFjnufExMTQr18/Bg8ezKhRo+zLs7OzycjI4K233qJv37506tSJzz77jH379rFs2bJCM0ydOhVfX1/7LSQk5Kq9X7MMuiEbgIUb1TFORORq2TB/J3fXW0ndFj5MWhZBbHYggU7HmRIZxeG/EvmfrRvXjGhmdsxKz9Q+PvHx8cTHx1+yTVhYGB4eHkBO0RMZGUnHjh2ZN28eTk7/1G1z587lnnvu4ciRIwQHB9uX16pVi5deeonRo0cXuP60tDTS0tLs95OTkwkJCakwfXwAEg6cpmb9qmThQvTKo4R2Db78k0REpEAXTiVRs0l1vnlmI2/Nq8raMy3sbTp5b+Oh4SncMrW9juyUkaL28TH1qq6AgAACAoo2ENOxY8eIjIykXbt2zJ07N0/RA9C1a1cA9uzZYy98EhISiI+PJzQ0tND1uru74+7uXsJ34Bj861Wjq+8Wlie1ZtGbBxinwkdEpEQunkqiKikk88/VWUOs63lwoj8dhrW4zJrELA7RxycmJoaIiAhCQkKYPn06J0+eJC4ujri4OHubhg0bcuONN/Lwww+zevVqtm/fzrBhw2jcuDGRkZEmpi8fBnVPBGDhUm9zg4iIOKijG2L/Lnr+mUoiGV9qWk7wwrVRHNmWxEcHutFhWFOTk8qlOEThs3jxYvbv38/SpUsJDg4mKCjIfrvQRx99RMeOHRkwYAA9e/bE1dWVX375BVdXV5OSlx+DHgoDICqhBUmHk8wNIyLiQLIzs/lpygYG9ki0Fz0X+viVIzz/ewS1mtcwIZ0Ul8OO43O1VLRxfC7UxP0Au9Pr8cUjq7ntjS5mxxERKdfOJ57nk0c2MOPzIHam1f97qQH8MzGoM5lErz+pUZXLgQo/jo8U38BWOWMaLfo+2+QkIiLl16l9CbzcO4qw6smMmt+dnWn1qUoyj7eP4rX+f2gqCQfnEFNWSOkYNNyfaRvgx+hmZJ7PxMVDH7+ISK4DSw/xxiPRzN3WnlQiAAhxPsbD1+9j1Dtt8K2bs+yODZpKwpHpVNdFKvKprqz0LGp5JHLKqE7Um1vo+XBrsyOJiJhuzXvbeP2FM3x7rCPG3ydC2nju4vERCQyedg2uXuon6gh0qkvycXZz5gbrLgAW/i/R3DAiIibKSs/i2yfX0tVnK13ua8E3xzpj4ET/Ghv4ffpmNp1pzJ3vdFXRUwGp8KlkBt6cc0XCwq1hGNk62CcilUtqfCr/N+QPGlU5yi3TOrE6pSVupHFPgxVs/24/P57owLXj22Bxslx+ZeKQVPhUMn0ebY4baezPCGPPLwfNjiMiUiaObz/JxB5R1K15jrFf9ORAZijVLKd5tmsUh/5KYs7e7jS7sf7lVyQOT4VPJVO1dlWuDdgKwMK3D5ucRkTk6tr1wwFGN15OaIuqvLgiglNGdawuh5h56x8ciXPjpZWaHb2yUeFTCQ3slQrAwpX+JicRESkdRzfEsmzGZo5uiMXINoh6cws31FpP04H1+GBPD9LwoJP3Nr5+fA17zwYz7queVKlZxezYYgJdz1wJDXy0AWO/gDUpzTi5K54aTYo2X5qISHl08fxZdZ2PciirNQAWsrkpaD2PT6pCl/s0f5boiE+lFNKxNm08d5GNMz/N2G12HBGREito/qxDWSG4c44Hmi1nz+LDfBvTSUWP2KnwqaQGXXMcgIU/6aCfiDimhAOnef6O/QXOn/XFM1v5v+09aNA7rOyDSbmmwqeSGjiqFgC/xjTnfOJ5k9OIiBRd7JbjPNEhitD6Lsw7kH/kZGcyaXdzXROSiSNQ4VNJtb2zMbWdYjmLN1Fvbzc7jojIZdmiDnN/0+WEtfFj+sYIzlCVlh57uLfxcs2fJUWmwqeSsjhZGNR4HwALP081OY2ISOG2fbOXf4etokFkHWbv6kE67nStupUfJm1gy9mGzN7Vg+j1J1n2xhai159k5DzNnyWFUwePSmzQ7Z7MmgQLdzXgnWxDI5WKSLmy9oPtTH3+LAvjOgINAegXsIEJk9zoPq5VnrbBHYJ0lEeKREd8KrHIh1pQhTMcyw5i82e6uktEzGdkGyz5zyaurbaZzqObszCuIxayGRy8hk0f7+Lnkx3yFT0ixaHCpxLz8POgT1BO/55FHxw3OY2IVGbZmdl8++RaOlTdRZ9n2rEssQ0uZHBPgxXs+imaL490pu2/m5gdUyoAFT6V3KD+OR0CF66rZXISEamMMlIzmD96Jc2rHOSWaZ3YlNoUT1J5qNUfHFh9gjl7u9PoeqvZMaUCUeFTyQ0Y3xgL2fx5rglHN8SaHUdEKolzCed4e/Af1Pc5zvAPurErvR6+JOVMGrozlf9u6UndznXMjikVkAqfSq5GkwA6V/37dNeMfSanEZGKLulwElP7RhEWcIYHv+7J4axgajmd4NXrozh8xMJLKyM0jY5cVbqqSxjUNYHVv8Ci3z15wOwwIlJhHN0Qy74VcTToHoiblwtvjt7BO2vakEwEAGEuR3jiFhsj/u8aPP0jTM0qlYcKH2HQ2BCe/gV+P9mSM3Fn8A70NjuSiDi4iycOdSaTzL8Lnqbu+3l6eBxDZnTE1SvE3KBS6ehUl9C4v5X6rtGk487iGRrFWUSuTEETh2biRiuP3Sx4eh3bzlgZOqsbrl6uJieVykiFj2BxsjCwRTQAC7/NNDeMiDi0/b8fYtT1RwucOPSNl89x09SOOLnop0fMo2+fADBoqB8AP9oak5WeZW4YEXE4+5ZEM6zeShpfV4dfT3XI97gzmTToHmhCMpG8VPgIAF3vbUY1y2nijQDWztlhdhwRcRB7frYx1LqSxn1C+MjWjSxcuL7GBiZ0WaaJQ6VcUuEjALh6udI/dCcAC+clmJxGRMq7XT8c4N9hq2jaP5SPD3YjG2cG1FzPug938NOJDry8KlITh0q5pKu6xG7gIAufvAULt4TwqtlhRKRc2rlwPy8+eIIvDnfCoB4AgwLXMXGaD+3uuiZPW00cKuWRjviIXb/xzXAhg93p9di3JNrsOCJSjmxfsI/b666m+Y1WPj/cBQMnbgpay5+f7ub72I60u0vzaIljUOEjdr51fYnw3wrAorejzQ0jIuXC1q/3cmvwGlr8qwFfHskpeP5Vey2bP9/DgphOtLmjsdkRRYpFhY/kMSjyDAAL//AzN4iImGrLF3v4V+21tBrckG+OdQbg1jpr+OurvXxzrBOtb29kckKRklHhI3kMfCTnnP3KpOYkHDhtchoRKWt/frKLm4LW0WZIIxbEdsJCNreFrGbbt/v46mhnWt7a0OyIIldEhY/kEdYtmBYee8nChZ9f32l2HBEpI5s+3sWgwHW0u6sJ38d1xEI2Q+quZvv3Nr443IXmNzcwO6JIqdBVXZLPoHYxbFvVkIU/WPj3/5mdRkRK24WTh8buPM2UJ8/w44mcK7KcyGJI6FqeezuQJjd0MTmpSOlT4SP5DBpZg5dXwc9HmpN+Jh03bzezI4lIKblw8lAwgJzLzZ3I4s7wNTz3Tm0aXd/V1IwiV5NOdUk+7Yc2IdDpOCn4sPz/NGmpSEVx8eShYAEMbq2zhl2/HOZ/tm40ut5qZkSRq06Fj+Tj5OLEDQ32ALDwkxST04hIafjryz3cfu3JAiYPtTD2cU8a9g03JZdIWVPhIwUaNNgDgIU76mFkGyanEZGS2rlwP7eFrKH17Y1YfaYlOae3/uFMJvW71jInnIgJVPhIgXo93BwPznEoK5ht3+4zO46IFNP+3w8x1LqS5jda+epozjg8Q+qu5sVeUZo8VCo1dW6WAnkFeNG71joWHe/IotkxGrtDxEEcWnWUF+85yLy9nckiFICbg9YyZaY/LW7JuUpr+IZY9q86Tv2utQjuoMlDpXJR4SOFGtQ3nUUfwcI1ATxrdhgRuaRjG2N5ZcRe3t/emQyCAehfYwMvzPCm3V2d8rTV5KFSmTnEqa7o6GhGjhxJeHg4np6e1KtXj0mTJpGenp6n3YYNG+jVqxd+fn5Uq1aNPn36sGXLFnNCVwA3jM8Zkn792ebEbjluchoRKciJHSd5rF0U9TpU4/+29yQDN3pV+5PVs7fx44kOmjxU5CIOUfjs3r2b7OxsZs+ezY4dO3jjjTeYNWsWEyZMsLdJSUmhb9++1K1bl3Xr1rFy5Up8fHzo27cvGRkZJqZ3XIEta3JNlZzL2X+YscfkNCJyoVP7EnimcxThzb14488I0vCgm89fLHtjC78ltKXzvS3MjihSLlkMw3DIS3amTZvGu+++i81mA2Djxo106NCBw4cPExISAsC2bdto2bIl+/fvp169ekVab3JyMr6+viQlJeHj43PV8juKl3tH8dxvEQystY6FcR3NjiNS6SUdTuKNuzcz44+2pJCzj+pQZQcvPnOOPs+0w+JkMTmhiDmK+vvtEEd8CpKUlIS/v7/9fqNGjQgICGDOnDmkp6dz7tw55syZQ7NmzQgNDS10PWlpaSQnJ+e5yT8G3V8bgCXHW5Ian2pyGpHK60zcGab2jSI8LJspf0SQgg+tPPaw8Ln1rEtuSt9n26voESkChyx8Dhw4wMyZM7n//vvty6pWrUpUVBQff/wxnp6eeHt78+uvv/LTTz/h4lJ4H+6pU6fi6+trv+UeLZIczW9uQJjLEc7jyW9vbDM7jkilcy7hHDNujMJa+xwTFkdw2qhGE7cDfPXYGv5MacDAF69RwSNSDKYWPpMnT8ZisVzytnHjxjzPiYmJoV+/fgwePJhRo0bZl587d4577rmHrl27snbtWlatWkWzZs3o378/586dKzTDM888Q1JSkv125MiRq/Z+HZHFycLApjmnExd+nX6Z1iJSWtKS03h78B/Uq5HE+IURnDRqUN81mo8fWMW2lDBufb0zTi4O+W9XEVOZ2scnPj6e+Pj4S7YJCwvDwyNnFOGYmBgiIyPp2LEj8+bNw8npnz/6OXPmMGHCBGJjY+3L09PTqVatGnPmzGHIkCFFyqQ+Pvn99tqf9H6qLbWcThCTFqCdrchVkDtjelj76vw+5xAvfGLlSFYdAEKdjzJx6EHufrczLh4ahUSkIEX9/Tb1LyggIICAgIAitT127BiRkZG0a9eOuXPn5il6AFJTU3FycsJi+eeQb+797OzsUs1d2fQY0xyfp5I4nl2TDfO303Fkc7MjiVQo+WdMrwtAbadYnrttHyPf74Sbd7CpGUUqCof4p3tMTAwRERGEhIQwffp0Tp48SVxcHHFxcfY2vXv35vTp04wdO5Zdu3axY8cORowYgYuLC5GRkSamd3xu3m5cH7IDgIUfXvoInYgUz5F1MYye3zXfjOmTei5j/0k/HvisB27ebmZGFKlQHKLwWbx4Mfv372fp0qUEBwcTFBRkv+Vq3LgxixYtYuvWrXTu3Jnu3bsTExPDL7/8kqedlMzAG3LOiC7cVMfkJCIVg5FtsOQ/m7iu+3mMfLtiCxE3VcPT39OUbCIVmcOO43O1qI9PwRIOnKZm/apk4YLtjyOE99DVbyIltea9bUx4KpOoxDZ/LzHIOdKTw5lMotef1LQSIsVQ4cfxkbLlX68a3f1yLmdf9F+byWlEHNPWr/cyKHAdXe5rQVRiG9xI45E2fzBj0B+aMV2kjOjyACmyQT2TiPoeFi6rykNmhxFxIPt/P8SkkUf57FBnDJxwJpMRjVbz/Nx61O3cE4DBmjFdpEzoVNdFdKqrcPt/P0SD60JxIYP4Q6n41vU1O5JIuXZsYywvDN3HnN1dyPr735m3hazmhdmBNLreanI6kYpFp7qk1NXvFUoTtwNk4sovr+8wO45IuRW/5xSPt4+ifgc/3tvdgyxc6F9jA39+upsvDndR0SNiIhU+UiyD2uSMbL3wex0oFLlY8tFkpkRGYW3syuubIjiPJ919/mLFO1v58UQH2tzR2OyIIpWeCh8plkEjqgPw0+FmZKRmmJxGpHywz6dVN4PJUTkTiLbx3MXPL27kj9Mt6TampdkRReRvKnykWDqOaEqAJZ5Ew4+Vs7abHUfEVBmpGbx/93Ia1Ehk/MIIThnVaeRm48tH17AxuRH9ntOM6SLljQofKRZnN2duqLcLgEWfJJmcRsQc2ZnZfPbgapr6HePe//XgWHYQIc7HmDN8BduT6jJ4hiYQFSmv9JcpxTboX64ALNwajpGtvj5SeRjZBj9MXE+bqvu48+0u7M8Io4blJG/e/Af7EgK4Z253TSIqUs7pcvaL6HL2yzsTd4aAIBfS8GDH9/tpOqi+2ZFErprcWdPjj57jzQ+8WZ2S01/HhySevG4zD/+vPd6B3ianFBGHmJ1dHJN3oDfX1tjAzyc78MaTMUwKqqJRZqVCyjtreg5PUnmo43qe/KQV/vUizAsnIiWiU11SItW90wD4YE8PQq+pyZzhK0xOJFK6ls3YzKj53S6YNR0sZLPif9H8Z20E/vWqmZhOREpKhY8U29ENsXx6sLP9fjbO3De/M0c3xJqYSqR0RK88yogGK+g1viUXThwKYOBESny6OcFEpFSo8JFi27ciLs+/ggGycGH/quMmJRK5cnFbT/BQqz9o2L0m8/Z3x8CZnFnT/+FMJvW71jInoIiUChU+UmwNugfiRFaeZU5k6QdBHNLpg4lM6BJFvVZVmLm1Jxm4cZ3/JtZ9uIMPhq3UrOkiFYw6N0uxBXcI4r1hK7hvfmf7xIsuZHI2/pzJyUSK7uyJs7w1dAOvLWlNohEBQMcq23llSgbXjm8HwDUjoO9YzZouUpHocvaL6HL2oju6IZY9UTFMfNGF1SmtaO25m7Vx4bj7uJsdTaRQaclpvD9yLS9924Tj2TUBaO6+j5fHJzDwxWs00rKIgyrq77cKn4uo8Cm+mD/jaNXehXgjgIda/cF/t/Q0O5JIPlnpWfxvzBomzwvjUFYwAFaXQ7ww+ihD3uyEs5vzZdYgIuVZUX+/1cdHrljttoHMn3QQgLf+6smi59ebnEjkH0a2wTdPrKFF1WhGzOnGoaxggpziePeO5exOqs2//6+rih6RSkSFj5SK/pM68GjbKABGvFyPYxt1abuYy8g2WDx1Ex2q7uLW6Z3ZlV4Pf0sCr/WPYv9xH+7/tAeuXq5mxxSRMqbCR0rN1GWdaeu5i1NGde7qHUdWetblnyRyFayevY1I/7/oO6Edm1Kb4k0KE3tEYYt25okfI/AK8DI7ooiYRIWPlBp3H3c+W+BBFc4QldiGqf01mrOUrb++3MPAWuvpen8L/khqjTvnebRtFLadaUz5IwLfur5mRxQRk6nwkVLVsG84/zd6CwCTf+/Gqne3mhtIKoV9S6K5M2wVrW9vxA8nrsGZTEY1Ws6+tQnM2BRBjSYBZkcUkXJCV3VdRFd1lY6h1pV8fLAbdZ2PsmWfN9XC/cyOJBVI7ozp3v5ufPDqKebs7mIfU+r2kNW88H4QDfuGm5xSRMqSZmcXU/3fylasDYtmf0YYo7qt5esjHTU+ipSKvDOmG+TOpzWg5npeesuX1rd3MTWfiJRvOtUlV0XV2lX5/MNzuJLOtzGdmH2X+vvIldv61R5Gz+96wVxxFiDncvUfjl9D69sbmRlPRByACh+5atrd1YT/DFoNwKOfdWD7gn0mJxJHdWpfAhO6RNHxthCMfLstC/61PU3JJSKOR4WPXFWPfNOD62ts4Dye3H6HhdT4VLMjiQNJOHCa57pFEd7QhalrIjiPF5oxXUSuhAofuaqcXJyYtyyMQKfj7Eyrz6PdN5odSRzA6YOJTOwRRXh9J15eFUEKPrT23M13z6zj/btXaMZ0ESkxXdV1EV3VdXX8Pu1Pej/ZGgMnvnx0DYNndDY7kpRDiYeSeHPYZt74ow3J5Iy509JjD5MfOs1NU//pIH90w4UzpqvoERFNUlpiKnyungldopi6JgJfktiyIoWwbsFmR5JyIulwEv8dvpkZy9qQ9HfB09x9H5MfjOfmqR1xctHBaRG5NE1SKuXOlN+60sl7G0n4cuf1p8k8n2l2JDFZ8tFkXrouirAwg0nLIkjCl2bu+/jy0TX8daYet0zrrKJHREqV9ihSZly9XPn0Rz98SGLNmRZMvm6l2ZHEJCkxKbzcO4qwulk8/3sEiYYfTd3388Ujq9l6ph6DZ6jgEZGrQ3sWKVPhPUJ4/5EdALyyqgdLX99sciIpSykxKUztG0VYcAbP/RbBaaMajd0O8NmDq9maHM5tb3RRwSMiV5X2MFLmbnujC6MaLcfAibueDOLkrnizI8lVdibuDK9eH0V4cDoTFkeQYPjTyM3GJ2NWsT0ljCFvdcHZzfnyKxIRuUIqfMQU/13ZniZuB4jNDmREhA0jW33sK6KzJ84ybUAU4bXP8/QvEZwyqtPA9SD/u38VO1JCufOdrip4RKRMqfARU3gFePHFp1m4c54fT1zDf29ZbnYkKQVHN8SybMZm9i05yOsDo7AGpvLkTxHEGwHUd41m/uiV7EwO4a53VfCIiDl0OftFdDl72fq/IX8w9oueuJLO2o8P0PbfTcyOJCX0z+Shzlw4eajV5RDPDzvCXW93wsVD8yKLyNWhy9nFITzwaQ9uDlpLBm4MGeFBSkyK2ZGkBP76cneBk4dOG7CM3Um1Gf5BNxU9IlIuqPARU1mcLHywohEhzsfYlxHOuO5/mR1JimHXDwe4r8lyOtweXuDkoe2vq4arl6sp2URECuIwhc+gQYOoW7cuHh4eBAUFMXToUGJiYvK0OXz4MAMHDqRKlSoEBATw0EMPkZ6eblJiKSr/etX45M14nMjiI1s3Pn5gldmR5BKMbIPfXvuT/jU30HRgPd7b3YMM3NHkoSLiCBym8ImMjOTLL79kz549fPPNNxw4cIBbb73V/nhWVhYDBgzg7NmzrFy5ks8//5xvvvmG8ePHm5haiqr7uFZMilwBwAOzWrJvSbS5gSSftOQ05o1aSasq++j9VFt+PtkBC9ncHLSW5TP/0uShIuIQHLZz88KFC7nppptIS0vD1dWVn3/+mRtuuIEjR45Qu3ZtAD7//HOGDx/OiRMnitxRWZ2bzZOVnsW1NbexPKk17bx2svp4fdy83cyOVemd3BXPrHHbeSeqKcezawJQhTPc03ITD78RRr1rQ+1tNXmoiJilQnduTkhI4JNPPqFLly64uub0H1izZg3Nmze3Fz0Affv2JS0tjU2bNhW6rrS0NJKTk/PcxBzObs58sqQW/pYENqU2ZULEarMjVWq5/XfqNq3CxKURHM+uSbBzDK/1j+KILZO3/uqZp+gBCO4QRMQjrVX0iEi55VCFz1NPPUWVKlWoXr06hw8f5vvvv7c/FhcXR61aefsTVKtWDTc3N+Li4gpd59SpU/H19bXfQkJCrlp+ubzgDkF8+PQ+AF7fFMHPL2wwOVHlktt/Z8AF/XfO40l7r518Om41tuQaPPFjBNXC/cyOKiJSIqYWPpMnT8ZisVzytnHjRnv7J554gs2bN7N48WKcnZ25++67ufBMncViyfcahmEUuDzXM888Q1JSkv125MiR0n2TUmw3vtKRcS3+AGDY5DBitxw3OVHFd3H/nZ8u6r+zPqUJd8zsoiu0RMThmTqwxrhx4xgyZMgl24SFhdn/PyAggICAABo2bEiTJk0ICQlh7dq1dO7cmcDAQNatW5fnuadPnyYjIyPfkaALubu74+7ufkXvQ0rftOUdWRG0h7/ON+LuXpv49XgNTV55FeTtv9MNuLj/TieTE4qIlC5TC5/cQqYkco/0pKWlAdC5c2defvllYmNjCQrK6V+wePFi3N3dadeuXekEljLj4efB59+40m7AWX5LaMez3ZfRZ7AfDboHqv9IKdj1wwHefOIYH+3uwHkiAAh2juGhvnsZ9XZrqoX3NDegiMhV4hBXda1fv57169fTrVs3qlWrhs1mY+LEicTGxrJjxw7c3d3JysqidevW1KpVi2nTppGQkMDw4cO56aabmDlzZpFfS1d1lS8fjljByHndyZ0CwYks3hu2+u9lUhRHN8Syb0Uc9bvWYu/yOGZMy+Knkx3sj7fz2sn4exK59dUOOpUlIg6rqL/fDlH4bNu2jYcffpi//vqLs2fPEhQURL9+/XjuueeoU6eOvd3hw4cZM2YMS5cuxdPTkzvvvJPp06cX61SWCp/y5ci6GOp2CiJ33icAJ7KwrYwltGuwecEcRGHzZ1nI5sag9Tw2wZNuY1picSq8H5yIiCOoUIVPWVLhU74sm7GZa8e3ybc8yCmO8f13c89brXWFUQESDyXx0aObeXhBTy4sGsFgRP2VTJhVl/q9Qgt7uoiIw6nQ4/hI5dGgeyBOZF201CA2O5DHf4igjtWNe5ssZ+vXe03JV57EbjnOrDuX0zdgIzXCvHh4QQR5ix4AC3ePraqiR0QqLRU+Uq4FdwjivWGr80yF8H9DVvDe0OW08NjLObx4f3cPWg1uSE+/LXw9fg0ZqRkmpy47+38/xLQBUXSpuo06bWrwwGc9WHyqPZm40tDVhoXsPO01f5aIVHY61XURneoqnwqaCsHINljxzlZmvprKgmMdyPr7IsU6TrE8cO0eRr/ZjJrNapgZu9QZ2QZbvtjDgnfjWLC+DtvTGuR5vJP3Nm7ufoqbHw2jQe8w5gxfwX3zO5OFi33+LHUMF5GKSH18SkiFj2M6uiGWWY/u4b3VzThp5BQ7bqRxu3UD456rxjUjmpmcsOSy0rNYNXs7C+Ym8t02K9GZ/4wu7kIGkf5bubnPWW58oiG12wbme77mzxKRykCFTwmp8HFsaclpfPnkBt7+2I/1Z5vbl19TZTvj7krkttc64O5T/gesTEtO4/c3tvLtp+dZuK+xvZgD8CSVfkFb+deNWQx4spk6d4uIoMKnxFT4VBzr5+7g7ZdO84WtA+nkFDs1LCe5t8sO7n+jkelHP3LH18kdlDH5aDI/T9vOggXw05HmpPDP96+a5TQDw3dw822u9BnfAq8ALxOTi4iUPyp8SkiFT8VzYsdJ3n9kB+8ubcSx7Jxix5lMbq6zgQef8qL72LIfx+bC8XUsZNPcfT970kLtBRrk9FW6qdlebh5alR5jW2hwQRGRS1DhU0IqfCqujNQMvn9+IzPneLI8qbV9eUuPPYwbfJw7p7fj9KHkPEdhiiMtOY3Ew8kkHkkhMSaV07HnSTyRTmJ8JokJ2Zw+DYnJTsQkuPPDiWvIf6k5NHQ9yM1tD/Gv+2rQfmgTzU8mIlJEKnxKSIVP5bD16728/XwcH+9uzzlyTht5cpbzeGLghBNZPNJ2BZ0j3HMKl1NZfxcuFhJTXEhMdeP0OQ8S071IzKxCYraPfT0lNfeeFQx7v5tGURYRKQEVPiWkwqdyOX0wkQ8f2sJ/f2rAkew6l3/CZVjIxteSjJ9zCn6uZ/FzP0c1zzT8vDPwq5qFny9gGEyO6olxwTBazmQSvf6k6f2OREQcVVF/v02dnV3EbNXC/Ri/KIJWr22i91P5C5/m7nsI90vEr0oGft5Z+PkaVPMHP39n/Gq44lfLnWp1vPCrUwW/kKr41KmKk4sf4HfJ161TwPg6wR00vo6IyNWmIz4X0RGfyunohlhCr6n592SeOa72URiNryMiUno0V5dIMRQ0NUbOUZirV5AEdwgi4pHWKnpERMqQjvhcREd8KjcdhRERcUzq4yNSAsEdglTwiIhUYDrVJSIiIpWGCh8RERGpNFT4iIiISKWhwkdEREQqDRU+IiIiUmmo8BEREZFKQ4WPiIiIVBoqfERERKTSUOEjIiIilYYKHxEREak0VPiIiIhIpaG5ui6SO2drcnKyyUlERESkqHJ/ty8397oKn4ukpKQAEBISYnISERERKa6UlBR8fX0LfdxiXK40qmSys7OJiYmhatWqWCyWAtskJycTEhLCkSNH8PHxKeOE5Ye2wz+0LXJoO+TQdsih7fAPbYscV3M7GIZBSkoKtWvXxsmp8J48OuJzEScnJ4KDg4vU1sfHp1J/gXNpO/xD2yKHtkMObYcc2g7/0LbIcbW2w6WO9ORS52YRERGpNFT4iIiISKWhwqcE3N3dmTRpEu7u7mZHMZW2wz+0LXJoO+TQdsih7fAPbYsc5WE7qHOziIiIVBo64iMiIiKVhgofERERqTRU+IiIiEilocJHREREKg0VPsD//d//ER4ejoeHB+3atWPFihWFtv3222/p3bs3NWrUwMfHh86dO/Prr7/maTNv3jwsFku+2/nz56/2W7lixdkWUVFRBb7P3bt352n3zTff0LRpU9zd3WnatCkLFiy42m/jihVnOwwfPrzA7dCsWTN7G0f8TixfvpyBAwdSu3ZtLBYL33333WWf88cff9CuXTs8PDywWq3MmjUrXxtH+z4UdztU5H1EcbdFRd1HFHc7VNR9xNSpU+nQoQNVq1alZs2a3HTTTezZs+eyzzN7P1HpC58vvviCRx55hGeffZbNmzfTvXt3rr/+eg4fPlxg++XLl9O7d29++uknNm3aRGRkJAMHDmTz5s152vn4+BAbG5vn5uHhURZvqcSKuy1y7dmzJ8/7bNCggf2xNWvWcPvttzN06FD++usvhg4dym233ca6deuu9tspseJuh//+97953v+RI0fw9/dn8ODBedo52nfi7NmztGrVirfffrtI7Q8ePEj//v3p3r07mzdvZsKECTz00EN888039jaO+H0o7naoyPuI4m6LXBVtH1Hc7VBR9xF//PEHY8eOZe3atSxZsoTMzEz69OnD2bNnC31OudhPGJXcNddcY9x///15ljVu3Nh4+umni7yOpk2bGlOmTLHfnzt3ruHr61taEctMcbfFsmXLDMA4ffp0oeu87bbbjH79+uVZ1rdvX2PIkCFXnPdqudLvxIIFCwyLxWJER0fblznqdyIXYCxYsOCSbZ588kmjcePGeZbdd999RqdOnez3HfH7cKGibIeCVJR9xIWKsi0q6j7iQiX5TlTEfYRhGMaJEycMwPjjjz8KbVMe9hOV+ohPeno6mzZtok+fPnmW9+nTh9WrVxdpHdnZ2aSkpODv759n+ZkzZwgNDSU4OJgbbrgh37/2ypsr2RZt2rQhKCiIXr16sWzZsjyPrVmzJt86+/btW+TtW9ZK4zsxZ84crrvuOkJDQ/Msd7TvRHEV9llv3LiRjIyMS7Ypr9+H0lBR9hFXoiLtI0pDRd1HJCUlAeT7rl+oPOwnKnXhEx8fT1ZWFrVq1cqzvFatWsTFxRVpHa+//jpnz57ltttusy9r3Lgx8+bNY+HChXz22Wd4eHjQtWtX9u3bV6r5S1NJtkVQUBDvvfce33zzDd9++y2NGjWiV69eLF++3N4mLi7uirZvWbvS70RsbCw///wzo0aNyrPcEb8TxVXYZ52ZmUl8fPwl25TX70NpqCj7iJKoiPuIK1VR9xGGYfDYY4/RrVs3mjdvXmi78rCf0OzsgMViyXPfMIx8ywry2WefMXnyZL7//ntq1qxpX96pUyc6depkv9+1a1fatm3LzJkzeeutt0ov+FVQnG3RqFEjGjVqZL/fuXNnjhw5wvTp0+nRo0eJ1llelDTzvHnz8PPz46abbsqz3JG/E8VR0Ha7eLkjfh9KqiLuI4qjIu8jSqqi7iPGjRvH1q1bWbly5WXbmr2fqNRHfAICAnB2ds5XRZ44cSJftXmxL774gpEjR/Lll19y3XXXXbKtk5MTHTp0KNeV+5Vsiwt16tQpz/sMDAy84nWWpSvZDoZh8OGHHzJ06FDc3Nwu2dYRvhPFVdhn7eLiQvXq1S/Zprx+H65ERdtHlBZH30dciYq6j3jwwQdZuHAhy5YtIzg4+JJty8N+olIXPm5ubrRr144lS5bkWb5kyRK6dOlS6PM+++wzhg8fzqeffsqAAQMu+zqGYbBlyxaCgoKuOPPVUtJtcbHNmzfneZ+dO3fOt87FixcXa51l6Uq2wx9//MH+/fsZOXLkZV/HEb4TxVXYZ92+fXtcXV0v2aa8fh9KqiLuI0qLo+8jrkRF20cYhsG4ceP49ttvWbp0KeHh4Zd9TrnYT5RKF2kH9vnnnxuurq7GnDlzjJ07dxqPPPKIUaVKFXtv+6efftoYOnSovf2nn35quLi4GO+8844RGxtrvyUmJtrbTJ482fjll1+MAwcOGJs3bzZGjBhhuLi4GOvWrSvz91ccxd0Wb7zxhrFgwQJj7969xvbt242nn37aAIxvvvnG3mbVqlWGs7Oz8Z///MfYtWuX8Z///MdwcXEx1q5dW+bvr6iKux1y3XXXXUbHjh0LXKcjfidSUlKMzZs3G5s3bzYAY8aMGcbmzZuNQ4cOGYaRfzvYbDbDy8vLePTRR42dO3cac+bMMVxdXY2vv/7a3sYRvw/F3Q4VeR9R3G1RUfcRxd0OuSraPuKBBx4wfH19jaioqDzf9dTUVHub8rifqPSFj2EYxjvvvGOEhoYabm5uRtu2bfNcijds2DCjZ8+e9vs9e/Y0gHy3YcOG2ds88sgjRt26dQ03NzejRo0aRp8+fYzVq1eX4TsqueJsi1dffdWoV6+e4eHhYVSrVs3o1q2b8eOPP+Zb51dffWU0atTIcHV1NRo3bpxnp1deFWc7GIZhJCYmGp6ensZ7771X4Poc8TuReylyYd/1grZDVFSU0aZNG8PNzc0ICwsz3n333XzrdbTvQ3G3Q0XeRxR3W1TUfURJ/jYq4j6ioG0AGHPnzrW3KY/7Ccvf4UVEREQqvErdx0dEREQqFxU+IiIiUmmo8BEREZFKQ4WPiIiIVBoqfERERKTSUOEjIiIilYYKHxEREak0VPiISIUUERHBI488AkBYWBhvvvmmqXlEpHxQ4SMiFd6GDRu49957i9RWRZJIxeZidgARkautRo0aZkcQkXJCR3xExOGdPXuWu+++G29vb4KCgnj99dfzPH7xUZzJkydTt25d3N3dqV27Ng899BCQc3rs0KFDPProo1gsFiwWS1m+DREpAyp8RMThPfHEEyxbtowFCxawePFioqKi2LRpU4Ftv/76a9544w1mz57Nvn37+O6772jRogUA3377LcHBwbzwwgvExsYSGxtblm9DRMqATnWJiEM7c+YMc+bM4aOPPqJ3794AzJ8/n+Dg4ALbHz58mMDAQK677jpcXV2pW7cu11xzDQD+/v44OztTtWpVAgMDy+w9iEjZ0REfEXFoBw4cID09nc6dO9uX+fv706hRowLbDx48mHPnzmG1Whk9ejQLFiwgMzOzrOKKiMlU+IiIQzMMo1jtQ0JC2LNnD++88w6enp6MGTOGHj16kJGRcZUSikh5osJHRBxa/fr1cXV1Ze3atfZlp0+fZu/evYU+x9PTk0GDBvHWW28RFRXFmjVr2LZtGwBubm5kZWVd9dwiYg718RERh+bt7c3IkSN54oknqF69OrVq1eLZZ5/Fyangf9fNmzePrKwsOnbsiJeXF//73//w9PQkNDQUyLkCbPny5QwZMgR3d3cCAgLK8u2IyFWmwkdEHN60adM4c+YMgwYNomrVqowfP56kpKQC2/r5+fGf//yHxx57jKysLFq0aMGiRYuoXr06AC+88AL33Xcf9erVIy0trdin0kSkfLMY+qsWERGRSkJ9fERERKTSUOEjIiIilYYKHxEREak0VPiIiIhIpaHCR0RERCoNFT4iIiJSaajwERERkUpDhY+IiIhUGip8REREpNJQ4SMiIiKVhgofERERqTRU+IiIiEil8f8CkjMlRA3yxgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "#the x axis is the atomic number\n",
    "\n",
    "ax.plot(distances,listEne,'.r-',label=\"UHF\")\n",
    "ax.plot(distances,listEneRHF,'.b-',label=\"RHF\")\n",
    "ax.set_xlabel(\"dist\")\n",
    "ax.set_ylabel(\"E(eV)\")\n",
    "ax.set_title(\"Dissociation curve for H2: Comparing UHF and RHF\")\n",
    "ax.legend();\n",
    "#if you see the two curves showing the same energy at dissoviactions it means that the UHF method has not converged into the lower unpaired singlet. \n",
    "# Need to perturb the guess density matrix in that case"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6f104f7",
   "metadata": {},
   "source": [
    "# Structure optimization\n",
    "\n",
    "A faster way to find the minimum of energy in the potential energy surface is by using an optimization algorithm. These algorithms will use the first derivative of the energy (the gradient) to move towards the minimum. It is an iterative process and it requires a threshold of gradient and energy for the computer to know when the minimum is found."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b77d977",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "e938a666",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimizer: Optimization complete!\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "-1.1313348258320226"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "psi4.core.clean()\n",
    "h2_geo = psi4.geometry(\"\"\"\n",
    "H 0.0 0.0 0.0 \n",
    "H 0.0 0.0 0.6\n",
    "\n",
    "units angstrom\n",
    "\"\"\")\n",
    "\n",
    "psi4.core.set_output_file('h2_opt.txt', False)\n",
    "psi4.set_options({'reference': 'rhf'})\n",
    "psi4.optimize('scf/6-31G**', molecule=h2_geo)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6788a1a",
   "metadata": {},
   "source": [
    "Your turn\n",
    "* Open the h2_opt.txt file and see the bond distance at which the algorithm has optimized the structure. Compare it with out previous scan.\n",
    "* What is the criteria for gradient and energy to decide when convergence is achieved?\n",
    "* Look at the code, what is the first loop doing? And the second?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "bc833ea5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    Total Energy =                         -1.1139361721990531\n",
      "\n",
      "    Total Energy =                         -1.1312540898766401\n",
      "\n",
      "    Total Energy =                         -1.1313243629427232\n",
      "\n",
      "    Total Energy =                         -1.1313334462445632\n",
      "\n",
      "    Total Energy =                         -1.1313348258320226\n",
      "\n",
      "**Optimization step\n",
      "   Step     Total Energy     Delta E     MAX Force     RMS Force      MAX Disp      RMS Disp    ~\n",
      "\n",
      "  --------------------------------------------------------------------------------------------- ~\n",
      "\n",
      "    Convergence Criteria    1.00e-06 *    3.00e-04 *             o    1.20e-03 *             o  ~\n",
      "\n",
      "  --------------------------------------------------------------------------------------------- ~\n",
      "\n",
      "      1      -1.11393617   -1.11e+00      1.62e-01      1.62e-01 o    2.31e-01      2.31e-01 o  ~\n",
      "\n",
      "**Optimization step\n",
      "   Step     Total Energy     Delta E     MAX Force     RMS Force      MAX Disp      RMS Disp   \n",
      "\n",
      "  ---------------------------------------------------------------------------------------------\n",
      "\n",
      "    Convergence Criteria    1.00e-06 *    3.00e-04 *             o    1.20e-03 *             o\n",
      "\n",
      "  ---------------------------------------------------------------------------------------------\n",
      "\n",
      "      2      -1.13125409   -1.73e-02      8.31e-03      8.31e-03 o    1.25e-02      1.25e-02 o  ~\n",
      "\n",
      "**Optimization step\n",
      "   Step     Total Energy     Delta E     MAX Force     RMS Force      MAX Disp      RMS Disp   \n",
      "\n",
      "  ---------------------------------------------------------------------------------------------\n",
      "\n",
      "    Convergence Criteria    1.00e-06 *    3.00e-04 *             o    1.20e-03 *             o\n",
      "\n",
      "  ---------------------------------------------------------------------------------------------\n",
      "\n",
      "      3      -1.13132436   -7.03e-05      2.95e-03      2.95e-03 o    4.53e-03      4.53e-03 o  ~\n",
      "\n",
      "**Optimization step\n",
      "   Step     Total Energy     Delta E     MAX Force     RMS Force      MAX Disp      RMS Disp   \n",
      "\n",
      "  ---------------------------------------------------------------------------------------------\n",
      "\n",
      "    Convergence Criteria    1.00e-06 *    3.00e-04 *             o    1.20e-03 *             o\n",
      "\n",
      "  ---------------------------------------------------------------------------------------------\n",
      "\n",
      "      4      -1.13133345   -9.08e-06      1.07e-03      1.07e-03 o    2.51e-03      2.51e-03 o  ~\n",
      "\n",
      "**Optimization step\n",
      "   Step     Total Energy     Delta E     MAX Force     RMS Force      MAX Disp      RMS Disp   \n",
      "\n",
      "  ---------------------------------------------------------------------------------------------\n",
      "\n",
      "    Convergence Criteria    1.00e-06 *    3.00e-04 *             o    1.20e-03 *             o\n",
      "\n",
      "  ---------------------------------------------------------------------------------------------\n",
      "\n",
      "      5      -1.13133483   -1.38e-06      3.45e-05 *    3.45e-05 o    8.33e-05 *    8.33e-05 o  ~\n",
      "\n"
     ]
    }
   ],
   "source": [
    "myOutput = open('h2_opt.txt','r').readlines()\n",
    "for line in myOutput:\n",
    "    if 'Total Energy =' in line:\n",
    "        print(line)\n",
    "for i in range(len(myOutput)):\n",
    "    if 'Step     Total Energy' in myOutput[i]:\n",
    "        print(\"**Optimization step\")\n",
    "        print(myOutput[i])\n",
    "        print(myOutput[i+1])\n",
    "        print(myOutput[i+2])\n",
    "        print(myOutput[i+3])\n",
    "        print(myOutput[i+4])\n",
    "        i = i + 4\n",
    "    i = i + 1\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "928b39eb",
   "metadata": {},
   "source": [
    "## Optimization of the molecule of water"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "648f6268",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimizer: Optimization complete!\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "-76.01072051843589"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "psi4.core.clean()\n",
    "psi4.core.set_output_file('h2o_opt.txt', False)\n",
    "h2o = psi4.geometry(\"\"\"\n",
    "0 1\n",
    "O\n",
    "H 1 0.96\n",
    "H 1 0.96 2 104.5\n",
    "\"\"\")\n",
    "\n",
    "psi4.set_options({'reference': 'rhf'})\n",
    "psi4.optimize('scf/6-31g*',molecule=h2o)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d616300",
   "metadata": {},
   "source": [
    "* See below how we can use psi4 to look at atomic distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "e7eb21e3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Distances between atoms in Bohr\n",
      "[[0.         1.79021278 1.79021278]\n",
      " [1.79021278 0.         2.85136018]\n",
      " [1.79021278 2.85136018 0.        ]]\n",
      "[[0.         0.9473398  0.9473398 ]\n",
      " [0.9473398  0.         1.50887483]\n",
      " [0.9473398  1.50887483 0.        ]]\n"
     ]
    }
   ],
   "source": [
    "#we can get the distances from the optimized structure directly\n",
    "myDistMat = psi4.core.Molecule.distance_matrix(h2o)\n",
    "print(myDistMat.name)\n",
    "print(myDistMat.np)\n",
    "\n",
    "print(myDistMat.np*psi4.constants.bohr2angstroms)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e108e4cf",
   "metadata": {},
   "source": [
    "# Questions\n",
    "\n",
    "* What is more stable for the molecule of F2, the triplet or the singlet?\n",
    "* What is more stable spin state for the molecule of NO? What is the Mulliken charge on these two atoms?\n",
    "* Calculate the bonding curve for H2, H2(+), H2(-) and plot them in the same graph. \n",
    "* Optimize the structure of HF and H2O and report the bond distances. (you will have to look at the output file or use psi4.core.Molecule.distance_matrix)\n",
    "* Optimize the structure of OH(-) and single point energy for F(-). Use these results to predict what is a better acid, H2O or HF?\n",
    "* Extra: you may also include in your acid comparison other molecules such as NH3, H2S, HCl, CH4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84ea3c82",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pchem",
   "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.12.2"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}