Car-Parrinello.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math as m\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References\n",
    "\n",
    "The referenced chapters, formulae and problems are from the book [Computational Physics, by Jos Thijssen](https://www.cambridge.org/core/books/computational-physics/BEE73B0139D4A9993193B57CDC62096E#fndtn-information).\n",
    "I also added some info on the blog: https://compphys.go.ro/car-parrinello-quantum-molecular-dynamics/\n",
    "\n",
    "What follows is copied directly from hatree-fock.ipynb, check that one for details. The interesting new part is after this portion. If you already looked over the Hartree-Fock, you can skip it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3.26, 4.16\n",
    "def Gaussian(alpha, r):\n",
    "    return m.exp(-alpha*r*r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha=(13.00773, 1.962079, 0.444529, 0.1219492)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see 4.114 and 4.116\n",
    "def F0(t):\n",
    "    if t==0:\n",
    "        return 1.\n",
    "    p = m.sqrt(t)\n",
    "    a = 1. / p\n",
    "    return a * m.sqrt(m.pi) / 2. * m.erf(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 4.98\n",
    "def Rp(alpha, beta, Ra, Rb):\n",
    "    return (alpha*Ra + beta*Rb) / (alpha + beta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 4.100\n",
    "def OverlapTwoCenters(alpha, beta, Ra, Rb):\n",
    "    difR = Ra - Rb\n",
    "    len2 = difR.dot(difR)\n",
    "    aplusb = alpha + beta\n",
    "    ab = alpha * beta / aplusb\n",
    "    return m.pow(m.pi / aplusb, 3./2.) * m.exp(-ab * len2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 4.103\n",
    "def KineticTwoCenters(alpha, beta, Ra, Rb):\n",
    "    difR = Ra - Rb\n",
    "    len2 = difR.dot(difR)\n",
    "    aplusb = alpha + beta\n",
    "    ab = alpha * beta / aplusb\n",
    "    O = m.pow(m.pi/aplusb, 3./2.) * m.exp(-ab * len2) # it's actually the overlap, check the OverlapTwoCenters\n",
    "    return ab * (3. - 2. * ab * len2) * O #this can be optimized with already computed overlap, see above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 4.115\n",
    "def Nuclear(alpha, beta, Ra, Rb, Rc, Z = 1.):\n",
    "    aplusb = alpha + beta\n",
    "    ab = alpha * beta / aplusb\n",
    "    difR = Ra - Rb\n",
    "    len2 = difR.dot(difR)\n",
    "    difRc = Rp(alpha, beta, Ra, Rb) - Rc\n",
    "    len2c = difRc.dot(difRc)\n",
    "    K = m.exp(-ab*len2)\n",
    "    return -2. * m.pi * Z / aplusb * K * F0(aplusb*len2c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 4.123\n",
    "def TwoElectronTwoCenter(alpha, beta, gamma, delta, Ra, Rb, Rc, Rd):\n",
    "    RP = Rp(alpha, gamma, Ra, Rc)\n",
    "    RQ = Rp(beta, delta, Rb, Rd)\n",
    "    alphaplusgamma = alpha + gamma\n",
    "    betaplusdelta = beta + delta\n",
    "    Rac = Ra - Rc\n",
    "    Rbd = Rb - Rd\n",
    "    Rpq = RP - RQ\n",
    "    Racl2 = Rac.dot(Rac)\n",
    "    Rbdl2 = Rbd.dot(Rbd)\n",
    "    Rpql2 = Rpq.dot(Rpq)\n",
    "    return 2. * m.pow(m.pi, 5./2.) / (alphaplusgamma * betaplusdelta * m.sqrt(alphaplusgamma+betaplusdelta)) * \\\n",
    "           m.exp(-alpha*gamma/alphaplusgamma*Racl2 - beta*delta/betaplusdelta*Rbdl2) *\\\n",
    "           F0(alphaplusgamma*betaplusdelta / (alphaplusgamma+betaplusdelta) * Rpql2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "basisSize = 4 # for each atom\n",
    "\n",
    "X = 1. # distance between atoms\n",
    "\n",
    "R0 = np.array([0, 0, 0])\n",
    "R1 = np.array([X, 0, 0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "H = np.zeros((basisSize * 2, basisSize * 2))\n",
    "O = np.zeros((basisSize * 2, basisSize * 2))\n",
    "\n",
    "for i in range(basisSize):\n",
    "    a = alpha[i]\n",
    "    for j in range(basisSize):\n",
    "        b = alpha[j]\n",
    "        \n",
    "        O[i, j] = OverlapTwoCenters(a, b, R0, R0)\n",
    "        O[i, 4 + j] = OverlapTwoCenters(a, b, R0, R1)\n",
    "        O[4 + i, j] = O[i, 4 + j]\n",
    "        O[4 + i, 4 + j] = O[i, j]   \n",
    "        \n",
    "        H[i, j] = KineticTwoCenters(a, b, R0, R0) + Nuclear(a, b, R0, R0, R0) + Nuclear(a, b, R0, R0, R1)\n",
    "        H[i, 4 + j] = KineticTwoCenters(a, b, R0, R1) + Nuclear(a, b, R0, R1, R0) + Nuclear(a, b, R0, R1, R1)\n",
    "        H[4 + i, j] =  H[i, 4 + j]\n",
    "        H[4 + i, 4 + j] = H[i, j]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "Q = np.zeros((basisSize*2, basisSize*2, basisSize*2, basisSize*2))\n",
    "\n",
    "for i in range(basisSize):\n",
    "    a = alpha[i]\n",
    "    basisSizei = basisSize + i\n",
    "    for j in range(basisSize):\n",
    "        b = alpha[j]\n",
    "        basisSizej = basisSize + j\n",
    "        for k in range(basisSize):\n",
    "            c = alpha[k]\n",
    "            basisSizek = basisSize + k\n",
    "            for l in range(basisSize):\n",
    "                basisSizel = basisSize + l\n",
    "                d = alpha[l]        \n",
    "                \n",
    "                Q[i, j, k, l]=TwoElectronTwoCenter(a, b, c, d, R0, R0, R0, R0)\n",
    "                Q[i, j, k, basisSizel]=TwoElectronTwoCenter(a, b, c, d, R0, R0, R0, R1)\n",
    "                Q[i, j, basisSizek, l]=TwoElectronTwoCenter(a, b, c, d, R0, R0, R1, R0)\n",
    "                Q[i, j, basisSizek, basisSizel]=TwoElectronTwoCenter(a, b, c, d, R0, R0, R1, R1)\n",
    "                Q[i, basisSizej, k, l]=TwoElectronTwoCenter(a, b, c, d, R0, R1, R0, R0)\n",
    "                Q[i, basisSizej, k, basisSizel]=TwoElectronTwoCenter(a, b, c, d, R0, R1, R0, R1)\n",
    "                Q[i, basisSizej, basisSizek, l]=TwoElectronTwoCenter(a, b, c, d, R0, R1, R1, R0)\n",
    "                Q[i, basisSizej, basisSizek, basisSizel]=TwoElectronTwoCenter(a, b, c, d, R0, R1, R1, R1)\n",
    "                Q[basisSizei, j, k, l]=TwoElectronTwoCenter(a, b, c, d, R1, R0, R0, R0)\n",
    "                Q[basisSizei, j, k, basisSizel]=TwoElectronTwoCenter(a, b, c, d, R1, R0, R0, R1)\n",
    "                Q[basisSizei, j, basisSizek, l]=TwoElectronTwoCenter(a, b, c, d, R1, R0, R1, R0)\n",
    "                Q[basisSizei, j, basisSizek, basisSizel]=TwoElectronTwoCenter(a, b, c, d, R1, R0, R1, R1)\n",
    "                Q[basisSizei, basisSizej, k, l]=TwoElectronTwoCenter(a, b, c, d, R1, R1, R0, R0)\n",
    "                Q[basisSizei, basisSizej, k, basisSizel]=TwoElectronTwoCenter(a, b, c, d, R1, R1, R0, R1)\n",
    "                Q[basisSizei, basisSizej, basisSizek, l]=TwoElectronTwoCenter(a, b, c, d, R1, R1, R1, R0)\n",
    "                Q[basisSizei, basisSizej, basisSizek, basisSizel]=TwoElectronTwoCenter(a, b, c, d, R1, R1, R1, R1)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "Qt = np.zeros((basisSize*2, basisSize*2, basisSize*2, basisSize*2))\n",
    "for p in range(2*basisSize):\n",
    "    for q in range(2*basisSize):\n",
    "        for r in range(2*basisSize):\n",
    "            for s in range(2*basisSize):\n",
    "                Qt[p, q, r, s] = 2. * Q[p, q, r, s] - Q[p, r, s, q]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Car-Parrinello\n",
    "\n",
    "From here is the new part, specific to Car-Parrinello method.\n",
    "\n",
    "#### First, the part from 9.3.1, Car-Parrinello only for the electronic part.\n",
    "\n",
    "The computation can take a while, maybe it's not a good idea to launch it in binder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "v = 1. / m.sqrt(O.sum())\n",
    "\n",
    "C = np.array([v, v, v, v, v, v, v, v])\n",
    "Cprev = C.copy()\n",
    "\n",
    "dt = 0.1 # time step\n",
    "dt2 = dt * dt\n",
    "dt4 = dt2 * dt2\n",
    "gamma = 1. # frictional constant\n",
    "\n",
    "mass = 1\n",
    "massMinusGamma = mass - 0.5*gamma*dt\n",
    "massPlusGamma = mass + 0.5*gamma*dt\n",
    "\n",
    "numPoints = 50\n",
    "energies = np.zeros((numPoints))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "F = np.zeros((2*basisSize, 2*basisSize))\n",
    "oldE = 100\n",
    "\n",
    "for cycle in range(numPoints):\n",
    "\n",
    "    # Fock matrix computation\n",
    "    \n",
    "    #for i in range(2*basisSize):\n",
    "    #    for j in range(2*basisSize):\n",
    "    #        F[i, j] = H[i, j]\n",
    "    #        for k in range(2*basisSize):\n",
    "    #            for l in range(2*basisSize):\n",
    "    #                F[i, j] += Qt[i, k, j, l] * C[k] * C[l]\n",
    "    \n",
    "    F = H + np.einsum('ikjl,k,l', Qt, C, C)\n",
    "        \n",
    "    # compute energy\n",
    "    Eg = C.dot(H + F).dot(C) + 1. / X\n",
    "\n",
    "    #print(Eg)\n",
    "    \n",
    "    energies[cycle] = Eg\n",
    "    if abs(oldE-Eg) < 1E-12:\n",
    "        break    \n",
    "        \n",
    "    # Verlet\n",
    "                    \n",
    "    # compute Ct - 9.31, but with friction force added\n",
    "    \n",
    "    Ct = (2. * mass * C - massMinusGamma * Cprev - 4. * F.dot(C) * dt2) / massPlusGamma\n",
    "                \n",
    "    # determine lambda - see 9.32 - but be careful, it's wrong! Get it from 9.28 by replacing C[r] = Ct[r] - lambda * h^2 * sum(S[r, s]*C[s]), h^4 and h^2 are missing (here h is dt)\n",
    "     \n",
    "    OC = O.dot(C)\n",
    "    OCt = O.dot(Ct)\n",
    "    OOC = O.dot(OC)\n",
    "    \n",
    "    a = OOC.dot(OC) * dt4\n",
    "    b = -2. * OC.dot(OCt) * dt2\n",
    "    c = OCt.dot(Ct) - 1.\n",
    "    \n",
    "    delta = b*b - 4.*a*c\n",
    "    if delta < 0:\n",
    "        print(\"Delta negative!\")\n",
    "        break\n",
    "        \n",
    "    sdelta = m.sqrt(delta)\n",
    "    lam1 = (-b-sdelta) / (2. * a)\n",
    "    lam2 = (-b+sdelta) / (2. * a)\n",
    "\n",
    "    if lam1 < 0:\n",
    "        lam = lam2\n",
    "    else:\n",
    "        lam = lam1\n",
    "    \n",
    "    # now adjust the newly computed Cs    \n",
    "    \n",
    "    Ct -= lam * dt2 * OC\n",
    "        \n",
    "    # switch to the next step\n",
    "    Cprev = C\n",
    "    C = Ct            \n",
    "    oldE = Eg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1.0781461092241256\n"
     ]
    }
   ],
   "source": [
    "print(Eg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD7CAYAAABnoJM0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAsH0lEQVR4nO3deXzc9X3n8ddnJI2OkXVbsg7LN/gEY4RxOcLphRCKA1lSyFHSJLBpQ3Nstwlttm22j0fTtE3a3SyExiQk3qSBnARInQRjE4cjGGzA2MbGlm35kGSd1n1rvvvHzBhZ6BrNjGTr934+HnrMzG9+x/fnBL31+57mnENERLzLN90FEBGR6aUgEBHxOAWBiIjHKQhERDxOQSAi4nEKAhERj4spCMwsz8y2mNmh8GvuKPt93sz2mdleM3vMzNKiOV5ERBIn1ieCB4CtzrklwNbw57OYWSnwGaDCObcSSALumujxIiKSWMkxHr8BuDb8fhPwW+CLo1wn3cz6gQygJsrjz1JQUODmz58/uRKLiHjUrl27Gp1zs4dvjzUIipxztQDOuVozKxy+g3Ou2sy+BhwHuoFnnHPPTPT4CDO7D7gPoLy8nJ07d8ZYdBERbzGzYyNtH7dqyMyeDdftD//ZMMEL5xL6y38BUAIEzOwj0RQewDm30TlX4ZyrmD37XYEmIiKTNO4TgXPuxtG+M7M6MysO/zVfDNSPsNuNwFHnXEP4mJ8DVwA/ACZyvIiIJFCsjcVPAfeE398DPDnCPseBdWaWYWYG3ADsj+J4ERFJoFiD4KvAejM7BKwPf8bMSsxsM4BzbgfwU+A1YE/4mhvHOl5ERKaOnY/TUFdUVDg1FouIRMfMdjnnKoZv18hiERGPUxCIiHicp4Jg6/46vvnbyukuhojIOcVTQfBCZSMPbqvkfGwXERFJFE8FQUl2Ol19g7R1D0x3UUREzhneCoKcdABqWrunuSQiIucOTwVBcU4aADUtCgIRkQhPBUFJduSJoGeaSyIicu7wVBDMnpVKss+o1ROBiMgZngqCJJ9RlJWmqiERkSE8FQQAJTlpqhoSERnCg0GQTq16DYmInOG5ICjOTudUaw/BoAaViYiAB4OgJCeN/kFHY0fvdBdFROSc4L0gUBdSEZGzeC4IIoPK1IVURCTEc0EQeSKoVhCIiAAeDIKcjBTSU5KoVdWQiAjgwSAwM4pz0tSFVEQkzHNBAKHqoeoWPRGIiIBHg6A4O02NxSIiYZ4MgpKcdBo6eukbCE53UUREpp1HgyAN56CuTdVDIiKeDILiyKAyVQ+JiHgzCCJLVqoLqYiIZ4MgvGSlupCKiHgzCDL8yWSnp6hqSEQEjwYBhNcl0FgCEZHYgsDM8sxsi5kdCr/mjrLf581sn5ntNbPHzCwtvP3LZlZtZm+Ef26JpTzRKMnWSmUiIhD7E8EDwFbn3BJga/jzWcysFPgMUOGcWwkkAXcN2eXfnHOrwz+bYyzPhBXnaO1iERGIPQg2AJvC7zcB7x9lv2Qg3cySgQygJsbrxqwkJ53W7n66+gamuygiItMq1iAocs7VAoRfC4fv4JyrBr4GHAdqgVbn3DNDdrnfzN40s0dHq1pKhDML1KidQEQ8btwgMLNnw3X7w382TOQC4V/uG4AFQAkQMLOPhL9+GFgErCYUEl8f4zz3mdlOM9vZ0NAwkUuPqTg73IVU1UMi4nHJ4+3gnLtxtO/MrM7Mip1ztWZWDNSPsNuNwFHnXEP4mJ8DVwA/cM7VDTnXI8AvxyjHRmAjQEVFRcwrz78zqExBICLeFmvV0FPAPeH39wBPjrDPcWCdmWWYmQE3APsBwuERcTuwN8byTNic7DTMVDUkIhJrEHwVWG9mh4D14c+YWYmZbQZwzu0Afgq8BuwJX3Nj+Ph/NrM9ZvYmcB3w+RjLM2EpST5mZ6bqiUBEPG/cqqGxOOeaCP2FP3x7DXDLkM9/B/zdCPt9NJbrx6okJ11PBCLieZ4dWQyhOYc035CIeJ2ng6A4OzTNhHMxtz2LiJy3PB4EaXT3D9LS1T/dRRERmTaeDoLScBdSVQ+JiJd5OgiKI2MJ1GAsIh7m6SAoydYCNSIing6CgsxUUpJMXUhFxNM8HQQ+nzEnO02DykTE0zwdBBDqQqqJ50TEyzwfBKUaXSwiHuf5ICjOTqOurYfBoAaViYg3KQhy0hkIOho7eqe7KCIi08LzQVCaE+pCWq12AhHxKM8HQXG2BpWJiLd5PggiaxerC6mIeJXngyArPZmAP0lVQyLiWZ4PAjOjOCddVUMi4lmeDwIIdSFV1ZCIeJWCgNCgsmo9EYiIRykICPUcauzopXdgcLqLIiIy5RQEQHF4LEF9mwaViYj3KAiAgkw/gEYXi4gnKQiAvEAqAM2dfdNcEhGRqacgAPIDoSeCJgWBiHiQggDID1cNNXUoCETEexQEQIY/mbQUH82daiMQEe9REITlB1JVNSQinqQgCMsL+NVYLCKeFFMQmFmemW0xs0Ph19xR9vusme01s31m9rloj58KCgIR8apYnwgeALY655YAW8Ofz2JmK4F7gbXAxcCtZrZkosdPlfxMvxqLRcSTYg2CDcCm8PtNwPtH2GcZ8LJzrss5NwBsB26P4vgpkR/w06TGYhHxoFiDoMg5VwsQfi0cYZ+9wHvMLN/MMoBbgLlRHD8l8gKp9PQH6eobmK4iiIhMi+TxdjCzZ4E5I3z1pYlcwDm338z+CdgCdAC7gah/25rZfcB9AOXl5dEePq4zg8o6+sjIG/efRURkxhj3N55z7sbRvjOzOjMrds7VmlkxUD/KOb4DfCd8zFeAk+GvJnR8+BwbgY0AFRUVbrxyRysyqKy5s4+5eRnxPr2IyDkr1qqhp4B7wu/vAZ4caSczKwy/lgN3AI9Fc/xUyDszzYTaCUTEW2INgq8C683sELA+/BkzKzGzzUP2+5mZvQU8DXzaOXd6rOOnQ3544jn1HBIRr4mpMtw51wTcMML2GkKNwpHPV0dz/HTIG1I1JCLiJRpZHBbwJ+FP9ikIRMRzFARhZkZBwE+jqoZExGMUBEPkZfo1A6mIeI6CYIi8QKqqhkTEcxQEQ4SmmVAQiIi3KAiG0AykIuJFCoIh8jP9dPUN0t03ON1FERGZMgqCIfI1ulhEPEhBMEReeHSxqodExEsUBEO8M9+QgkBEvENBMERB5jtTUYuIeIWCYIjIE4EGlYmIlygIhshMTcaf5FPVkIh4ioJgCDMLjSVQ1ZCIeIiCYBgNKhMRr1EQDJOf6adRQSAiHqIgGCY/oBlIRcRbFATD5AVS1UYgIp6iIBgmP9NPZ98gPf2ab0hEvEFBMEy+RheLiMcoCIY5M6hM1UMi4hEKgmHyMzUDqYh4i4JgGM1AKiJeoyAY5p35hhQEIuINCoJhstKSSUkyGtVGICIeoSAY5sx8Q2ojEBGPUBCMIC+QqqohEfEMBcEI8gN+jSMQEc+IKQjMLM/MtpjZofBr7ij7fdbM9prZPjP73JDtXzazajN7I/xzSyzliZe8gH/KVimrb+vhgZ+9SXtP/5RcT0RkuFifCB4AtjrnlgBbw5/PYmYrgXuBtcDFwK1mtmTILv/mnFsd/tkcY3niIj9z6qai/vW+Uzz+6gl+8Xr1lFxPRGS4WINgA7Ap/H4T8P4R9lkGvOyc63LODQDbgdtjvG5C5Qf8dPQO0DuQ+PmG9la3AvDjnScTfi0RkZHEGgRFzrlagPBr4Qj77AXeY2b5ZpYB3ALMHfL9/Wb2ppk9OlrV0lSbykFle6rbSPIZe6pb2V/blvDriYgMN24QmNmz4fr94T8bJnIB59x+4J+ALcCvgd3AQPjrh4FFwGqgFvj6GOW4z8x2mtnOhoaGiVx60iKDyhLdTtDTP8ihunY+WDGXlCTjJ3oqEJFpMG4QOOdudM6tHOHnSaDOzIoBwq/1o5zjO865Nc659wDNwKHw9jrn3KBzLgg8QqgdYbRybHTOVTjnKmbPnh39nUahIHNqZiA9WNfOQNDxniUFrF9exC/eqKZvIJjQa4qIDBdr1dBTwD3h9/cAT460k5kVhl/LgTuAx8Kfi4fsdjuhaqRp9840E4kdVLa3OlQVtLI0mzsvnUtzZx/bDtQl9JoiIsMlx3j8V4Efm9kngOPAnQBmVgJ82zkX6Q76MzPLB/qBTzvnToe3/7OZrQYcUAX8txjLExf54TaCRFcN7a1pJSstmbLcdIqz0yjKSuUnO09y88ri8Q8WEYmTmILAOdcE3DDC9hpCjcKRz1ePcvxHY7l+omSlJ5Pss4Q3Fu+rbmVlaTZmRnKScceaMr61/TD1bT0UZqUl9NoiIhEaWTwCMyM3wYPK+geD7D/VzqrS7DPb7ry0jKCDn2tMgYhMIQXBKBI9zcShug76BoKsGBIEC2dnUjEvlx/vPIFzLmHXFhEZSkEwitDo4sQ1Fu+tCQ0kW1mSddb2D1bM5UhDJ68db0nYtUVEhlIQjCLRM5Duq24l4E9ifn7grO23XFRMekoSP9l5ImHXFhEZSkEwikRXDe2taWNFSTY+n521PTM1mfddVMwv36ylq29glKNFROJHQTCK/ICf9p7EzDc0GHS8VdPGitKsEb+/89IyOnoH+NWeU3G/tojIcAqCUeSFRxef7oz/9NBHGzvo7h9kZUn2iN+vXZDHvPwMfrJL1UMikngKglHkR+YbSkCD8Z7wjKMrS0cOAjPjzkvLePlIM8ebuuJ+fRGRoRQEo0jkDKR7q9tIS/GxaHZg1H02rC4F4Lm3R5y+SUQkbhQEo3hnvqFEBEEry4qzSE4a/Z+/LDedDH8Sx/REICIJpiAYRWQG0sY4jy4OhhuKR2sfiDAzyvMyON7cGdfri4gMpyAYRVZaCkk+i/ugsuPNXbT3DrBylB5DQ5XnZeiJQEQSTkEwCp/PyM2I/9rFkRHFK8Z5IgCYl5/B8eYuTTchIgmlIBhDfgImnttb3UZKknFB0axx9y3Py6B3IEh9e2LXRRARb1MQjCE/M/6ji/dWt3LhnFn4k8f/py8PTz+h6iERSSQFwRjyAvGtGnLOsbemddyG4oh5eRkAHGtSg7GIJI6CYAyhqqH4VctUt3TT0tV/1tTTYynJScdncKJZTwQikjgKgjHkBVJp6xmgfzA+C8pH1iheNcEg8Cf7KMlJ55iCQEQSSEEwhnfmG4pP9dC+mlaSfMbSOeM3FEfMy1cXUhFJLAXBGAoC8R1Utre6lSWFmaSlJE34mPK8gKqGRCShFARjiPc0E5E1CKJRnpdBU2cfHb1am0BEEkNBMIb8zNDEcw0dPTGfq66th4b23gmNKB5qXr56DolIYikIxlCakw5ATUvsQbB3nKmnR1Me7kKq6ahFJFEUBGNI9ydRkOnn5OnYfwnvq2nDDJYVR/dEUB5+IjiudgIRSRAFwThKczM4ebo75vMcbuigJDudzNTkqI7LSkshNyNFXUhFJGEUBOMoy02PSxBUNXUxvyBjUseW5wdUNSQiCaMgGEdZbjrVp7sJBmObAbSqsZP5+aOvSDaW0LoECgIRSQwFwTjKcjPoGwzSEMNUE6c7+2jt7mdBweSCYF5eBtUt3XEb4SwiMlRMQWBmd5rZPjMLmlnFGPvdbGZvm1mlmT0wZHuemW0xs0Ph19xYypMIZbmhnkOxNBhXhbt+TvqJID+DwaCjpiX2KioRkeFifSLYC9wB/G60HcwsCXgIeC+wHLjbzJaHv34A2OqcWwJsDX8+p8w9EwST/yV8Jggm2UYQmYVU1UMikggxBYFzbr9z7u1xdlsLVDrnjjjn+oDHgQ3h7zYAm8LvNwHvj6U8iVCaE/olHEsQHG3swmcwN2+yjcWRQWUKAhGJv6loIygFTgz5fDK8DaDIOVcLEH4tnILyRCUeYwmONXVSkpNOavLE5xgaqmhWGv5kn54IRCQhxu3UbmbPAnNG+OpLzrknJ3ANG2Fb1F1wzOw+4D6A8vLyaA+PSaxjCWLpMQSh9ZNDC9lrmgkRib9xg8A5d2OM1zgJzB3yuQyoCb+vM7Ni51ytmRUD9WOUYyOwEaCiomJKV3Mvy03nrZq2SR3rnONoYye3rS6JqQyhLqRqLBaR+JuKqqFXgSVmtsDM/MBdwFPh754C7gm/vweYyBPGlCvLmfxYgpauftp6BmJ6IoBwEDR14tyUZqCIeECs3UdvN7OTwB8A/2lmvwlvLzGzzQDOuQHgfuA3wH7gx865feFTfBVYb2aHgPXhz+ecstz0SY8lOBpj19GIefkZdPYN0hTHNZRFRGACVUNjcc49ATwxwvYa4JYhnzcDm0fYrwm4IZYyTIWy3EjPoS6KstKiOraqMdJ1NPYggFAX0oLw9NgiIvGgkcUTUBbDWIKqxs5w19H0mMpwLkxH3djRy7e2H+alw43TVgYRib+Yngi8ojSWIGjqiqnraERZbgZm0zOWYF9NK999sYqndtfQNxBkfn4G2/7iWny+kTqEicj5RkEwARn+ZPID/kkGQeek5xgaKi0liTlZaRxrnpoupINBx5a3TvHoi1W8crSZ9JQk/qhiLoWzUvn6loO8fKSJKxYXTElZRCSxFAQTFJqOOrq/xiNdR9+/unT8nSdgbl7GlCxk3zswyPsfeon9tW2U5qTzpVuW8cHL5pKdnkJP/yCPPH+Ex149oSAQmSEUBBNUlpvB/troxhKc7uqnvWcg5obiiHl5GWw/2BCXc43llaPN7K9t43++bxkfu2I+yUnvNCWlpSRxx5oyfrjjOM2dfeQF/Akvj4gklhqLJ6gsN52TLdGNJTga6TGUP7k5hoabl59BfXsv3X2DcTnfaLYdqCc12ceHL593VghE3LV2Ln2DQX7+2smElkNEpoaCYILKctPpGwjSGMVYgnh1HY2ITFp3Ig5rKI/GOce2A/VcsSifdP/IDdxL52RxSXkOj71yXAPcRGYABcEERcYSnIiiwfhYU7jraG68nggC4fMmLgiONHZyrKmL65eOPf/f3ZeVc7ihk53HTiesLCIyNRQEEzSZBWqONnVRmpuOPzk+/8yRdQkSOfnccwdC0z1dN04Q3HpxMZmpyTy243jCyiIiU0NBMEGTGUsQ66yjw+VkpDArLTmh01Fv3V/PhUWzzjwBjSbDn8yG1SX8555aWrv6E1YeEUk8BcEERTuWwDkXtzEEEWaW0IXs23r6ebWqedyngYi715bTOxDkF29UJ6Q8IjI1FARRiGYsQXNnH+09A2fq9eNlXn5GwqaZeOFQIwNBN277QMTK0mxWlWar0VjkPKcgiEJZbgbVE3wiiKxTvGCS6xSPpjwvwInTXQxOYkrs8WzdX092egprynMmfMxda+dy4FQ7b5xoiXt5RGRqKAiiUBrFWIKqxtBf7fFsI4DQ5HP9g45TbT1xPW8w6Nh+sJ5rLpg94tiB0dx2cQnpKUk89ooajUXOVwqCKEQzlqAq3HV0vEbXaM3LT0zPoTerW2ns6JtwtVDErLQUbru4hKd319Leo0ZjkfORgiAKkS6kExlLcLSxk7LcjLh1HY1I1HTU2/bX4TO45oLZUR9719q5dPcP8tTumvF3FpFzjoIgCkMXqBlPVVNn3EYUD1WcnUayzzgW555D296uZ015LrmTmDto9dwcls6ZxY9ePRHXMonI1FAQRKE0Z2JjCZxzHGvsYkGc5hgaKjnJx7z8DCrrO+J2zrq2HvZWt0242+hwZsaG1aW8ebKVU63xbbsQkcRTEEQhkJpMXsBPdcvYQdDU2Ud7b/y7jkYsL8nmrZroZkIdS2Q08Q3LJhcEQ4/dFj6XiJw/FARRCo0lGDsIIpPNxXMw2VArS7KobunmdJwWst92oJ6S7DQuLJo16XMsKcykLDedbQfq4lImEZk6CoIoTWRQWVW4ITcRbQQAK0qyAdgXh6eC3oFBXqhs5LqlhZhNfulJM+OGpYW8UNlIT39ip8kWkfhSEEQpMqhsrJG0VY2dJPnsTC+jeFtRkgWE1hKO1Y4jzXT1DUbdbXQk1y8roqc/yO8PN8V8LhGZOgqCKJXlptM7EKRhjLEER5s6KctNJyWKgVnRyA34Kc1Jj8sTQWQRmisWxb7s5OUL8sjwJ6mdQOQ8oyCIUtkEZiE91hTfWUdHsrwkK+YngoksQhONtJQkrlxcwLYD9VM691BlfTub99RO2fVEZhoFQZTeGUswchA456hq7Irb8pSjWVGSxZHGTjp7ByZ9jiONnRxvHn8RmmjcsLSQ6pZu3q5rj9s5x7KvppUPPPx7/uw/XuPftx+ekmuKzDQKgii9M5Zg5Abjxo4+Onrjt2D9aFaUZOMcHDg1+eqh3x1sAODaC+MXBJGxCFv3J7566MCpNj7y7R0E/EmsX17EV391gO//virh1xWZaRQEUYqMJRjtiSAyB1Cig2BlaaTBePJB8GJlI/PyM86shRwPRVlprCrNPjM2IVEO1bXz4Ud24E/28cN71/HND6/hxmWF/M2T+/jprpMJvbbITKMgmISxxhIcjSxYn+A2gjlZaeQF/OyrnlwQ9A8GeflIM1cujr2ReLjrlxby2vHTNMdpnMNwlfUd3P3IDnw+47F71zG/IEBKko8HP7SGqxYX8IWf7labgUgUFASTUJoz+liCqqbEdh2NMDNWlGSxd5INxrtPtNDRO8DVCQqCoIPtB+P/VHC0sZMPPfIy4Hjs3nUsnJ155ru0lCQ2/vGlrCnP5TOPva7BbSITFFMQmNmdZrbPzIJmVjHGfjeb2dtmVmlmDwzZ/mUzqzazN8I/t8RSnqlSlps+4liCgcEgu0+0MjeBXUeHWl6SxcG6dvoGglEf+/yhRszgDxblx71cq0qzKchMjXs7wbGmTu7e+DIDQccP713H4sLMd+2T4U/m0T+5jGXFWXzqB6/xUmVjXMsgMhPF+ttqL3AH8LvRdjCzJOAh4L3AcuBuM1s+ZJd/c86tDv9sjrE8U6IsN+NdYwkO1rVzx8Mv8UJlI7deVDIl5VhZkk3/oONQffQ9dF6sbOSi0mxyMqKfbXQ8Pp9x/dLZbD/YQP9g9CE1kmDQ8akfvEbvwCD/8cnLuWCM6TCy0lLY9PG1zM/P4JP/b+eElxcV8aqYgsA5t9859/Y4u60FKp1zR5xzfcDjwIZYrjvdho4lGBgM8tBzldz6jReoPt3NNz+8hv9x04VTUo53RhhH107Q3tPP6ydaEtI+EHH90iLaewbYdex0XM737P469te28Te3LmdZcda4++cF/Dz6scvoHwzyre1H4lIGkZlqKtoISoGhE9WfDG+LuN/M3jSzR80sd7STmNl9ZrbTzHY2NDQkqqwTEhlL8NsD9Xzg4Zf4l9+8zfrlRTzz+fdwy6riKSvH/PwAAX8S+6qjayfYcaSZwaDjqiWJC4KrlhTgT/LFZZSxc46HnqukPC+D2y6e+NNWWW4GH1hTxo92nqA+zkt7iswk4waBmT1rZntH+JnoX/UjzWQWqVx/GFgErAZqga+PdhLn3EbnXIVzrmL27OhX0Yqn0vATwTe2VXLidDcPfWgND314DfmZqVNaDp/PwiOMo3sieKGykbQUH5fOGzV3Y5aZmszlC/PYuj/2BtvnDzWy+2Qrf3rtoqjWUwb402sXMTAY5JHn9VQgMppx/6tyzt3onFs5ws+TE7zGSWDukM9lQE343HXOuUHnXBB4hFA10jkvMzWZqxYXcOtFxTzz+ffwvoum7ilguBUl2bxV20YwOPEpHV6obGTtgnxSk2OfVmIs1y8t5HBD55lpuSfrwecqKc5O4441pePvPMy8/AC3XVzCf+w4nrDurCLnu6moGnoVWGJmC8zMD9wFPAVgZkN/g95OqPH5vPCDT17Ogx9aQ8EUPwUMt7wki66+QY5OcDH7U609VNZ3cNXi+PcWGi4ydUUs1UM7jjTxytFm7nvPwkkH16evW0xX3yDfffHopMshMpPF2n30djM7CfwB8J9m9pvw9hIz2wzgnBsA7gd+A+wHfuyc2xc+xT+b2R4zexO4Dvh8LOXxomgbjF8Id6e8anHiq9fm5QdYXJjJc29PPggefK6Sgkw/d11WPulzLCmaxc0r5vC9l6po6+mf9HlEZqpYew094Zwrc86lOueKnHM3hbfXOOduGbLfZufcBc65Rc65fxiy/aPOuVXOuYucc7c55zQcNEpLCmfhT/JNeCbSFysbyQ/4WTpn8quRReOGpYW8fKSJpjGm7R7N7hMtPH+okU9evTDm2VHvv34x7T0DfP/3x2I6j8hMpJHF5zl/so8L5mROaA1j5xwvVDZy5eICfL7Jr0YWjTsr5hJ08PUtB6M+9sHnKslOT+Ej6+bFXI6Vpdlce+FsvvPCUbr6Jj9jq8hMlDzdBZDYrSjO5pm3TuGcG3O5yYN1HTS093JVAscPDLe4MJM//oN5bHqpio+umzehMQAQmll0y1t1fO7GJWSmxuf/pn9+/WI+8PDv+eGO43zy6oVxOedInHNseqmKr/76AP2DjiSfkeyzM69pKUn8yZXzuffqhTEtDyoSL3oimAFWlGZxuquf2tax+8o/fyg0/uLKBI4fGMnnbriA7PQU/v7ptya8YM1Dzx0mMzWZj10xP27luHReHusW5vHI80cStq7ywGCQv31yH19++i3WLsjnT69ZxMevXMCHLy/nA2vKuPWiEhbNzuQrmw/wqR/sUpuFnBP0RDADDF3MviRn9MnuXqxsZGFB4MyaClMlOyOF/77+Av7myX38Zt8pbl45dnfbIw0d/PLNGj51zaK4T4Hx59cv4cPf3sFPd52MS5XTUK3d/dz/w9d4/lAj/+2ahXzxpqUjVsE55/jOC0f5x18d4Lb/+wIPf+TSCT8piSSCnghmgGXFszAbezH7voEgO44mZtrpibh7bTkXFs3iHzbvH/ev8Yd/e5jUZB+fuGpB3MtxxaJ8Vs/N4d+3H47bPEgAx5u6+MDDL/H7w0380wdW8VfvXTZqO4yZ8cmrF/L4fevo6hvk9m++yM+0hoJMIwXBDJDhT2ZhQYC9Y6xN8Prx03T1DSZ0WomxJCf5+Ns/XM6J5m4eHaU/v3OOb2w9xE92neRDa+clZIyGmfHn1y/m5OluNr1UFZdzvlrVzPu/+SIN7b18/xOX80cT7Op62fw8fvmZq1g9N4e/+Mlu/vqJPQmrshIZi6qGZoiVpdm8erR51O9frGzEZ7BuYeIHko3mysUFrF9exIPbKvmva8oozEo7811P/yB/+dM3eXp3DXdcUsoXbk7cxH3XLy3kxmVF/NOvD3D5gnxWlWVP+lyb99TyucffoCw3ne987DIWRLkyXeGsNH7wicv52jMH+ffth9lf28a3/7hiSqYraeroZeex0+ysaua14y0EUpNZv6yQG5YVjVnFKDOPTbTx7lxSUVHhdu7cOd3FOKds/N1hvrL5AK//zXpyA++uV7/9my8C8MSfXTnVRTtLVWMn6/9tOxtWl/K1Oy8GoL6th3v/307erG7lCzct5VPXJL43zenOPm75xvOkJvv45WeunlTPpF/vPcWnf/gaa8pzeOSPK2Juz/jVnlo+96M3KMpK43t/ctlZi+7Ew8BgkP/cU8sLhxrZdew0R8JTf/iTfKwqy6apo5eqptCU3StKsrhxWRHrlxexoiRLvZtmCDPb5Zx719oxeiKYIYY2GA+v/mnr6Wf3iRY+fd3i6SjaWeYXBPj4VQv41vYjfHTdPJJ8xic37aStp59//8il3LRizpSUIzfg5xt3X8Iffev3fOmJPfzvP1od1S+7rfvr+PPHXuPismy++ydr49LF9b2riinKTuPeTTu54+GX2PjRCtYuyIv5vAAvVTby9798iwOn2snNSOHSeXl88LK5XDY/l5Wl2aQmJ+Gc43BDJ8/ur+PZt+r4xrZD/J+th1g4O8AXblrKTSuKFAgzlJ4IZoiWrj5W//0WHnjvUj51zSIgVN3yVm0bv9pTyyPPH+VH963j8mmsGopo7+nnuq9tJystmdrWHnIzUvj2PZexvGTqe878362H+PqWg/zzBy7ig5fNHf8AYPvBBu7dtJNlxbP4/icvJystJa5lOt7Uxce+9wonm7v5lzsvYsPq6Cfbi6hq7OQrm/fzzFt1lOWm86VblnHzyjkT+oXe1NHL1gP1PPK7Ixyq76BiXi5fet8yLilP3Ky1klijPREoCGaQK7+6jbyAn2XFs3jzZCuH6jsYDM9Kurgwk82fuRp/8rnRP+BHrx7niz/bwyXlOXzro5dSOCtt/IMSYDDo+Oh3dvDa8dM8ff9VLBlj5TMItbV8/Huvsmh2Jo/du47sjPiGQERLVx/3fX8Xrxxt5i9vupA/u3ZRVH+Nt/f08+C2Sr77YhXJScanr1vMJ65aQFpK9FN1DAwG+cmuk/zrloM0tPfyvouK+eJNSynPz4j6XBPV0TvAnpOtVLd0k5rsC/2kJJ15PystmQUFmSRN0Qj5mUJB4AGfffx1nnyjhryAn1Wl2VxUls2q0mxWlWUzJyvtnHqsd87x+yNNrCnPndQvp3iqb+vhvf/neQoyU3ny/itHLc/LR5r42HdfYX5+gB/eu468Edpi4ql3YJAv/vRNfvFGDf9leRF3Vszl6iUFo5bPOcebJ1t5encNT7xeTXNXH/91TRl/edOFZzXMT1Zn7wAbf3eEjb87wkAwyEfWzePDl5ezuDC2eav6B4O8faqdN060sPtEC7tPtnCovoPxfjVlpiZzSXkOl83Po2JeLqvLc8jwq7Z7LAoCD+juG6Slu++c+6V/Pvjt2/V87Luv8qHLy/nK7auA0NNCfXsPNS09HK7v4MtP76MkJ53H71s3ZdOPO+d4cFsl337hKK3d/WT4k7huaSE3r5jDdUsLCfiTeLuunad31/D07lqON3eRkmRcc0Ehn7lhMReV5cS9THVtPfzrMwf56WsnGQw6VpVmc/slpdy2umRC/y7OOQ7WdfBCZSMvVjay40gTnX2hbrN5AT8Xl2Vz8dwcLp6bw4L8AAPBID39QXoHgvQODNI7EKS5o4/XT5xmZ9Vp3q5rxzlI8hkrSrK4bH4el83PpWJ+Xtz+d2rq6KWyvoND9R1UNXbS2t1Pe88AHb0DtPe88z4zNZn8TD8FmankZ/rJD6RSMCuV0pw05ucHmJuXQUqUiyvFk4JAZBz/+Kv9fGv7ES6em0Njey+n2nrOVK0BLJwd4PF718Xlr+to9Q8GeflIE7/ae4pn9p2isaMPf7KP4uw0jjV1keQzrliUzx9eXMJNy+ckrMpqqPr2Hp7eXcsTr59kb3UbST7j6iUF3HZxCdnpKfQPBukbdPQPBOkfDNLdP8juEy28eLiJhvbQbLQLCgJcuTiftQvyuWRuDmW56VH/EdPa3c9rx0PdYF+tOs0bJ1roGwgNFlw4O8Da+XlcNj+P+QUBcjJSyElPISs95axfyP2DQeraeqht7aGmpZualh5OnO6isr6DyvqOsxY1SkvxkZPuZ1ZacvgnhVlpyQT8yXT0DdDU0UtTRx9NnX2c7uo768kmyWfMzU1nQUGA+QUB5uZmUDArlYKAP/SamUpOesqoI9KDLrTk42QnjVQQiIyjfzDIX/18DyeauyjNSac4J42SnHRKstMpyQn9x3sutLEMBh27jp3mV3trOdrYyQ1LC3nvquJpXSTpUF07P3+9midfr6ZmjDmvCjL9XLGogKsWF3DF4vwz63/HU+/AIHurW3nlaCQcmmnrefeMs5mpyWSnpzAQDNLQ3svwRf5yMlJYPDuTJUWZLJqdyZKiWSwpzKQ4e+JP3AODQZo7+zhxupuqxk6ONnZytKmTow2dVDV10tX37gGEST5jVloywaBjMOgYGPIKsOnja7nmgsmtJ6IgEJGECwYdb9W2MRh0pCT58Cf78Cf5SEk2UpJ85GX4p2wK9KFlqmzooKalm9buflq6+s+8tnT3kWRGcU46JdlpFOekU5qTRnF2OoE4zXo7GuccLV39NHX20tDeR2NHL00dvTR29NHW04/PQrPVJif5zpq99g8vLmF+lAMXIzSOQEQSzuczVpZOfqR2Ivh8xgVFs7hgnB5hU83MyA34yQ34WVw4vWWZ/udcERGZVgoCERGPUxCIiHicgkBExOMUBCIiHqcgEBHxOAWBiIjHKQhERDzuvBxZbGYNwLFJHl4ANMaxOOcD3bM36J69IZZ7nuece9f8FOdlEMTCzHaONMR6JtM9e4Pu2RsScc+qGhIR8TgFgYiIx3kxCDZOdwGmge7ZG3TP3hD3e/ZcG4GIiJzNi08EIiIyhIJARMTjPBUEZnazmb1tZpVm9sB0lycRzOxRM6s3s71DtuWZ2RYzOxR+zZ3OMsaTmc01s+fMbL+Z7TOzz4a3z+R7TjOzV8xsd/ie/1d4+4y95wgzSzKz183sl+HPM/qezazKzPaY2RtmtjO8Le737JkgMLMk4CHgvcBy4G4zWz69pUqI7wE3D9v2ALDVObcE2Br+PFMMAH/hnFsGrAM+Hf7fdSbfcy9wvXPuYmA1cLOZrWNm33PEZ4H9Qz574Z6vc86tHjJ2IO737JkgANYClc65I865PuBxYMM0lynunHO/A5qHbd4AbAq/3wS8fyrLlEjOuVrn3Gvh9+2EfkmUMrPv2TnnOsIfU8I/jhl8zwBmVga8D/j2kM0z+p5HEfd79lIQlAInhnw+Gd7mBUXOuVoI/eIEpnmF1MQws/nAJcAOZvg9h6tI3gDqgS3OuRl/z8D/Br4ABIdsm+n37IBnzGyXmd0X3hb3e/bS4vU2wjb1nZ0hzCwT+BnwOedcm9lI/3PPHM65QWC1meUAT5jZymkuUkKZ2a1AvXNul5ldO83FmUpXOudqzKwQ2GJmBxJxES89EZwE5g75XAbUTFNZplqdmRUDhF/rp7k8cWVmKYRC4D+ccz8Pb57R9xzhnGsBfkuoXWgm3/OVwG1mVkWoWvd6M/sBM/uecc7VhF/rgScIVXHH/Z69FASvAkvMbIGZ+YG7gKemuUxT5SngnvD7e4Anp7EscWWhP/2/A+x3zv3rkK9m8j3PDj8JYGbpwI3AAWbwPTvn/so5V+acm0/ov91tzrmPMIPv2cwCZjYr8h74L8BeEnDPnhpZbGa3EKpnTAIedc79w/SWKP7M7DHgWkJT1dYBfwf8AvgxUA4cB+50zg1vUD4vmdlVwPPAHt6pO/5rQu0EM/WeLyLUSJhE6I+5Hzvn/t7M8pmh9zxUuGrofzjnbp3J92xmCwk9BUCoGv+Hzrl/SMQ9eyoIRETk3bxUNSQiIiNQEIiIeJyCQETE4xQEIiIepyAQEfE4BYGIiMcpCEREPO7/A1xIepS06k47AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(0, numPoints, numPoints)\n",
    "plt.plot(x, energies)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.09420978 0.16749818 0.1130246  0.02398223 0.09420978 0.16749818\n",
      " 0.1130246  0.02398223]\n"
     ]
    }
   ],
   "source": [
    "print(C)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Now, the nuclear motion from 9.3.2\n",
    "\n",
    "Now some derivatives for the nuclear motion molecular dynamics. Again, those could be optimized by passing overlap and so on. For example a lot of things are computed here and also re-computed in the calls of the functions inside. \n",
    "Also the computation would benefit from using symmetries to avoid computing many of the electron-electron integrals and also their derivative (it offers a roughly 8x speed improvement from this).\n",
    "I'm not going to bother, this is done in the C++ Hartree-Fock project, check that out for the details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 9.39\n",
    "def F0Deriv(t):\n",
    "    if t == 0:\n",
    "        return -1./3.\n",
    "    return (m.exp(-t) - F0(t))/(2. * t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 9.34\n",
    "def OverlapTwoCentersDeriv(alpha, beta, Ra, Rb, X):\n",
    "    difR = Ra - Rb\n",
    "    X2 = difR.dot(difR)\n",
    "    if X2 == 0:\n",
    "        return 0        \n",
    "    return -2. * alpha * beta / (alpha + beta) * X * OverlapTwoCenters(alpha, beta, Ra, Rb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 9.36\n",
    "def KineticTwoCentersDeriv(alpha, beta, Ra, Rb, X):\n",
    "    difR = Ra - Rb\n",
    "    len2 = difR.dot(difR)\n",
    "    if len2 == 0:\n",
    "        return 0\n",
    "    aplusb = alpha + beta\n",
    "    ab = alpha * beta / aplusb # sigma in the book\n",
    "    ab2 = ab * ab\n",
    "    return -4. * ab2 * X * OverlapTwoCenters(alpha, beta, Ra, Rb) + (3 * ab - 2 * ab2 * X * X) * OverlapTwoCentersDeriv(alpha, beta, Ra, Rb, X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 9.40 and 9.41\n",
    "def NuclearDeriv(alpha, beta, Ra, Rb, X, Z = 1.):\n",
    "    aplusb = alpha + beta\n",
    "    theta = 2. * m.sqrt(aplusb/m.pi)\n",
    "    X2 = X * X\n",
    "    \n",
    "    difR = Ra - Rb\n",
    "    len2 = difR.dot(difR)\n",
    "    if len2 == 0:\n",
    "        return -2. * Z * theta * OverlapTwoCenters(alpha, beta, Ra, Rb) * F0Deriv(aplusb * X2) * X * aplusb # 9.40\n",
    "\n",
    "    alpha2 = alpha * alpha\n",
    "    beta2 = beta * beta\n",
    "    \n",
    "    t1 = alpha2 * X2 / aplusb # 9.42a\n",
    "    t2 = beta2 * X2 / aplusb  # 9.42b\n",
    "    \n",
    "    # 9.41\n",
    "    return -Z * theta * (OverlapTwoCentersDeriv(alpha, beta, Ra, Rb, X) * (F0(t1) + F0(t2)) + 2. / aplusb * X * OverlapTwoCenters(alpha, beta, Ra, Rb) * (F0Deriv(t1) * alpha2 + F0Deriv(t2) * beta2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 9.47 - there is some renaming to have the same names as in TwoElectronTwoCenter\n",
    "def TwoElectronTwoCenterDeriv(alpha, beta, gamma, delta, Ra, Rb, Rc, Rd, X):   \n",
    "    RP = Rp(alpha, gamma, Ra, Rc) # 9.38\n",
    "    RQ = Rp(beta, delta, Rb, Rd) # 9.45        \n",
    "    alphaplusgamma = alpha + gamma\n",
    "    betaplusdelta = beta + delta\n",
    "    PQ = RP - RQ\n",
    "    PQ2 = PQ.dot(PQ)\n",
    "    cdiv = alphaplusgamma * betaplusdelta / (alphaplusgamma+betaplusdelta)\n",
    "\n",
    "    t = cdiv * PQ2 # 9.44\n",
    "    rho = 2. * m.sqrt(cdiv / m.pi) # 9.46\n",
    "    F0t = F0(t)\n",
    "        \n",
    "    return rho * (\n",
    "        (OverlapTwoCentersDeriv(alpha, gamma, Ra, Rc, X) * OverlapTwoCenters(beta, delta, Rb, Rd) + \n",
    "        OverlapTwoCenters(alpha, gamma, Ra, Rc) * OverlapTwoCentersDeriv(beta, delta, Rb, Rd, X)) * F0t +\n",
    "        OverlapTwoCenters(alpha, gamma, Ra, Rc) * OverlapTwoCenters(beta, delta, Rb, Rd) * F0Deriv(t) * cdiv * 2 * PQ2 / X\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "HDeriv = np.zeros((basisSize * 2, basisSize * 2))\n",
    "ODeriv = np.zeros((basisSize * 2, basisSize * 2))\n",
    "QDeriv = np.zeros((basisSize*2, basisSize*2, basisSize*2, basisSize*2))\n",
    "QtDeriv = np.zeros((basisSize*2, basisSize*2, basisSize*2, basisSize*2))\n",
    "FDeriv = np.zeros((2*basisSize, 2*basisSize))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is the interesting part:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = 1.35\n",
    "Xprev = X\n",
    "R1 = np.array([X, 0, 0])\n",
    "\n",
    "for i in range(basisSize):\n",
    "    a = alpha[i]\n",
    "    for j in range(basisSize):        \n",
    "        b = alpha[j]\n",
    "        O[i, j] = OverlapTwoCenters(a, b, R0, R0)\n",
    "        O[i, 4 + j] = OverlapTwoCenters(a, b, R0, R1)\n",
    "        O[4 + i, j] = O[i, 4 + j]\n",
    "        O[4 + i, 4 + j] = O[i, j]  \n",
    "            \n",
    "# reinitialize Cs\n",
    "\n",
    "v = 1. / m.sqrt(O.sum())\n",
    "\n",
    "C = np.array([v, v, v, v, v, v, v, v])\n",
    "Cprev = C.copy()    \n",
    "\n",
    "M = 1836.5 * mass\n",
    "\n",
    "numPoints = 43\n",
    "\n",
    "Nucleardt = numPoints * dt # 4.3 (atomic units)\n",
    "Nucleardt2 = Nucleardt * Nucleardt\n",
    "\n",
    "numNuclearPoints = 300\n",
    "distances = np.zeros((numNuclearPoints))\n",
    "\n",
    "# friction stuff for nuclei\n",
    "\n",
    "Nucleargamma = 15. # set it to zero if you want no 'friction'\n",
    "\n",
    "MassMinusGamma = M - 0.5 * Nucleargamma * Nucleardt\n",
    "MassPlusGamma = M + 0.5 * Nucleargamma * Nucleardt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "for nuclearCycle in range(numNuclearPoints):\n",
    "       \n",
    "    # recompute each time since the atoms positions change each nuclear cycle\n",
    "        \n",
    "    for i in range(basisSize):\n",
    "        a = alpha[i]\n",
    "        for j in range(basisSize):        \n",
    "            b = alpha[j]\n",
    "            O[i, j] = OverlapTwoCenters(a, b, R0, R0)\n",
    "            O[i, 4 + j] = OverlapTwoCenters(a, b, R0, R1)\n",
    "            O[4 + i, j] = O[i, 4 + j]\n",
    "            O[4 + i, 4 + j] = O[i, j]  \n",
    "        \n",
    "            H[i, j] = KineticTwoCenters(a, b, R0, R0) + Nuclear(a, b, R0, R0, R0) + Nuclear(a, b, R0, R0, R1)\n",
    "            H[i, 4 + j] = KineticTwoCenters(a, b, R0, R1) + Nuclear(a, b, R0, R1, R0) + Nuclear(a, b, R0, R1, R1)\n",
    "            H[4 + i, j] =  H[i, 4 + j]\n",
    "            H[4 + i, 4 + j] = H[i, j] \n",
    "            \n",
    "            for k in range(basisSize):\n",
    "                c = alpha[k]\n",
    "                for l in range(basisSize):\n",
    "                    d = alpha[l]                                    \n",
    "                    Q[i, j, k, l]=TwoElectronTwoCenter(a, b, c, d, R0, R0, R0, R0)\n",
    "                    Q[i, j, k, basisSize+l]=TwoElectronTwoCenter(a, b, c, d, R0, R0, R0, R1)\n",
    "                    Q[i, j, basisSize+k, l]=TwoElectronTwoCenter(a, b, c, d, R0, R0, R1, R0)\n",
    "                    Q[i, j, basisSize+k, basisSize+l]=TwoElectronTwoCenter(a, b, c, d, R0, R0, R1, R1)\n",
    "                    Q[i, basisSize+j, k, l]=TwoElectronTwoCenter(a, b, c, d, R0, R1, R0, R0)\n",
    "                    Q[i, basisSize+j, k, basisSize+l]=TwoElectronTwoCenter(a, b, c, d, R0, R1, R0, R1)\n",
    "                    Q[i, basisSize+j, basisSize+k, l]=TwoElectronTwoCenter(a, b, c, d, R0, R1, R1, R0)\n",
    "                    Q[i, basisSize+j, basisSize+k, basisSize+l]=TwoElectronTwoCenter(a, b, c, d, R0, R1, R1, R1)\n",
    "                    Q[basisSize+i, j, k, l]=TwoElectronTwoCenter(a, b, c, d, R1, R0, R0, R0)\n",
    "                    Q[basisSize+i, j, k, basisSize+l]=TwoElectronTwoCenter(a, b, c, d, R1, R0, R0, R1)\n",
    "                    Q[basisSize+i, j, basisSize+k, l]=TwoElectronTwoCenter(a, b, c, d, R1, R0, R1, R0)\n",
    "                    Q[basisSize+i, j, basisSize+k, basisSize+l]=TwoElectronTwoCenter(a, b, c, d, R1, R0, R1, R1)\n",
    "                    Q[basisSize+i, basisSize+j, k, l]=TwoElectronTwoCenter(a, b, c, d, R1, R1, R0, R0)\n",
    "                    Q[basisSize+i, basisSize+j, k, basisSize+l]=TwoElectronTwoCenter(a, b, c, d, R1, R1, R0, R1)\n",
    "                    Q[basisSize+i, basisSize+j, basisSize+k, l]=TwoElectronTwoCenter(a, b, c, d, R1, R1, R1, R0)\n",
    "                    Q[basisSize+i, basisSize+j, basisSize+k, basisSize+l]=TwoElectronTwoCenter(a, b, c, d, R1, R1, R1, R1)       \n",
    "    \n",
    "    \n",
    "    for p in range(2*basisSize):\n",
    "        for q in range(2*basisSize):\n",
    "            for r in range(2*basisSize):\n",
    "                for s in range(2*basisSize):\n",
    "                    Qt[p, q, r, s] = 2. * Q[p, q, r, s] - Q[p, r, s, q]                    \n",
    "                    \n",
    "                 \n",
    "    # compute derivatives, similar for the other ones, but with removed computations for the same center for overlap and h derivatives, for the values that are zero\n",
    "    \n",
    "    for i in range(basisSize):        \n",
    "        a = alpha[i]\n",
    "        for j in range(basisSize):        \n",
    "            b = alpha[j]\n",
    "            ODeriv[i, 4 + j] = OverlapTwoCentersDeriv(a, b, R0, R1, X)\n",
    "            ODeriv[4 + i, j] = ODeriv[i, 4 + j]\n",
    "            \n",
    "            HDeriv[i, j] = NuclearDeriv(a, b, R0, R0, X)            \n",
    "            HDeriv[i, 4 + j] = NuclearDeriv(a, b, R0, R1, X) + KineticTwoCentersDeriv(a, b, R0, R1, X)\n",
    "            HDeriv[4 + i, j] = HDeriv[i, 4 + j]\n",
    "            HDeriv[4 + i, 4 + j] = HDeriv[i, j] \n",
    "    \n",
    "            for k in range(basisSize):\n",
    "                c = alpha[k]\n",
    "                for l in range(basisSize):\n",
    "                    d = alpha[l]                                    \n",
    "                    QDeriv[i, j, k, l]=TwoElectronTwoCenterDeriv(a, b, c, d, R0, R0, R0, R0, X)\n",
    "                    QDeriv[i, j, k, basisSize+l]=TwoElectronTwoCenterDeriv(a, b, c, d, R0, R0, R0, R1, X)\n",
    "                    QDeriv[i, j, basisSize+k, l]=TwoElectronTwoCenterDeriv(a, b, c, d, R0, R0, R1, R0, X)\n",
    "                    QDeriv[i, j, basisSize+k, basisSize+l]=TwoElectronTwoCenterDeriv(a, b, c, d, R0, R0, R1, R1, X)\n",
    "                    QDeriv[i, basisSize+j, k, l]=TwoElectronTwoCenterDeriv(a, b, c, d, R0, R1, R0, R0, X)\n",
    "                    QDeriv[i, basisSize+j, k, basisSize+l]=TwoElectronTwoCenterDeriv(a, b, c, d, R0, R1, R0, R1, X)\n",
    "                    QDeriv[i, basisSize+j, basisSize+k, l]=TwoElectronTwoCenterDeriv(a, b, c, d, R0, R1, R1, R0, X)\n",
    "                    QDeriv[i, basisSize+j, basisSize+k, basisSize+l]=TwoElectronTwoCenterDeriv(a, b, c, d, R0, R1, R1, R1, X)\n",
    "                    QDeriv[basisSize+i, j, k, l]=TwoElectronTwoCenterDeriv(a, b, c, d, R1, R0, R0, R0, X)\n",
    "                    QDeriv[basisSize+i, j, k, basisSize+l]=TwoElectronTwoCenterDeriv(a, b, c, d, R1, R0, R0, R1, X)\n",
    "                    QDeriv[basisSize+i, j, basisSize+k, l]=TwoElectronTwoCenterDeriv(a, b, c, d, R1, R0, R1, R0, X)\n",
    "                    QDeriv[basisSize+i, j, basisSize+k, basisSize+l]=TwoElectronTwoCenterDeriv(a, b, c, d, R1, R0, R1, R1, X)\n",
    "                    QDeriv[basisSize+i, basisSize+j, k, l]=TwoElectronTwoCenterDeriv(a, b, c, d, R1, R1, R0, R0, X)\n",
    "                    QDeriv[basisSize+i, basisSize+j, k, basisSize+l]=TwoElectronTwoCenterDeriv(a, b, c, d, R1, R1, R0, R1, X)\n",
    "                    QDeriv[basisSize+i, basisSize+j, basisSize+k, l]=TwoElectronTwoCenterDeriv(a, b, c, d, R1, R1, R1, R0, X)\n",
    "                    QDeriv[basisSize+i, basisSize+j, basisSize+k, basisSize+l]=TwoElectronTwoCenterDeriv(a, b, c, d, R1, R1, R1, R1, X)       \n",
    "     \n",
    "    \n",
    "    for p in range(2*basisSize):\n",
    "        for q in range(2*basisSize):\n",
    "            for r in range(2*basisSize):\n",
    "                for s in range(2*basisSize):\n",
    "                    QtDeriv[p, q, r, s] = 2. * QDeriv[p, q, r, s] - QDeriv[p, r, s, q]\n",
    "\n",
    "                    \n",
    "    # the electronic loop - it's identical with the one above that's done only for electronic part, with the energy computation moved out of it\n",
    "\n",
    "    for cycle in range(numPoints):\n",
    "\n",
    "        # Fock matrix computation\n",
    "        #for i in range(2*basisSize):\n",
    "        #    for j in range(2*basisSize):\n",
    "        #        F[i, j] = H[i, j]\n",
    "        #        for k in range(2*basisSize):\n",
    "        #            for l in range(2*basisSize):\n",
    "        #                F[i, j] += Qt[i, k, j, l] * C[k] * C[l]\n",
    "        \n",
    "        F = H + np.einsum('ikjl,k,l', Qt, C, C)\n",
    "                \n",
    "        # Verlet for electrons\n",
    "                    \n",
    "        # compute Ct - 9.31, but with friction force added\n",
    "            \n",
    "        Ct = (2. * mass * C - massMinusGamma * Cprev - 4. * F.dot(C) * dt2) / massPlusGamma\n",
    "                \n",
    "        # determine lambda - see 9.32 - but be careful, it's wrong! Get it from 9.28 by replacing C[r] = Ct[r] - lambda * h^2 * sum(S[r, s]*C[s]), h^4 and h^2 are missing (here h is dt)\n",
    "     \n",
    "        OC = O.dot(C)\n",
    "        OCt = O.dot(Ct)\n",
    "        OOC = O.dot(OC)\n",
    "    \n",
    "        a = OOC.dot(OC) * dt4\n",
    "        b = -2. * OC.dot(OCt) * dt2\n",
    "        c = OCt.dot(Ct) - 1.\n",
    "        \n",
    "        delta = b*b - 4.*a*c\n",
    "        if delta < 0:\n",
    "            print(\"Delta negative!\")\n",
    "            break\n",
    "        \n",
    "        sdelta = m.sqrt(delta)\n",
    "        lam1 = (-b-sdelta) / (2. * a)\n",
    "        lam2 = (-b+sdelta) / (2. * a)\n",
    "\n",
    "        if lam1 < 0:\n",
    "            lam = lam2\n",
    "        else:\n",
    "            lam = lam1\n",
    "    \n",
    "        # now adjust the newly computed Cs    \n",
    "    \n",
    "        Ct -= lam * dt2 * OC\n",
    "        \n",
    "        # switch to the next step\n",
    "        Cprev = C\n",
    "        C = Ct            \n",
    "                        \n",
    "        # end electronic loop\n",
    "        \n",
    "\n",
    "    # compute energy\n",
    "    #Eg = C.dot(H + F).dot(C) + 1. / X\n",
    "    #print(Eg)\n",
    "                        \n",
    "    # verlet for nuclear\n",
    "                         \n",
    "    #for i in range(2*basisSize):\n",
    "    #    for j in range(2*basisSize):\n",
    "    #        FDeriv[i, j] = HDeriv[i, j]\n",
    "    #        for k in range(2*basisSize):\n",
    "    #            for l in range(2*basisSize):\n",
    "    #                FDeriv[i, j] += QtDeriv[i, k, j, l] * C[k] * C[l]\n",
    "        \n",
    "    FDeriv = HDeriv + np.einsum('ikjl,k,l', QtDeriv, C, C)\n",
    "        \n",
    "    lam *= massPlusGamma * 0.5 # correction - there is a 2 factor in the lambda, also needs to be adjusted because of the massPlusGamma division\n",
    "                \n",
    "    # C.dot(ODeriv).dot(C) can be reduced to a quarter, knowing that diagonal sectors are zero\n",
    "    LambdaFactor = 0\n",
    "    for i in range(basisSize):\n",
    "        for j in range(basisSize):\n",
    "            LambdaFactor += ODeriv[basisSize + i, j] * C[basisSize + i] * C[j]\n",
    "            LambdaFactor += ODeriv[i, basisSize + j] * C[i] * C[basisSize + j]\n",
    "        \n",
    "    FNuclear = C.dot(HDeriv + FDeriv).dot(C) - 1. / (X * X) + lam * LambdaFactor\n",
    "    \n",
    "    Xnew = (2. * M * X - MassMinusGamma * Xprev - 2. * FNuclear * Nucleardt2) / MassPlusGamma\n",
    "        \n",
    "    # switch to the next step\n",
    "    Xprev = X\n",
    "    X = Xnew\n",
    "    \n",
    "    # also update the nuclei position\n",
    "    R1 = np.array([X, 0, 0]) \n",
    "    \n",
    "    distances[nuclearCycle] = X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABIcAAAEvCAYAAADfBqG/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABJ+klEQVR4nO3dd3hc1aH14bWnqPduy3KTuw1ucsWAARtMiyH0FkgIJKGk3Zub3C+FkOQmIYV0IIQOAUIg9GoTOthGLrjIvUtWl9XblP39oXENLtiSzozm9z6Pnpk5Z8qyORx7lvfex1hrBQAAAAAAgOjkcjoAAAAAAAAAnEM5BAAAAAAAEMUohwAAAAAAAKIY5RAAAAAAAEAUoxwCAAAAAACIYpRDAAAAAAAAUczjdIBPk5WVZQcPHux0DAAAAAAAgD5j6dKlNdba7IO3h2U5NHjwYBUXFzsdAwAAAAAAoM8wxmz/tO1MKwMAAAAAAIhilEMAAAAAAABRjHIIAAAAAAAgih2xHDLGPGCMqTLGrD7C86YYYwLGmIs/62sBAAAAAADgjKMZOfSQpHmHe4Ixxi3pDkmvf9bXAgAAAAAAwDlHLIeste9KqjvC026V9IykqmN4LQAAAAAAABxy3GsOGWPyJV0o6Z7jjwMAAAAAAIDe1B0LUv9e0nettYHjeRNjzI3GmGJjTHF1dXU3xAIAAAAAAMCReLrhPYokPWmMkaQsSecYY/zW2uc+y5tYa++VdK8kFRUV2W7IBQAAAAAAgCM47nLIWjtkz31jzEOSXvqsxRAQThpafVq9q0E76lolSUaSy2U0JCtRY/qlKDG2OzpVAAAAAADCwxG/5RpjnpA0W1KWMaZU0m2SvJJkrT3sOkOf9lpr7f3HmRnoVtZaLd2+W09+vFOLt9ZqZ13bIZ9rjFSYnaTpQzN0yeQCnTggVaFRcwAAAAAARCRjbfjN4CoqKrLFxcVOx0Af1+EP6KmPd+qxRTu0vrJJSbEenTIiS+PyUzWuf6oKc5LkMpK1kj9gtbGqSavKGrSqtEEfbK5Ruy+oUXnJurSoQFdMHaj4GLfTvyQAAAAAAA7JGLPUWlv0H9sphxCNPthUox8+t1pbalp0Qn6qrpo2UOeP73/UU8Ya23168ZNdeqq4VJ/srFe/1Dj9z7yRmj8+Xy4XI4kAAAAAAOGHcgiQVNPcoZ+9VKLnVuzSoMwE/WT+OJ06Ivu43nPxllr97OW1WlXWoPEDUvWT+eM0viCtewIDAAAAANBNKIcQ9VaW1uuGR4q1u8Wnr84u1E2zCxXn7Z6pYMGg1XMrynTHa+tU29yp/z5rpG48eSijiAAAAAAAYeNQ5RCXXUJUeGnlLv33Pz9RZmKsnrv5JI3pn9Kt7+9yGX1+0gCdMSpX3/vXSv3y1XX6YFONfnvpeOUkx3XrZwEAAAAA0J1cTgcAepK1Vn9YuFG3PL5c4/qn6vlbur8Y2l9qgld3XTVJP7/wBC3ZWqdz//i+VpU29NjnAQAAAABwvCiH0KfduWCDfrdwgz4/KV9/v2GaspJie/wzjTG6ctpAvXDLLMW4Xbrs3o/09vqqHv9cAAAAAACOBeUQ+qy73t6kP/17ky6fUqDfXjJesZ7evdT8yLxkPXvTTA3OTNT1DxfrqeKdvfr5AAAAAAAcDcoh9EkPfbBVv3ptveZP6K//u/AEGePMwtA5KXH6x1ema2Zhpv7n6ZW6770tjuQAAAAAAOBQKIfQ5zy/okw/frFEZ47J1W8uGS+3w1cMS47z6oHrpuicE/L0s5fX6tGPtjmaBwAAAACA/XG1MvQpa3Y16LvPrNTUIRn605UT5XWHR//pdbv0h8snqtO/TD98fo1iPC5dNmWg07EAAAAAAGDkEPqO+tZOffWxpUqLj9FfrpzU62sMHYnX7dJfrpqoU0dk63v/WqVnl5c6HQkAAAAAAMoh9A2BoNXXn1yhyoYO3X31JGUn9/xVyY5FrMetv14zWdOHZOo7/1ypDzfVOB0JAAAAABDlKIfQJ/x+4Qa9u6Fat88fq4kD052Oc1hxXrf++oXJGpqdqK88tlQbK5ucjgQAAAAAiGKUQ4h4xdvq9Oe3NumSyQN0xdTIWMcnJbRIdazHrS8+9LGqmzqcjgQAAAAAiFKUQ4horZ1+/fc/P1F+Wrxu+9xYp+N8JgPSE3T/tUWqae7Qlx8pVrsv4HQkAAAAAEAUohxCRLvj1XXaVtuqX188XkmxkXfxvfEFafr9ZRP1yc56/fC51bLWOh0JAAAAABBlKIcQsT7cVKOHP9qu62YO1ozCTKfjHLN54/L09dOH6Z9LS/XkxzudjgMAAAAAiDKUQ4hIzR1+fefplRqSlajvzhvldJzj9o05I3TKiGzd9vwafbKz3uk4AAAAAIAoQjmEiPSnNzeqrL5Nv7nkRMXHuJ2Oc9zcLqM/XDZB2cmx+tpjS1XX0ul0JAAAAABAlKAcQsTZVNWs+9/fqkuLBmjyoAyn43Sb9MQY3XP1ZNW0dOrbT61g/SEAAAAAQK+gHEJEsdbq9hfXKD7Grf/pA9PJDnbCgFT94NzRent9tR78YJvTcQAAAAAAUYByCBHl9TWVem9jjb49d4SykmKdjtMjrpk+SHNG5+qXr65Tya5Gp+MAAAAAAPo4yiFEjHZfQD99qUQjc5N1zfRBTsfpMcYY/eriE5WW4NWtTyxTW2fA6UgAAAAAgD6McggR46/vbFFZfZtunz9WHnffPnQzEmP0u8smaEtNi37yUonTcQAAAAAAfVjf/oaNPqO2uUP3vrtZZ4/L0/ShmU7H6RUnDcvSV04p1BNLduitdVVOxwEAAAAA9FGUQ4gIf3lrs9p8Af3XmSOdjtKrvjV3uEblJeu7z6xUfSuXtwcAAAAAdD/KIYS9svo2PbZouy6ePEDDcpKcjtOrYj1u/eaS8apr6dSPX1jjdBwAAAAAQB9EOYSw94eFGyRJ35gzwuEkzhiXn6pbTh+m51bs0murK5yOAwAAAADoYyiHENY2VTXr6aWlumbGIOWnxTsdxzE3nzZM4/JT9P1nV6m2ucPpOAAAAACAPoRyCGHtzgXrFe9166bZhU5HcZTX7dJvL5mgxnafbn+Rq5cBAAAAALoP5RDC1rqKRr2yqkLXnzxUmUmxTsdx3Mi8ZN182jC98Mku/XtdpdNxAAAAAAB9BOUQwtZdb21WYoxb1580xOkoYeOm2cM0IjdJP3h2tZo7/E7HAQAAAAD0AZRDCEvbalr00spdunrGIKUmeJ2OEzZiPC794vMnqryxXb9+bZ3TcQAAAAAAfcARyyFjzAPGmCpjzOojPG+KMSZgjLl4v23zjDHrjTGbjDHf647AiA5/fXezPG6Xrp/FqKGDTR6UrmtnDNYji7Zr6fY6p+MAAAAAACLc0YwcekjSvMM9wRjjlnSHpNcP2vYXSWdLGiPpCmPMmGNOiqhR0dCup5eW6rKiAuUkxzkdJyx956yR6p8ar+89s0q+QNDpOAAAAACACHbEcsha+66kIw1PuFXSM5Kq9ts2VdIma+0Wa22npCclzT/WoIge9767RUEr3XjKUKejhK3EWI9u/9xYbaxq1gPvb3U6DgAAAAAggh33mkPGmHxJF0q656Bd+ZJ27ve4NLQNOKTa5g49sWSH5k/or4KMBKfjhLU5Y3I1Z3SOfr9wo3bVtzkdBwAAAAAQobpjQerfS/qutTZw0HbzKc+1h3oTY8yNxphiY0xxdXV1N8RCJHp00Xa1+QK6aXah01Eiwm3nj5WV1U9fKnE6CgAAAAAgQnVHOVQk6UljzDZJF0u6yxhzgbpGChXs97wBknYd6k2stfdaa4ustUXZ2dndEAuRpsMf0GOLtuv0UTkalpPsdJyIUJCRoFtPH65XV1forfVVR34BAAAAAAAHOe5yyFo7xFo72Fo7WNLTkm6y1j4n6WNJw40xQ4wxMZIul/TC8X4e+q4XPylXTXOnvnQSVyj7LL588hANzU7Ubc+vUbvv4AF8AAAAAAAc3tFcyv4JSR9JGmmMKTXGXG+M+aox5quHe5211i/pFnVdwWytpKestWu6IzT6Hmut7n9/q0bmJuukYZlOx4kosR63fjp/nHbUterutzc7HQcAAAAAEGE8R3qCtfaKo30za+11Bz1+RdIrnz0Wos2iLXVaW96oOy46QcZ82nJVOJyThmXp/PH9dfc7m3XhxHwNzkp0OhIAAAAAIEJ0x5pDwHF74IOtykiM0fwJXNDuWP3g3NGKcbv0oxfWyNpDrv0OAAAAAMABKIfguO21LVq4tlJXTRuoOK/b6TgRKzclTt+aO0LvbqjWa6srnI4DAAAAAIgQlENw3EMfbpPHZXTN9EFOR4l4184YpFF5ybr9xRK1dPidjgMAAAAAiACUQ3BUa6dfTxeX6twT+iknJc7pOBHP43bp/y4cp4rGdv3x3xudjgMAAAAAiACUQ3DUi5/sUlOHX1czaqjbTB6UoYsmDdAD72/V1poWp+MAAAAAAMIc5RAc9ffFOzQyN1mTB6U7HaVP+e68kYpxu/Szl0qcjgIAAAAACHOUQ3DMytJ6rSxt0FXTB3L5+m6WkxKnW88YrjfXVent9VVOxwEAAAAAhDHKITjm8cU7FO9164KJXL6+J3zxpMEanJmgn7xUIl8g6HQcAAAAAECYohyCIxrbfXp+xS7Nn9BfKXFep+P0SbEet3543hhtqW7Rwx9uczoOAAAAACBMUQ7BEc8uK1ObL6CrprEQdU86fVSOTh2RrT8s3Kia5g6n4wAAAAAAwhDlEHqdtVZ/X7xdJw5I1QkDUp2O06cZY/TD88aozRfQb15f73QcAAAAAEAYohxCr1u2Y7c2VDbryqkDnY4SFYblJOnamYP1j+KdWlXa4HQcAAAAAECYoRxCr3vq41IlxLh13vj+TkeJGl8/Y7gyEmJ0+4trZK11Og4AAAAAIIxQDqFXtXb69dLKXTrnhH5KivU4HSdqpMZ79Z2zRqp4+2698Mkup+MAAAAAAMII5RB61SurKtTSGdClRQVOR4k6lxQVaFx+in7xyjq1dvqdjgMAAAAACBOUQ+hV/yzeqcGZCZoyON3pKFHH7TL68fljVdHYrrvf3ux0HAAAAABAmKAcQq/ZXtuixVvrdElRgYwxTseJSkWDM/S58f3113e3aGddq9NxAAAAAABhgHIIvebppaVyGenzk/KdjhLV/vecUXIbo5+/stbpKAAAAACAMEA5hF4RCFo9s7RUs4Znq19qvNNxolq/1HjdNLtQr66u0Ieba5yOAwAAAABwGJeLQq/4cHONdjW06/+dO9rpKJB0wylD9Y/inbr9hRK9/PVZ8rjpidE3WWvV2OZXXWun2n0BdfqD6gwEZSTFetyK8bgU73UrIylGiTFuprwCAAAgKlEOoVc8vbRUqfFezRmd63QUSIrzuvWDc0frq48t0+NLdugLMwY7HQk4ZtZa7axr06qyBm2qata22hZtrWnRrvo21bV0yh+0R/U+sR6XspJi1T8tTgMzEjUwI0FDshM1pl+KhmQlyu2iOAIAAEDfRDmEHtfc4dfrayp00aQBivO6nY6DkLPG5mlmYaZ++8YGnX9if6UnxjgdCTgqgaDV6rIGvb+pRou21GpVWYPqW3179/dPjdPgrETNHpmtrKRYZSbFKiPRq3hv10ghr9sla6VOf1Ad/qBaOv3a3dKp2pZO1TR1qHR3mz7YVKNnGtv3vmec16WReSmaWJCmosHpKhqUobzUOCd++QAAAEC3oxxCj3tjTYXafUFdOJGFqMOJMUY/On+Mzv3j+7pzwQb99IJxTkcCDqm106+31lXr1dXlem9jjRrausqgkbnJmjc2TycMSNWJ+WkanpvUbSV0uy+gzdXNWlvepLXljVpd1qAnP96hhz7cJkkanJmgk4dn6+ThWZpRmKnkOG+3fC4AAADQ2yiH0OOeXV6mAenxmjwo3ekoOMiovBRdPW2gHl20XVdOG6jR/VKcjgTs5Q8E9c6Gaj29tFRvra9Suy+ozMQYnTkmV7OGZ2lmYZayk2N77PPjvG6N7Z+qsf1T927zBYJas6tRxdvq9NHmWj2zrFSPLtour9to+tBMnTk2T3NH5zKqCAAAABHFWHt0azH0pqKiIltcXOx0DHSDqsZ2Tf/Fm7r5tGH6rzNHOh0Hn6K+tVOn/eZtjcxL1hM3TGdBXjiudHernliyQ08vLVVlY4eykmJ0zgn9dPa4fpo6JCOs1v7p9Ae1bMduvbWuSm+UVGprTYskaergDF04KV/nnNBPqfGMKAIAAEB4MMYstdYW/cd2yiH0pPve26KfvbxWC799qoblJDkdB4fw6KLt+uFzq3XXVZN0zgn9nI6DKLViZ73+9t4Wvba6QtZazR6Zo8umFOj0UTnyRsAV9ay12lzdrNdWV+jZ5WXaXN2iGI9Lc0bn6MKJA3TqiGzFeML/1wEAAIC+i3IIjjjvT+/JZYxeuGWW01FwGIGg1bl/fE9N7X69+V+nsnA4etWHm2v0+4UbtWRrnZLjPLpy6kBdO3Ow+qfFOx3tmFlrtaqsQf9aVqYXP9ml2pZOpSd4df74/rpy2kCNymMKJwAAAHof5RB63cbKJs393bv60Xlj9KVZQ5yOgyNYtKVWl9+7SN+aM0LfmDPc6TiIAku21unOBeu1aEudclNideMphbpsSoGSYvvWcni+QFDvbazWv5aV6Y2SSnX6g5o6OEPXzBiks8bmMZoIAAAAveZQ5VDf+hs4wsqzy8vkdhmdP76/01FwFKYPzdS5J/TT3e9s0sVFA5QfwaM2EN42VTXrF6+s1ZvrqpSdHKvbzh+jK6YO7LMj1rxul04flavTR+Vqd0un/rl0px5btEO3PrFc2cmxumLqQF05dSCLWAMAAMAxjBxCjwgGrU7+1VsalpOkh7801ek4OEqlu1s15853dMboXP3lyklOx0Efs7ulU394c6MeW7Rd8V63bj59mK6bObjPlkKHEwxavbOhWo98tE1vb6iWyxjNG5unL588RBMHcmVHAAAA9AxGDqFXLd+5W2X1bfqvM0c4HQWfwYD0BH311EL9fuFGXTO9VtOHZjodCX2AtVb/XFqqX7yyVg1tPl05baC+OWeEspJ67jL04c7lMjptVI5OG5Wj7bUt+vviHXpyyQ69vKpcUwdn6MZThur0UTlyhdGV2QAAANB3MXIIPeLHL6zR40t2aOkP5ig5jss4R5K2zoDm3PmOUuK9eunWWWF12XBEng2VTfrBs6u1ZFudpgxO108vGMdizIfQ3OHXPz7eqQfe36qy+jYVZifqxlOG6oKJ+Yr1RN/oKgAAAHS/Q40cYhVMdLtA0OqlleU6fWQOxVAEio9x6/+dM1pryxv1+JIdTsdBhPIFgvrjmxt17h/f04aqJv3qohP1jxtnUAwdRlKsR9fPGqK3vzNbf7h8gmI9bn33mVWadcdbuuvtTWps9zkdEQAAAH3UEcshY8wDxpgqY8zqQ+yfb4xZaYxZYYwpNsbM2m/fN4wxq40xa4wx3+zG3Ahji7fUqqa5g4WoI9g5J+RpxtBM/eb19apt7nA6DiLMuopGXXjXB7pzwQadPa6f/v1fs3XplAKmSB0lr9ul+RPy9fLXZ+mx66dpVF6yfvXaep30i3/rV6+tUw3/TwIAAKCbHc3IoYckzTvM/jcljbfWTpD0JUn3SZIxZpykGyRNlTRe0nnGGK6PHQVeXFmuhBi3Th+V43QUHCNjjH4yf6xaOvy647V1TsdBhAgGre55Z7PO/9P7Kq9v1z1XT9Ifr5iojMQYp6NFJGOMZg3P0qPXT9NLt87SKSOydfc7m3XSL/+t255frdLdrU5HBAAAQB9xxHLIWvuupLrD7G+2+xYuSpS05/5oSYusta3WWr+kdyRdeJx5EeZ8gaBeXV2uuWNyFR/DGhmRbHhusq4/eYieKi7V0u2HPAUAkqSKhnZdff9i/fLVdZozOldvfOsUzRvXz+lYfca4/FT95apJWvjtUzV/Qn89vmSHZv/6bX37qRXaVNXkdDwAAABEuG5Zc8gYc6ExZp2kl9U1ekiSVks6xRiTaYxJkHSOpILu+DyEr/c31ai+1afzT2RKWV/w9dOHq19qnH7w3Br5A0Gn4yBMLSip1Lw/vKvlO+p1x0Un6K6rJikziq9E1pMKs5P0q4vH653vnKYvzBisV1dVaO7v3tVXHi3WJzvrnY4HAACACNUtl7K31j4r6VljzCmSfippjrV2rTHmDkkLJDVL+kSS/1DvYYy5UdKNkjRw4MDuiAUHvPjJLiXHeXTyiCyno6AbJMZ69MPzxuimvy/TY4u267qThjgdCWHEHwjqN29s0D3vbNa4/BT94fKJKsxOcjpWVOifFq8fnT9Gt5w+TA99sFUPfbhNr6+p1KxhWbppdqFmFGbKGNZ4QmTq8AfU1O5Xc7tfzR1+Nbb71NIRkC8QlC8QVKc/KH/Qhh5bBYJBuYwJ/Uhul5ExRm6XkdftUrzXrfgYl+K87tD9rtukWI9S473yuLk+CwAAR3Upe2PMYEkvWWvHHcVzt0qaYq2tOWj7zyWVWmvvOtJ7cCn7yNTuC2jKzxZq3rg8/fqS8U7HQTex1uoLDyzRih31evO/T1VOcpzTkRAGapo79PUnluvDzbW6atpA/ej8MVxu3UFN7T49vniH7nt/q6qbOjShIE03zS7UnNG5LASOsNDhD2hXfbtKd7eqsrFDNc0dqmkK3TZ3hm471NjmV2cvj1RNjHErNd6rlHivUkM/Wcmxyk6KVXbyfj+hx3FeznUAgMh1qEvZH/fIIWPMMEmbrbXWGDNJUoyk2tC+HGttlTFmoKTPS5pxvJ+H8PXOhmo1dfi5Slkf07U49Tid9bt39YtX1ul3l01wOhIctnzHbt3092Wqa+nUry8+UZcUMWPYaclxXn3l1EJdO3Ownl5aqr++u1k3PrpUI3KT9LXZhTr/xP6MjkCPa2j1aVN1szZXN2t7bYtKd7eFfroKoYPFe93KSo5RVlKsCjISNHFgulLjvUqO8yg5zqOk2K6f5DivEmPdivW45XEbxbhd8rpd8ri7Rga5XUZBa2WDUsBaBa1VMGgVtFKnP6h2f0CtnQG1dQbU7guozdf1uKndp8Y2vxrafHt/Gtt82lbboqXbd6uutVOf9m+oyXGevWVRv9Q49U+LV/+0eOWHbvunxSk5ztsLv+MAAHSfI5ZDxpgnJM2WlGWMKZV0mySvJFlr75F0kaQvGGN8ktokXbbfAtXPGGMyJfkk3Wyt3d39vwSEi1dXlSstwasZhZlOR0E3G5KVqBtPGao/v7VJl08p0LSh/DeORtZa/X3xDt3+4hrlpcbpma/N1Lj8VKdjYT9xXreunj5Il08p0MurynXXW5v1rX98ot++sUFfOWWoLikqYNQDjltdS6fWljdqbXmjNle3aHN1s7ZUN6umuXPvczwus7cwOWV4tgakJ2hAerzy0+PVLzVOWUmxSoztltUNeowvEFRdS6eqmzr2/TTvu1/V1K7i7btVsbJc/uCBLVJynEf5+xVG+enxGpSRoIKMBA3MTFAK5REAIMwc1bSy3sa0ssjT4Q+o6KdMKevL2joDmnPnO0qMdevlr58sL6MQokpbZ0Dff26V/rWsTLNHZuv3l01QWgKXqA93waDVv9dV6S9vb9LyHfXKSorV9bOG6OrpAxnZgCOy1mpHXavW7Ooqgkp2NaqkvFHlDe17n5OW4NWw7CQVZiepMCdRhdlJGpaTpPy0+KgZrRYIWlU3daisvk279vspq28P3bapoc13wGvSE7waGCqLBmUm7Hc/UXkpcXIzHRQA0EMONa2Mcgjd4q11VfriQx/rweum6LRROU7HQQ9ZUFKpGx4p1vfPGa0bThnqdBz0ksrGdn354WKt3tWgb5wxXF8/fTjr2EQYa60WbanTXW9v0nsba5Qc59G1MwbriycN5spy2KuhzaeVpfVavqNey3fs1oqd9drd2lVquF1Gw7KTNLpfssb0T9GYfqka1S9ZWRw/R6Wx3acdta3aWdeqHQf9lO1uO2DkkddtNCC9qzDa+5O57364j7hC3xcIWrX79k3TbPcF1e4LqMMfUKffKhC08gWDCgSs/MGuBeT9ARu6DR54G+yaCro/o31/xzj42gr7P3S7jDwuI7fbJa9r3yL0e7Z73K7Q7YH7Yj0uxXrcoduu+zGh+/z9BtGgx9YcAiTp1dXlSo71aOYwphv1ZXNG5+j0UTn6/cINOm98P/VLjXc6EnrYml0Nuv6hYjW2+/S3a4o0Z0yu05FwDIwxmlGYqRmFmVpV2qC73t6kv7y9Sfe9v0WXTxmoG04Zqvw0/n+OJv5AUBsqm7ViZ1cRtHxnvTZVNUvq+jI2PCdJZ47J04SBaRrXP1XDc5OYkngcUuK8Gpef+qlTcf2BoMob2g8sjWq7blfsrP+PUUdZSTEqyEjQgPQEFaTHqyAjQQXpCSrI6JrCxshe7C8QtGru6Lry354rAO6539LhV9Pe7T41dwRC+3xqDa3R1e4LhgqgriKowxfs9UXje5M3tK5ZrLerPIr5lAJpz+M4b9dVEGM9odv97sd593uOxx3av+f+vn17XhPrcXGVUTiOkUM4bv5AUFP+b6FOHZGt318+0ek46GE7alt15u/f0akjsvXXa/6jcEYf8ubaSt36xHKlxnt1/7VTNKZ/itOR0I02VTXrnnc267nlZZKkCybm66unFmpYTpLDydATfIGgVpU16KPNtVq0pVZLt+9Wa2dAkpSRGKOJBWmaODBNEwrSdWJBKmvihJGGVt9Bo41atKOuVTvruqav7T/qyGWkvJQ4DdivMCoIrfdUkJGgXKasRaRg0Kqpw6/GNp/qW/ctoF7f1rlvQfXQ9v33N7T51NzhP6rPiPe6lRTnUXKsR4mxHiXEdBUa8aGiIz7GHSo09m3bcz/W6wrduuXdb6F4ryt0e9DIHY/byONyhW6NXIcoRQ7+mmplD9gXsFaBQGiU0n6jkwIHjVYKBIPyBbpGNHUGgur0d/10+IPq8O8rvDp8ocd79+9/f99z99zfN2Jq3+2xMkaHLJ3iPF2/x7EHFUtxhyqj9hRO+702br8yKm7va9ycD46StVbWqs+MLGNaGXrMB5tqdNV9i3XP1ZM1b1ye03HQC+5+e7PueG0d/837KGutHvxgm372conG9k/V/dcWKSclzulY6CFl9W3627tb9OTHO9ThD+qsMXm6duZgTR+awb9iRjB/qAxatKVOH22pVfG2ur1l0MjcZE0dkqHJg9I1cWCaBmYk8N86QvkDQVU2dWhnXdeUtZ2721Ra16qdu1tVurtNFY3tB3zBdruMcpNj1S8tXnmpceqXEqd+aV2LhOelxql/aryyk2P5wthDOvwBNbT6VN+2r8ipb+3cd7+tUw1tftW3dnYVQftdRS94mK9sMW6XUhO8So3v+kkL3aaEbvddAbDryn977ifFeZQU41FirDtq1gjrSdbaruLI13WVxP2Lo3bfpxdKe+53+AJq33Mbev1/vk/X/oNfe/CC+J+Fx2UOKJf2TMHzuIzcLtfeYs/j2nO7/7S9Ax/vKf4kI2O6pgDu+aPFHLCta2PXY7Pfc0LbjNl7xcmgtbLWhq5E2fV7HAx2bQ+ECptgaF8wNEUxuN9zA/u9jy8QKhEPUSR27Qvut+/Ax1dOG6ifX3jC8RwiYYNyCD3mB8+t0jNLy7Tsh3MVH8OQ82jgCwR1/p/e1+7WTi389qksbNuH+ANB3f5iiR5dtF1njc3V7y6boIQYZiBHg9rmDj34wTY9tni76lt9GpGbpGtmDNbnJ+azxkkE8AeCWrOrUR9t6RoZ9PHWOrWEyqDhOUmaPrRrWuHUIRmsExRFOvwB7apvDxVHrSqvb9euhjZVNLSrvKFd5Q1tavcdONrB7TLKSopRVlLsfj+hx8kxykyM3Xs/IyEm6kqFDn9g7/Sspnb/vlE8e8qd1n33Dxjl0+pTmy9wyPd1GSktIWZvqbOn4EmN9yotwfuf2xO8Sovven6clylJ0cwfCO4rlvaMYjqoWOo4bEG1r4wKBIPyBW1ovah95cmeUmVfuRI8aFswNIrLhgppu7eYtgqNvFHXiK8992W1dyzYAftl5TZdI8qM6Rqt4zJ7frT31piuQmrPNhO6dbuMzEHP3VNsHa782rNG1ac+drt0Yn5qn1legXIIPSIYtJr2izc1ZXC67rpqstNx0ItW7KzXhXd9oGumD9JP5o9zOg66QVO7Tzc/vlzvbqjWV04dqu+eNarPDJ/F0Wv3BfTCJ7v0yEfbtLqsUUmxHl08eYCunj6IKWdhxB8IqqS8UYu21OqjzbX6eNvuvdNHhuUkafrQDE0fmqlpQzKVnUwZhE9nrVVDm0+76ttV0djWVRjVt6uysV21LZ2qae5QTVOHapo7D7nOTFKsRylxHqWEyouUOK9S4j1dZUZc17bkOE9oKtKBU5L2TYXZt80T+mJ3rL+eA77IBqw6Al1fevesm9N+wP3QWjqdAbX6Dix9mkLr8DS171ujp6nDr84jTB2K8biUFip00uJjQiVO6HFCzN6CZ8/+tISuoicpxsOfuQB6BQtSo0cs3bFb1U0dmjeun9NR0MsmFKTp2hmD9fBH2zR/Qr4mD0p3OhKOQ1VTu6574GNtqGzSLz9/gi6fOtDpSHBInNetS4sKdMnkAVq+s16PfrRdjy/eoYc+3KYZQzN1SdEAzRuXx4iyXhYIWpXs6iqDFm2p1ZKtdWoKlUFDsxM1f0L/rjJoaIZykpkGiqNjjFFaQozSEmIOu66ctV1r3uwpimqaO1Tb3KHalk41tvnV2L5v+lNZfZvWlvvU2N5VrBxbLh00SiB0PzSCwO0yXdNKDpr+cRyzayR1TbFJjvN0TbeK9So51qO8lLjQY8/eNXmS47xKCq3N01X67Ct6WLgdQKRi5BCOy09eLNFji7dr2Q/nKolpB1GnucOvM+98R4mxHr309VmK9fAXoki0vbZF19y/RDXNHbr76sk6dUS205EQZmqaO/SPj3fqHx/v1I66ViXGuHXuif108eQCTRmcznSGHhAIWq0t31cGLd5at/eL9tCsRE0LTRObPiSDNcEQtgJBq+b2rvJo/+kre0br7JsCs29ETyDYtcZIMLTgcDC0rkhgz3oiwa61Rtyma1qI5z/WROna7nZJbpdLMe6uNVXiY7oW542P2X/kknvviKaEGDdXjAIQFZhWhm5nrdVJv/y3xvRP0X3XTnE6Dhzy1voqffHBj3XzaYX6zlmjnI6Dz2h1WYOue3CJAkGrB784VRMK0pyOhDBmrdXH23br6aU79fLKcrV0BjQoM0EXTRqg+RP6a1BmotMRI9aeNYMWb+0aFbRka50aQ2XQkKzEA6aJ5aVSBgEAgGPDtDJ0u09KG7SroV3fPnOk01HgoNNG5uiiSQN0zztbdPa4fhqXn+p0JBylDzfV6MZHlyo13qtHrp+qwmzWk8HhGWM0dUiGpg7J0I8/N1avr6nQ00tLdeeCDbpzwQaNykvW2eP6ad64PI3ITeJf4A+jwx/QqtIGLd5ap8Vb67R0274FpIdkJerscf00o7Brmli/1HiH0wIAgL6OkUM4Zr94da3uf2+rlv5grlITuFpVNGto9Wnu795RRmKMXrhllmI80XXlkkj0yqpyffPJFRqSlaiHvzSVkQg4LmX1bXptdYVeX12hj7fXydqugmPeuDzNG5unE/JTo36h1YY2nz7ZWa+l23drydY6LduxWx2hhW1H5CZp2pCuK4lNY5oYAADoQUwrQ7ey1mr2b97WwIwEPXr9NKfjIAwsKKnUDY8U65tzhuubc0Y4HQeH8eii7frR86s1eWC67r92CuUuulV1U4feKKnQa6sr9NHmWvmDVukJXs0ozNTMwiydNCxLgzMT+vSookDQamNVk5bvqNfyHbu1bEe9NlU1S+paaHdMv5S9ZdDUIRnKSIxxODEAAIgWTCtDt1pb3qTtta366qmFTkdBmJg7JlfzJ/TXn/+9SXNG5zK9LAxZa/W7hRv1xzc3as7oHP35yklcVQXdLjs5VldNG6Srpg1SQ6tPb62v0vubavThphq9sqpCktQ/NU4zh2VpZmGmThyQpqFZiRE7sqjDH9CGimaVlDdoza5Glexq1Nryxr1TxNITvJo4MF3zx/fXpEHpOnFAqpLjKGQBAEB4YeQQjsmdb6zXn9/apI+/P0eZSbFOx0GYqG/t1Jm/e1ep8V69eOssiocwEgha/fD51Xp88Q5dWjRAP7/wBHncTP9D77HWamtNiz7YXKsPN9Xow821amjzSZISY9wa2z9VY/NTdEJ+qk7IT9XQ7CS5w6gwau30a2tNi7ZUt4Rum7WuokmbqprlD10/OynWo9H9kjW2f6rGF6RqYkG6BvXxUVIAACCyMHII3erV1RWaOiSDYggHSEuI0a8vGa9rH1iiX7++Xj88b4zTkSCp3RfQN59codfWVOim2YX6zlkj+bKKXmeM0dDsJA3NTtI10wcpELTaUNmkVWUNWl3WoFVlDXp88Y696/DEuF0akBGvgRkJGpSRoIKMBA3MSNDAzATlpcQpOc7bbeWRtVZ1LZ2qaGxXZWO7Kho6uu43tKu0vlVbqltU3tB+wGvy0+I1PDdJZ4zO0dj+qRrTL0UDMxIidgQUAACIbpRD+Mw2VTVrY1Wzrp4+1ukoCEOnjsjWF2YM0v3vb9UZo3I0c1iW05GiWmO7Tzc+UqxFW+r0o/PG6EuzhjgdCZAkuV1Go/ulaHS/FF1aVCCp63Lum6qbtaq0QRurmrWzrlXba1u1dNtuNXX4D3i9MVJyrEdpCTFKjfcqLcGrlHivYt0uyUguY2TUdesKDZJr6QiopcOvpg6/Wjr8ag7dNrb51RkI/sf7ZyXFqn9avGYMzdSQrMRQuZWowZmJio9hZCQAAOg7KIfwmb22ulySdNbYPIeTIFz979mj9f7GGv33Pz/Rq988RanxrK/hhKqmdl37wMfaWNmkP1w+QfMn5DsdCTgsj9ulUXkpGpWXcsB2a63qW33aUdeq7XWtqmnqUEObTw1tPtW3dqo+dL9sd5t8waCCoZ4naK2sDd2qa/paYqxHSbEe5aXEdd2P8yglzqvclFjlpcQpNzVOeSlxyk6OlZeplwAAIEpQDuEze3V1hSYNTOPS1zik+Bi37rxsgi66+0N9/9lV+tMVE5nG1Mu21bToCw8sUU1zh+6/bopOHZHtdCTgmBljlJ4Yo/TEGI0vSHM6DgAAQJ/DP4nhM9lZ16o1uxp19rh+TkdBmJtQkKZvzx2hl1aW66ninU7HiSqryxp08T0fqqndp8dvmE4xBAAAAOCwKIfwmbxRUilJOnNsrsNJEAm+dmqhZg3L0m0vrNGGyian40SFDzfV6PJ7FynW49bTX5upCYyyAAAAAHAElEP4TBaWVGpEbpIGZSY6HQURwOUyuvOy8UqK9eiWx5ep3RdwOlKf9vLKcl334MfKT4vXM1+bqcLsJKcjAQAAAIgAlEM4avWtnVqyrU5zxzBqCEcvJzlOd146QRsqm3X7i2ucjtNnPfLRNt3yxDKdOCBVT31lBmuCAQAAADhqlEM4am+vr1YgaDVnNOUQPptTRmTrptmFemLJTtYf6mbWWv3qtXX60fNrdMaoXD325WlKTeDqcAAAAACOHuUQjtqCkkrlJMdq/IA0p6MgAn177gjNGpalHzy3WqtKG5yO0yf4AkH91z8/0V1vb9YVUwfqnqsnKc7rdjoWAAAAgAhDOYSj0uEP6O31VTpjdK5cLi5Jjs/O43bpj1dMVHZSrL762FLVtXQ6HSmitXT4df3DxfrXsjJ9a84I/fzCcfK4OaUDAAAA+Oz4JoGjsmhLnVo6A5o7JsfpKIhgGYkxuvvqSapu7tCtTyyTPxB0OlJEqm7q0OX3LtIHm2r0y8+foG/MGS5jKG0BAAAAHBvKIRyVBSUVive6NbMwy+koiHAnDkjTzy4Ypw821epnL691Ok7E2VbToovv+VAbq5p07zWTdfnUgU5HAgAAABDhPE4HQPiz1mphSZVOGZHFeiboFpcWFWhDRZPue3+rhmQl6tqZg52OFBGWbt+tGx8pVtBaPX7DdE0amO50JAAAAAB9ACOHcESryxpV0diuuWPynI6CPuR/zxmtOaNzdfuLa/TWuiqn44S951eU6Yq/LVJirEdPf20mxRAAAACAbkM5hCNaUFIhl5FOH8V6Q+g+bpfRHy6foNH9UnTL48u0trzR6UhhKRi0uvON9frGkys0oSBNz918kgqzk5yOBQAAAKAPoRzCES1YW6WiQRnKSIxxOgr6mMRYj+6/doqS47y69oEl2lHb6nSksNLWGdCtTyzXH/+9SZdMHqDHrp/G/4cAAAAAuh3lEA5rZ12r1pY3au6YXKejoI/KS43TI9dPVWcgqKvuX6TKxnanI4WFqsZ2XX7vR3pldbn+9+xR+tXFJyrGwykbAAAAQPfjmwYO6821lZKkOZRD6EEjcpP18Benqq65U9fcv1i7WzqdjuSo1WUNmv+XD7Sxqln3XlOkr5xayKXqAQAAAPSYI5ZDxpgHjDFVxpjVh9g/3xiz0hizwhhTbIyZtd++bxlj1hhjVhtjnjDGxHVnePS8BWsrNSwnSUOyEp2Ogj5ufEGa7rt2irbVtuq6B5eooc3ndCRHPL+iTJfc85GMpKe/OpNRewAAAAB63NGMHHpI0rzD7H9T0nhr7QRJX5J0nyQZY/IlfV1SkbV2nCS3pMuPJyx6V0ObT4u31GnOaL6confMKMzU3VdNUkl5o6782yLVRdEIok5/ULc9v1rfeHKFThiQquduOUlj+qc4HQsAAABAFDhiOWStfVdS3WH2N1trbehhoiS7326PpHhjjEdSgqRdx5EVvezt9VXyBy0jF9Crzhidq799oUibqpp1+b0fqaqp769BVN7Qpsvv/UgPf7RdX541RH//8jTlJDPQEgAAAEDv6JY1h4wxFxpj1kl6WV2jh2StLZP0G0k7JJVLarDWvtEdn4fesaCkUllJMZpYkOZ0FESZ2SNz9OB1U1S6u02X/XWRyurbnI7UY95YU6Gz//Ce1lU06c9XTtQPzhsjr5vl4AAAAAD0nm75BmKtfdZaO0rSBZJ+KknGmHRJ8yUNkdRfUqIx5upDvYcx5sbQmkXF1dXV3RELx6HTH9Q766t1xqhcuVwshIveN3NYlh69fqpqmjs0/88f6JOd9U5H6lbtvoB+9Pxq3fjoUuWnxeulW2fpvBP7Ox0LAAAAQBTq1n+eDk1BKzTGZEmaI2mrtbbaWuuT9C9JMw/z2nuttUXW2qLs7OzujIVjsHhrrZo6/Ewpg6MmD8rQM1+bqTivS5fd+5FeXVXudKRusaq0QZ/78/t6JDSN7F83zdTQ7CSnYwEAAACIUsddDhljhpnQNZaNMZMkxUiqVdd0sunGmITQ/jMkrT3ez0PvWFBSqTivSycNy3I6CqLciNxkPXfzSRrTL0Vf+/sy/enNjQoG7ZFfGIY6/UH99o31uuCuD9TQ5tODX5yiH5w3RrEet9PRAAAAAEQxz5GeYIx5QtJsSVnGmFJJt0nySpK19h5JF0n6gjHGJ6lN0mWhBaoXG2OelrRMkl/Sckn39sQvAt3LWquFJZU6eXi24mP40grnZSXF6vEbpuu7z6zUbxds0JJtdbrz0gnKTo51OtpRW7p9t77/7Cqtq2jSRZMG6EfnjVFqgtfpWAAAAAAgs+9CY+GjqKjIFhcXOx0jaq0ua9B5f3pfv7roRF06pcDpOMBe1lo9vmSHfvJiiZLjvPr9ZRM0a3h4j26rbe7QHa+t01PFpcpLidP/XThOZ4xmuiYAAACA3meMWWqtLTp4+xFHDiH6LFxbKWOk00fnOB0FOIAxRldNG6TJg9J1y+PLdfX9i3XF1IH63rxRYTcKp8Mf0BOLd+h3CzeqpcOvr5wyVF8/Y7gSYzntAgAAAAgvfEvBf1hQUqnJA9OVlRQ5U3YQXUblpejFW2bpN2+s14MfbNWCkgr94Nwxmj+hv0JLoDkmELR6fkWZ7lywQaW72zSzMFM//txYjchNdjQXAAAAABwK5RAOUFbfpjW7GvW9s0c5HQU4rPgYt3543hhdODFf339utb75jxV65KNt+tbcEZo1LKvXS6IOf0DPL9+lv723RRurmjUuP0U/v/AEnTy897MAAAAAwGdBOYQDvLm2UpK4hD0ixrj8VP3razP1VPFO/fHNjbrm/iUqGpSuW88YrpOHZcnl6tliprKxXU99vFMPf7RdNc0dGpWXrD9fOVHnjOvX458NAAAAAN2BcggHWFBSqaFZiSrMTnI6CnDU3C6jK6YO1Ocn5eup4lLd9dYmXfvAEg1Ij9clkwt0cdEA5afFd9vnNXf49caaCj27vEwfbKpR0EqnjsjWDScP1UnDMhkpBAAAACCiUA5hr8Z2nxZtqdWXThridBTgmMR63Lpm+iBdWjRAr62u0FPFO/W7hRv0+zc3aFz/VM0anqWTh2Vp0qB0xXndR/2+zR1+rStv1AebavX+pmot31Evf9CqICNet5w2TBdMzNdQClUAAAAAEYpyCHu9s75avoDVHKaUIcLFetyaPyFf8yfka2ddq55bXqZ3N1brb+9u0d1vb5YxUn5avIZmJ2lwZoKSYj1KiHErzutWhz+oxnafGtt8qmho14bKZpXVt0mSjJHG9U/VDacM1RmjcjR5UDqjhAAAAABEPMoh7LWgpFKZiTGaNDDd6ShAtynISNCtZwzXrWcMV3OHX4s212r1rgZtqW7RlppmrdixWy2dAQWCdu9rYtwupcR7lZ0cq8mD0nXltIEanpOkKYMzlJ4Y4+CvBgAAAAC6H+UQJEm+QFBvra/SvLF5crOILvqopFiP5ozJ/Y/RcdZa+QJWbZ0BxXpdn2nKGQAAAABEOsohSJKWbK1TU7ufKWWISsYYxXiMYjwup6MAAAAAQK/jmxAkdU0pi/W4dPLwLKejAAAAAACAXkQ5BFlrtaCkUrOGZSkhhsFkAAAAAABEE8ohaG15k8rq2zSXKWUAAAAAAEQdyiFo4dpKGSOdMZpyCAAAAACAaEM5BC0oqdSEgjRlJ8c6HQUAAAAAAPQyyqEoV97QplVlDUwpAwAAAAAgSlEORbmFa6skSXOZUgYAAAAAQFSiHIpyC0sqNTgzQcNykpyOAgAAAAAAHEA5FMWaO/z6aHOt5ozOlTHG6TgAAAAAAMABlENR7J311eoMBFlvCAAAAACAKEY5FMUWrq1UeoJXkwelOx0FAAAAAAA4hHIoSvkCQf17XZVOG5Ujj5vDAAAAAACAaEUrEKWKt+1WQ5tPZzKlDAAAAACAqEY5FKUWlFQqxuPSycOznY4CAAAAAAAcRDkUhay1WrC2QicVZiox1uN0HAAAAAAA4CDKoSi0obJZO+vaNIcpZQAAAAAARD3KoSi0oKRCkjRnNOUQAAAAAADRjnIoCi1YW6XxBWnKTYlzOgoAAAAAAHAY5VCUqWxs1yc76zV3dI7TUQAAAAAAQBigHIoyC9dWSpLmjslzOAkAAAAAAAgHlENRZmFJpQoy4jUiN8npKAAAAAAAIAxQDkWRlg6/Pthcq7mj82SMcToOAAAAAAAIA5RDUeSdDdXq9Ac1l0vYAwAAAACAEMqhKPL6mgqlJ3g1ZXC601EAAAAAAECYOGI5ZIx5wBhTZYxZfYj9840xK40xK4wxxcaYWaHtI0Pb9vw0GmO+2c35cZQ6/UH9e12V5ozOlcdNJwgAAAAAALocTUvwkKR5h9n/pqTx1toJkr4k6T5Jstaut9ZOCG2fLKlV0rPHExbH7qMttWpq9+ussVylDAAAAAAA7HPEcsha+66kusPsb7bW2tDDREn2U552hqTN1trtx5QSx+31NRVKiHFr1vAsp6MAAAAAAIAw0i3zi4wxFxpj1kl6WV2jhw52uaQnjvAeN4ampRVXV1d3RyyEBINWC0oqNXtktuK8bqfjAAAAAACAMNIt5ZC19llr7ShJF0j66f77jDExkj4n6Z9HeI97rbVF1tqi7Ozs7oiFkOU7d6u6qYMpZQAAAAAA4D9068rEoSlohcaY/ecunS1pmbW2sjs/C0fv9TWV8rqNThuV43QUAAAAAAAQZo67HDLGDDPGmND9SZJiJNXu95QrdIQpZeg51lq9vqZCMwqzlBLndToOAAAAAAAIM54jPcEY84Sk2ZKyjDGlkm6T5JUka+09ki6S9AVjjE9Sm6TL9ixQbYxJkDRX0ld6JD2OaH1lk7bXturGU4Y6HQUAAAAAAIShI5ZD1torjrD/Dkl3HGJfq6TMY4uG7vD66koZI80dk+t0FAAAAAAAEIa6dc0hhJ83Sio0aWC6cpLjnI4CAAAAAADCEOVQH7azrlVrdjXqrLGMGgIAAAAAAJ+OcqgPe6Ok6wJxXMIeAAAAAAAcCuVQH/b6mgqNykvWoMxEp6MAAAAAAIAwRTnUR9U2d6h4W53OZNQQAAAAAAA4DMqhPmrh2koFrVhvCAAAAAAAHBblUB/1+ppKDUiP15h+KU5HAQAAAAAAYYxyqA9q7vDr/Y01OmtsnowxTscBAAAAAABhjHKoD3p7fZU6A0GdOYYpZQAAAAAA4PAoh/qg19dUKjMxRkWDM5yOAgAAAAAAwhzlUB/T7gvorXVVmjM6V24XU8oAAAAAAMDhUQ71Me9trFFzh1/nnNjP6SgAAAAAACACUA71Ma+sKldaglczCzOdjgIAAAAAACIA5VAf0u4LaGFJpc4ckyuvm/+0AAAAAADgyGgQ+pD3NtaoqcOvc05gShkAAAAAADg6lEN9yCurypUa79VJw7KcjgIAAAAAACIE5VAf0eHvmlJ21limlAEAAAAAgKNHi9BHvLeBKWUAAAAAAOCzoxzqI5hSBgAAAAAAjgXlUB/Q4Q9oAVcpAwAAAAAAx4AmoQ/YO6XsRKaUAQAAAACAz4ZyqA944ZNdSk/wahZTygAAAAAAwGdEORThWjv9WlBSqXNO6MeUMgAAAAAA8JnRJkS4BSWVavMF9Lnx/Z2OAgAAAAAAIhDlUIR7YcUu9UuN05TBGU5HAQAAAAAAEYhyKILtbunUOxuqdf74/nK5jNNxAAAAAABABKIcimCvrq6QP2iZUgYAAAAAAI4Z5VAEe+GTMg3NStTY/ilORwEAAAAAABGKcihCVTS0a/HWOp0/vr+MYUoZAAAAAAA4NpRDEeqllbtkrfS5CUwpAwAAAAAAx45yKEI9u7xM4/JTVJid5HQUAAAAAAAQwSiHItC6ikat2dWoiyYNcDoKAAAAAACIcJRDEeiZpaXyuAxXKQMAAAAAAMftiOWQMeYBY0yVMWb1IfbPN8asNMasMMYUG2Nm7bcvzRjztDFmnTFmrTFmRneGj0b+QFDPLt+l00blKDMp1uk4AAAAAAAgwh3NyKGHJM07zP43JY231k6Q9CVJ9+237w+SXrPWjpI0XtLaY4uJPd7bWKOa5g6mlAEAAAAAgG7hOdITrLXvGmMGH2Z/834PEyVZSTLGpEg6RdJ1oed1Suo8jqyQ9PTSUqUneHX6qBynowAAAAAAgD6gW9YcMsZcaIxZJ+lldY0ekqShkqolPWiMWW6Muc8Yk9gdnxetGlp9WlBSqfkT8hXjYbkoAAAAAABw/LqlYbDWPhuaOnaBpJ+GNnskTZJ0t7V2oqQWSd871HsYY24MrVlUXF1d3R2x+pwXV+5SZyDIlDIAAAAAANBtunX4ibX2XUmFxpgsSaWSSq21i0O7n1ZXWXSo195rrS2y1hZlZ2d3Z6w+45llpRqRm6Rx+SlORwEAAAAAAH3EcZdDxphhxhgTuj9JUoykWmtthaSdxpiRoaeeIankeD8vWm2qatLyHfW6aNIAhX67AQAAAAAAjtsRF6Q2xjwhabakLGNMqaTbJHklyVp7j6SLJH3BGOOT1CbpMmutDb38Vkl/N8bESNoi6Yvd/iuIEo8v3imv2+iiyUwpAwAAAAAA3edorlZ2xRH23yHpjkPsWyGp6JiSYa92X0DPLCvVmWPzlJUU63QcAAAAAADQh3DJqwjwyqpyNbT5dNXUgU5HAQAAAAAAfQzlUAR4fPEODclK1IzCTKejAAAAAACAPoZyKMxtqGxS8fbdumJqAQtRAwAAAACAbkc5FOYeX7xDMW6XLp5c4HQUAAAAAADQB1EOhbF2X0D/WlaqeePylJEY43QcAAAAAADQB1EOhbGXVparsd2vK6exEDUAAAAAAOgZlENhylqrB97fquE5SZo2JMPpOAAAAAAAoI+iHApTi7bUqaS8UV8+eQgLUQMAAAAAgB5DORSm7n9/izITYzR/Qr7TUQAAAAAAQB9GORSGtlQ3a+HaKl01fZDivG6n4wAAAAAAgD6McigMPfjBNsW4Xbpm+iCnowAAAAAAgD6OcijM1Ld26p9Ld2r+hP7KTo51Og4AAAAAAOjjKIfCzN8X71C7L6jrTx7idBQAAAAAABAFKIfCSLsvoIc/3KZZw7I0Ki/F6TgAAAAAACAKUA6FkSeW7FBVU4duPm2Y01EAAAAAAECUoBwKE+2+gO5+e7OmDcnQjMJMp+MAAAAAAIAoQTkUJh5f3DVq6JtzRjgdBQAAAAAARBHKoTDQ7gvo7nc2a/pQRg0BAAAAAIDeRTkUBh5fvEPVTR36xhmMGgIAAAAAAL2LcshhjBoCAAAAAABOohxy2H3vbVF1U4e+xVpDAAAAAADAAZRDDtpV36a/vLVZ88bmadpQRg0BAAAAAIDeRznkoF+8uk5Ba/X9c0c7HQUAAAAAAEQpyiGHLNlapxc/2aWvnFqogowEp+MAAAAAAIAoRTnkgEDQ6rYX1ig/LV5fO7XQ6TgAAAAAACCKUQ454PHF27W2vFHfP3e04mPcTscBAAAAAABRjHKol22tadEvXl2nk4Zl6uxxeU7HAQAAAAAAUY5yqBd1+oP6xpPL5XW79JtLxssY43QkAAAAAAAQ5TxOB4gmv1u4QStLG3TP1ZPULzXe6TgAAAAAAACMHOotH26u0T3vbNYVUws0b1w/p+MAAAAAAABIohzqFbvq2/Stf6zQkKxE/fC8MU7HAQAAAAAA2ItyqIdVN3Xo6vsWq7UjoL9cOUkJMczkAwAAAAAA4YOmogfVt3bqmvsXq7yhXY99eapG90txOhIAAAAAAMABjjhyyBjzgDGmyhiz+hD75xtjVhpjVhhjio0xs/bbt80Ys2rPvu4MHu6a2n269sGPtaWmRfddW6TJgzKcjgQAAAAAAPAfjmZa2UOS5h1m/5uSxltrJ0j6kqT7Dtp/mrV2grW26JgSRqi/vbtFa8oadNeVk3TSsCyn4wAAAAAAAHyqI04rs9a+a4wZfJj9zfs9TJRkuyFXxLv1jOE6eUS2pgxmxBAAAAAAAAhf3bIgtTHmQmPMOkkvq2v00B5W0hvGmKXGmBu747MihdftohgCAAAAAABhr1vKIWvts9baUZIukPTT/XadZK2dJOlsSTcbY0451HsYY24MrVlUXF1d3R2xAAAAAAAAcATdeil7a+27kgqNMVmhx7tCt1WSnpU09TCvvddaW2StLcrOzu7OWAAAAAAAADiE4y6HjDHDjDEmdH+SpBhJtcaYRGNMcmh7oqQzJX3qFc8AAAAAAADgjCMuSG2MeULSbElZxphSSbdJ8kqStfYeSRdJ+oIxxiepTdJl1lprjMmV9GyoN/JIetxa+1qP/CoAAAAAAABwTI7mamVXHGH/HZLu+JTtWySNP/ZoAAAAAAAA6GnduuYQAAAAAAAAIgvlEAAAAAAAQBSjHAIAAAAAAIhilEMAAAAAAABRjHIIAAAAAAAgihlrrdMZ/oMxplrSdqdzdIMsSTVOhwAOg2MU4Y5jFJGA4xThjmMU4Y5jFOGuLx2jg6y12QdvDMtyqK8wxhRba4uczgEcCscowh3HKCIBxynCHccowh3HKMJdNByjTCsDAAAAAACIYpRDAAAAAAAAUYxyqGfd63QA4Ag4RhHuOEYRCThOEe44RhHuOEYR7vr8McqaQwAAAAAAAFGMkUMAAAAAAABRjHKohxhj5hlj1htjNhljvud0HkCSjDHbjDGrjDErjDHFoW0ZxpgFxpiNodt0p3MiehhjHjDGVBljVu+37ZDHpDHmf0Pn1fXGmLOcSY1ocohj9MfGmLLQuXSFMeac/fZxjKJXGWMKjDFvGWPWGmPWGGO+EdrOuRRh4TDHKOdShAVjTJwxZokx5pPQMXp7aHtUnUeZVtYDjDFuSRskzZVUKuljSVdYa0scDYaoZ4zZJqnIWluz37ZfSaqz1v4yVGSmW2u/61RGRBdjzCmSmiU9Yq0dF9r2qcekMWaMpCckTZXUX9JCSSOstQGH4iMKHOIY/bGkZmvtbw56Lscoep0xpp+kftbaZcaYZElLJV0g6TpxLkUYOMwxeqk4lyIMGGOMpERrbbMxxivpfUnfkPR5RdF5lJFDPWOqpE3W2i3W2k5JT0qa73Am4FDmS3o4dP9hdf1hDfQKa+27kuoO2nyoY3K+pCettR3W2q2SNqnrfAv0mEMco4fCMYpeZ60tt9YuC91vkrRWUr44lyJMHOYYPRSOUfQq26U59NAb+rGKsvMo5VDPyJe0c7/HpTr8CRDoLVbSG8aYpcaYG0Pbcq215VLXH96SchxLB3Q51DHJuRXh5BZjzMrQtLM9w8w5RuEoY8xgSRMlLRbnUoShg45RiXMpwoQxxm2MWSGpStICa23UnUcph3qG+ZRtzN9DODjJWjtJ0tmSbg5NlwAiBedWhIu7JRVKmiCpXNJvQ9s5RuEYY0ySpGckfdNa23i4p37KNo5T9LhPOUY5lyJsWGsD1toJkgZImmqMGXeYp/fJY5RyqGeUSirY7/EASbscygLsZa3dFbqtkvSsuoY/Vobmgu+ZE17lXEJA0qGPSc6tCAvW2srQXyKDkv6mfUPJOUbhiNAaGc9I+ru19l+hzZxLETY+7RjlXIpwZK2tl/S2pHmKsvMo5VDP+FjScGPMEGNMjKTLJb3gcCZEOWNMYmgRQBljEiWdKWm1uo7Na0NPu1bS884kBPY61DH5gqTLjTGxxpghkoZLWuJAPkS5PX9RDLlQXedSiWMUDggtpHq/pLXW2jv328W5FGHhUMco51KEC2NMtjEmLXQ/XtIcSesUZedRj9MB+iJrrd8Yc4uk1yW5JT1grV3jcCwgV9KzXX8+yyPpcWvta8aYjyU9ZYy5XtIOSZc4mBFRxhjzhKTZkrKMMaWSbpP0S33KMWmtXWOMeUpSiSS/pJsj/aoQCH+HOEZnG2MmqGsI+TZJX5E4RuGYkyRdI2lVaL0MSfp/4lyK8HGoY/QKzqUIE/0kPRy66rhL0lPW2peMMR8pis6jXMoeAAAAAAAgijGtDAAAAAAAIIpRDgEAAAAAAEQxyiEAAAAAAIAoRjkEAAAAAAAQxSiHAAAAAAAAohjlEAAAAAAAQBSjHAIAAAAAAIhilEMAAAAAAABR7P8DLXe9yCRerMgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1440x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(0, numNuclearPoints, numNuclearPoints)\n",
    "\n",
    "plt.rcParams[\"figure.figsize\"] = (20,5)\n",
    "plt.plot(x, distances)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.3879589486243356\n"
     ]
    }
   ],
   "source": [
    "print(X) # the experimental value given in the book is 1.401, the equilibrium with Hartree-Fock using only s-type orbitals is 1.3881"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}