notebooks/nanoDFT-demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7cb2dc26",
   "metadata": {},
   "source": [
    "Copyright (c) 2023 Graphcore Ltd. All rights reserved.\n",
    "\n",
    "# NanoDFT notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35622c3d",
   "metadata": {},
   "source": [
    "## Dependencies and configuration\n",
    "\n",
    "Install the JAX experimental for IPU (and addons)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "6dcb855a",
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install -e \"..[ipu]\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "aba6a829",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "os.environ['OMP_NUM_THREADS'] = \"16\"\n",
    "\n",
    "if \"JAX_IPU_USE_MODEL\" in os.environ:\n",
    "    # executable caching is not supported with IPU model\n",
    "    os.environ[\"TF_POPLAR_FLAGS\"] = \"--show_progress_bar=true\"\n",
    "else:\n",
    "    os.environ[\"TF_POPLAR_FLAGS\"] = \"\"\"\n",
    "      --executable_cache_path=/tmp/ipu-ef-cache\n",
    "      --show_progress_bar=true\n",
    "    \"\"\"\n",
    "\n",
    "import numpy as np\n",
    "import jax\n",
    "from tqdm import tqdm\n",
    "\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import py3Dmol\n",
    "\n",
    "from pyscf_ipu.nanoDFT import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ca1b593b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Platform=ipu\n",
      "Number of devices=4\n",
      "IpuDevice(id=0, num_tiles=1472, version=ipu2)\n",
      "IpuDevice(id=1, num_tiles=1472, version=ipu2)\n",
      "IpuDevice(id=2, num_tiles=1472, version=ipu2)\n",
      "IpuDevice(id=3, num_tiles=1472, version=ipu2)\n"
     ]
    }
   ],
   "source": [
    "# Check IPU hardware configuration\n",
    "assert jax.default_backend() == \"ipu\"\n",
    "print(f\"Platform={jax.default_backend()}\")\n",
    "print(f\"Number of devices={jax.device_count()}\")\n",
    "\n",
    "devices = jax.devices()\n",
    "print(\"\\n\".join([str(d) for d in devices]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c7475f5",
   "metadata": {},
   "source": [
    "## Dissociation of Molecular Hydrogen"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8313857c",
   "metadata": {},
   "source": [
    "Density functional theory (DFT) is a widely used approximation for quantum-mechanics.\n",
    "We have created a simplified application `nanoDFT` to facilitate exploring different\n",
    "accuracy - compute tradeoffs in DFT simulations. \n",
    "\n",
    "Below we use `nanoDFT` to calculate the dissociation curve for the hydrogen molecule.\n",
    "`nanoDFT` allows us to simulate varying the distance between a pair of hydrogen atoms \n",
    "while keeping all other parameters the same.  The energy calculated by DFT will\n",
    "tend to:\n",
    "\n",
    "* Reach a minimum which corresponds to the experimentally measurable bond length\n",
    "* Increases very slowly as the interatomic distance is increased \n",
    "* Rapidly increases as the two hydrogen atoms are brought very close together.\n",
    "\n",
    "This nicely captures the simple intuition that atoms tend to attract one another\n",
    "to form chemical bonds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "5b60342b",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 16/16 [00:52<00:00,  3.26s/it]\n"
     ]
    }
   ],
   "source": [
    "opts, _ = nanoDFT_options(backend=\"ipu\", float32=True)\n",
    "\n",
    "def h2_geom(r):\n",
    "    return [\n",
    "        [\"H\", (0.0, 0.0, 0.0)], \n",
    "        [\"H\", (  r, 0.0, 0.0)]\n",
    "    ]\n",
    "    \n",
    "\n",
    "def h2_energy(r: float) -> float:\n",
    "    mol = build_mol(h2_geom(r), opts.basis)\n",
    "    energies,_dbginfo = nanoDFT(mol, opts)\n",
    "    Etotal, E_core, E_J2, negE_K4, E_xc, E_nuc = energies[-1, :]\n",
    "    return Etotal\n",
    "\n",
    "\n",
    "r = np.linspace(0.2, 2.5, 16)\n",
    "energy = np.array([h2_energy(v) for v in tqdm(r)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "3d9a786a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj4AAAHHCAYAAAC/R1LgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABXlUlEQVR4nO3deXhTVf4G8Ddps3RNW7rvC1Ioa6mAIKsgRQRFQcUZEQQBHRwVEQRxWH8Oq6CjjKCyOTiCC6KDCxRlp7IIVdZigbZQukC3dN9yfn+UBEIX0tLmpsn7eZ77QG7Ovfmmkeb13HPOlQkhBIiIiIhsgFzqAoiIiIjMhcGHiIiIbAaDDxEREdkMBh8iIiKyGQw+REREZDMYfIiIiMhmMPgQERGRzWDwISIiIpvB4ENEREQ2g8GHiEySnJwMmUyGDRs2SF1Kvfr374/+/fs36TnnzZsHmUzWpOckImkw+BBZgA0bNkAmk+HYsWO1Pt+/f3906NDB8Li4uBirVq3C4MGD4efnBxcXF0RHR+PDDz9EVVWVSa8pk8kMm729PTw8PBATE4NXXnkFZ86caZL31ZIUFxdj3rx52LNnj9Sl1FBaWoqVK1eiR48e0Gg0UKvVaNOmDV566SWcP39e6vKIWhQZ79VFJL0NGzbgueeew9GjR3HvvffWeL5///64fv06Tp06BQA4deoUOnXqhIEDB2Lw4MFwdXXFjh078M033+DZZ5/Fxo0b7/iaMpkMDz74IJ599lkIIZCfn4/ff/8dX375JYqKirBkyRK89tprhvZCCJSVlUGhUMDOzq7p3nwTKy8vBwAolcoGHXf9+nV4eXlh7ty5mDdvntFzlZWVqKyshFqtbqoyG1TXkCFD8Ntvv2HYsGEYNGgQnJ2dkZiYiM2bNyMjI8PwnonozuylLoCIGs7X1xcnT55E+/btDfsmT56M8ePHY/369fjHP/6B1q1b3/E8bdq0wTPPPGO0b/HixRg+fDimTZuGtm3bYujQoQCqg5IUX/wN1dDAYwp7e3vY20vz63LcuHE4ceIEvvrqK4wcOdLouYULF2L27NlN8jqVlZXQ6XTN8vMjsiS81EXUAnl6ehqFHr3HHnsMAHD27NlGn7tVq1bYvHkz7O3t8fbbbxv21zbGJyMjA8899xwCAwOhUqng5+eHRx99FMnJyYY2x44dQ2xsLDw9PeHg4ICwsDCMHz/e6DWLioowbdo0BAUFQaVSITIyEsuXL0dtHdKbNm1C9+7d4ejoCHd3d/Tt2xc7d+40PH/7GJ/y8nLMmTMHMTEx0Gg0cHJyQp8+fbB7926j9+bl5QUAmD9/vuESoL7np7YxPpWVlVi4cCEiIiKgUqkQGhqKN998E2VlZUbtQkNDMWzYMBw4cADdu3eHWq1GeHg4Pv300/o/CACHDx/G999/jwkTJtQIPQCgUqmwfPnyOt+73rhx4xAaGmr0fmUyGZYvX453333X8B5OnDgBe3t7zJ8/v8Y5EhMTIZPJ8MEHHxj25eXl4dVXXzV8bq1bt8aSJUug0+nu+N6IpMIeHyILkp+fj+vXr9fYX1FRYdLxGRkZAKqD0d0IDg5Gv379sHv3bmi1Wri6utbabuTIkTh9+jT+/ve/IzQ0FFlZWYiLi0Nqaqrh8eDBg+Hl5YWZM2fCzc0NycnJ2Lp1q+EcQgg88sgj2L17NyZMmIAuXbpgx44dmD59OtLS0rBy5UpD2/nz52PevHno1asXFixYAKVSicOHD+OXX37B4MGDa61Rq9Xik08+wdNPP42JEyeioKAAa9euRWxsLI4cOYIuXbrAy8sLH374IV588UU89thjePzxxwEAnTp1qvNn9Pzzz2Pjxo0YNWoUpk2bhsOHD2PRokU4e/YsvvnmG6O2SUlJGDVqFCZMmICxY8di3bp1GDduHGJiYmoNsHrfffcdAGDMmDF1trkb69evR2lpKSZNmmQIrv369cMXX3yBuXPnGrXdsmUL7Ozs8MQTTwCoHhPVr18/pKWlYfLkyQgODsahQ4cwa9YspKen4913322WmonumiAiya1fv14AqHdr3759vecoKysTUVFRIiwsTFRUVNzxNQGIKVOm1Pn8K6+8IgCI33//XQghxKVLlwQAsX79eiGEELm5uQKAWLZsWZ3n+OabbwQAcfTo0TrbbNu2TQAQ//d//2e0f9SoUUImk4mkpCQhhBB//vmnkMvl4rHHHhNVVVVGbXU6neHv/fr1E/369TM8rqysFGVlZUbtc3NzhY+Pjxg/frxh37Vr1wQAMXfu3Bo1zp07V9z66zIhIUEAEM8//7xRu9dff10AEL/88othX0hIiAAg9u3bZ9iXlZUlVCqVmDZtWl0/FiGEEI899pgAIHJzc+ttp3f7e9cbO3asCAkJMTzWf5aurq4iKyvLqO2aNWsEAHHy5Emj/VFRUeKBBx4wPF64cKFwcnIS58+fN2o3c+ZMYWdnJ1JTU02qmcjceKmLyIKsWrUKcXFxNbb6eh70XnrpJZw5cwYffPBBk4xHcXZ2BgAUFBTU+ryDgwOUSiX27NmD3NzcWtu4ubkBALZv315nr9UPP/wAOzs7vPzyy0b7p02bBiEEfvzxRwDAtm3boNPpMGfOHMjlxr+66ptqbmdnZxi3otPpkJOTg8rKStx77704fvx4ncfV54cffgAAo8Hf+poB4PvvvzfaHxUVhT59+hgee3l5ITIyEhcvXqz3dbRaLQDAxcWlUXXeyciRIw2X+PQef/xx2NvbY8uWLYZ9p06dwpkzZ/DUU08Z9n355Zfo06cP3N3dcf36dcM2aNAgVFVVYd++fc1SM9Hd4qUuIgvSvXv3Wmd16b9c6rJs2TJ8/PHHWLhwoWEw8t0qLCwEUPeXrkqlwpIlSzBt2jT4+Pjgvvvuw7Bhw/Dss8/C19cXANCvXz+MHDkS8+fPx8qVK9G/f3+MGDECf/nLX6BSqQAAKSkp8Pf3r/E67dq1MzwPABcuXIBcLkdUVFSD38vGjRvxzjvv4Ny5c0YBLCwsrMHn0tckl8trDCD39fWFm5uboWa94ODgGudwd3evMzDq6S8xFhQUGEJkU6rt/Xt6emLgwIH44osvsHDhQgDVl7ns7e0NlwAB4M8//8Qff/xRIzjpZWVlNXm9RE2BPT5ELdyGDRvwxhtv4IUXXsBbb73VZOc9deoU7Ozs6g0Hr776Ks6fP49FixZBrVbjH//4B9q1a4cTJ04AqO6J+eqrrxAfH4+XXnoJaWlpGD9+PGJiYgzBqrlt2rQJ48aNQ0REBNauXYuffvoJcXFxeOCBB+56EK6pixrWNf1f3GE1kbZt2wIATp48eVf11LW2k4ODQ637R48ejfPnzyMhIQEA8MUXX2DgwIFGY8d0Oh0efPDBWnso4+Liah2MTWQJGHyIWrBvv/0Wzz//PB5//HGsWrWqyc6bmpqKvXv3omfPnne8zBIREYFp06Zh586dOHXqFMrLy/HOO+8Ytbnvvvvw9ttv49ixY/jss89w+vRpbN68GQAQEhKCq1ev1rikdu7cOcPz+tfR6XQNXlzxq6++Qnh4OLZu3YoxY8YgNjYWgwYNQmlpqVG7hqzMHBISAp1Ohz///NNof2ZmJvLy8gw1363hw4cDqA5vpnB3d0deXl6N/bf3QN3JiBEjoFQqsWXLFiQkJOD8+fMYPXq0UZuIiAgUFhZi0KBBtW619XIRWQIGH6IWat++fRg9ejT69u2Lzz77rMa4l8bKycnB008/jaqqqnrXiCkuLq4RHiIiIuDi4mKY0p2bm1ujV6NLly4AYGgzdOhQVFVVGU2TBoCVK1dCJpPhoYceAlD9ZSyXy7FgwYIaPTX19Zzoe1tubXP48GHEx8cbtXN0dASAWoPD7fSXE2+fubRixQoAwMMPP3zHc5iiZ8+eGDJkCD755BNs27atxvPl5eV4/fXXDY8jIiJw7tw5XLt2zbDv999/x8GDBxv0um5uboiNjcUXX3yBzZs3Q6lUYsSIEUZtnnzyScTHx2PHjh01js/Ly0NlZWWDXpPIXDjGh6gFSklJwSOPPAKZTIZRo0bhyy+/NHq+U6dOJg2IPn/+PDZt2gQhBLRarWHl5sLCQqxYsQJDhgyp99iBAwfiySefRFRUFOzt7fHNN98gMzPT0DuwceNG/Pvf/8Zjjz2GiIgIFBQU4OOPP4arq6shPAwfPhwDBgzA7NmzkZycjM6dO2Pnzp349ttv8eqrryIiIgIA0Lp1a8yePRsLFy5Enz598Pjjj0OlUuHo0aPw9/fHokWLaq1z2LBh2Lp1Kx577DE8/PDDuHTpElavXo2oqCijy20ODg6IiorCli1b0KZNG3h4eKBDhw5GtwrR69y5M8aOHYuPPvoIeXl56NevH44cOYKNGzdixIgRGDBgwB1/9qb69NNPMXjwYDz++OMYPnw4Bg4cCCcnJ/z555/YvHkz0tPTDWv5jB8/HitWrEBsbCwmTJiArKwsrF69Gu3btzcMlDbVU089hWeeeQb//ve/ERsbW2OM0fTp0/Hdd99h2LBhhqn5RUVFOHnyJL766iskJyff9bIKRM1CyillRFRNP529rmnf/fr1M5rOvnv37nqnvtc2Jft2t7aXy+XCzc1NREdHi1deeUWcPn26Rvvbp7Nfv35dTJkyRbRt21Y4OTkJjUYjevToIb744gvDMcePHxdPP/20CA4OFiqVSnh7e4thw4aJY8eOGZ27oKBATJ06Vfj7+wuFQiHuuecesWzZMqNp6nrr1q0T0dHRQqVSCXd3d9GvXz8RFxdn9LO6dUq3TqcT//znP0VISIhQqVQiOjpabN++vcYUbyGEOHTokIiJiRFKpdLo53j7dHYhhKioqBDz588XYWFhQqFQiKCgIDFr1ixRWlpq1C4kJEQ8/PDDNd5HXVPPa1NcXCyWL18uunXrJpydnYVSqRT33HOP+Pvf/26Y7q+3adMmER4eLpRKpejSpYvYsWNHndPZ61uKQKvVCgcHBwFAbNq0qdY2BQUFYtasWaJ169ZCqVQKT09P0atXL7F8+XJRXl5u0nsjMjfeq4uIiIhsBsf4EBERkc1g8CEiIiKbweBDRERENoPBh4iIiGwGgw8RERHZDAYfIiIishlcwPA2Op0OV69ehYuLS4OWsCciIiLpCCFQUFAAf3//eleyZ/C5zdWrVxEUFCR1GURERNQIly9fRmBgYJ3PM/jcRn9DxsuXL8PV1VXiaoiIiMgUWq0WQUFBd7yxMoPPbfSXt1xdXRl8iIiIWpg7DVPh4GYiIiKyGQw+REREZDMYfIiIiMhmMPgQERGRzWDwISIiIpvB4ENEREQ2g8GHiIiIbAaDDxEREdkMBh8iIiKyGQw+REREZDMYfIiIiMhmMPgQERGRzWDwMZP8kgqcTddCCCF1KURERDaLwccMSiuq0Hn+Tjz03n7kl1RIXQ4REZHNYvAxA7XCDp7OSgDAldwSiashIiKyXQw+ZuLv5gAAuJrH4ENERCQVBh8zCbgRfNIYfIiIiCTD4GMm7PEhIiKSHoOPmbDHh4iISHoMPmbibwg+pRJXQkREZLsYfMwk0P1G8OGsLiIiIskw+JiJvsfnemEZSiuqJK6GiIjINlll8Fm1ahVCQ0OhVqvRo0cPHDlyROqS4O6ogIPCDgCQns/LXURERFKwuuCzZcsWvPbaa5g7dy6OHz+Ozp07IzY2FllZWZLWJZPJ4O+mBsCZXURERFKxuuCzYsUKTJw4Ec899xyioqKwevVqODo6Yt26dVKXhgB3RwAc50NERCQVqwo+5eXl+O233zBo0CDDPrlcjkGDBiE+Pr7WY8rKyqDVao225hJwo8eHU9qJiIikYVXB5/r166iqqoKPj4/Rfh8fH2RkZNR6zKJFi6DRaAxbUFBQs9XHtXyIiIikZVXBpzFmzZqF/Px8w3b58uVmey2u3kxERCQte6kLaEqenp6ws7NDZmam0f7MzEz4+vrWeoxKpYJKpTJHeezxISIikphV9fgolUrExMTg559/NuzT6XT4+eef0bNnTwkrq6bv8UnPK4VOJySuhoiIyPZYVY8PALz22msYO3Ys7r33XnTv3h3vvvsuioqK8Nxzz0ldGnw1ashlQHmVDtcLy+Dtqpa6JCIiIptidcHnqaeewrVr1zBnzhxkZGSgS5cu+Omnn2oMeJaCwk4OH1c10vNLcSWvhMGHiIjIzKzqUpfeSy+9hJSUFJSVleHw4cPo0aOH1CUZBHCAMxERkWSsMvhYsgDerJSIiEgyDD5mxintRERE0mHwMTNOaSciIpIOg4+Z3Qw+vEM7ERGRuTH4mNnNMT7FEldCRERkexh8zEw/xkdbWomC0gqJqyEiIrItDD5m5qyyh8ZBAQC4ystdREREZsXgIwHO7CIiIpIGg48E9AOcrzD4EBERmRWDjwQC3KpvVcEeHyIiIvNi8JEAV28mIiKSBoOPBDjGh4iISBoMPhLg6s1ERETSYPCRgD74ZGpLUVGlk7gaIiIi28HgIwFPZxWUdnLoBJCRz7V8iIiIzIXBRwJyuQx+nNlFRERkdgw+EuE4HyIiIvNj8JEIZ3YRERGZH4OPRNjjQ0REZH4MPhK5GXw4uJmIiMhcGHwkcnP15mKJKyEiIrIdDD4SuTnGpxRCCImrISIisg0MPhLx01RPZy+pqEJucYXE1RAREdkGBh+JqBV28HRWAeDMLiIiInNh8JGQfpzPFd6lnYiIyCwYfCQUcGP1Zk5pJyIiMg8GHwkFcBFDIiIis2LwkZBhLR9e6iIiIjILBh8JGaa05zP4EBERmQODj4RuLmLI4ENERGQODD4S0l/qyi4qR2lFlcTVEBERWT8GHwlpHBRwUtoB4MwuIiIic2DwkZBMJrvl1hUMPkRERM2NwUdiHOdDRERkPgw+EmOPDxERkfkw+EhMP8D5CoMPERFRs2PwkRhXbyYiIjIfBh+JGcb4MPgQERE1OwYfienH+GTkl6JKJySuhoiIyLox+EjMx0UFO7kMFVUC1wrKpC6HiIjIqllV8AkNDYVMJjPaFi9eLHVZ9bK3k8PXVQ2Al7uIiIiam73UBTS1BQsWYOLEiYbHLi4uElZjmgA3B6TllSAtrwQxIe5Sl0NERGS1rC74uLi4wNfXV+oyGsTfrbrHhzO7iIiImpdVXeoCgMWLF6NVq1aIjo7GsmXLUFlZWW/7srIyaLVao83cuHozERGReVhVj8/LL7+Mrl27wsPDA4cOHcKsWbOQnp6OFStW1HnMokWLMH/+fDNWWRNXbyYiIjIPmRDCoudQz5w5E0uWLKm3zdmzZ9G2bdsa+9etW4fJkyejsLAQKpWq1mPLyspQVnZzNpVWq0VQUBDy8/Ph6up6d8WbaE9iFsatP4q2vi746dW+ZnlNIiIia6LVaqHRaO74/W3xPT7Tpk3DuHHj6m0THh5e6/4ePXqgsrISycnJiIyMrLWNSqWqMxSZi371Zs7qIiIial4WH3y8vLzg5eXVqGMTEhIgl8vh7e3dxFU1Lf2lroLSSmhLK+CqVkhcERERkXWy+OBjqvj4eBw+fBgDBgyAi4sL4uPjMXXqVDzzzDNwd7fsKeJOKnu4OSqQV1yBq3klcPVl8CEiImoOVhN8VCoVNm/ejHnz5qGsrAxhYWGYOnUqXnvtNalLM0mAmwPyiiuQlluCtr7mGVtERERka6wm+HTt2hW//vqr1GU0mr+bA05f1XJmFxERUTOyunV8Wir9AOcrDD5ERETNhsHHQgQY1vIplbgSIiIi68XgYyFurt5cLHElRERE1ovBx0JwLR8iIqLmx+BjIfRr+WQVlKG8UidxNURERNaJwcdCeDorobSXQwggI5/jfIiIiJoDg4+FkMlkvNxFRETUzBh8LAiDDxERUfNi8LEg/m5qAOAihkRERM2EwceCBLg5AgDSchl8iIiImgODjwUx9PjkM/gQERE1BwYfC3JzEUMGHyIioubA4GNBbh3cLISQuBoiIiLrw+BjQXw1ashkQFmlDtlF5VKXQ0REZHUYfCyIyt4OXs4qAJzZRURE1BwYfCwMx/kQERE1HwYfC+PPRQyJiIiaDYOPhQlk8CEiImo2DD4WRt/jwzE+RERETY/Bx8Lwfl1ERETNh8HHwtzs8SmVuBIiIiLrw+BjYfSzunKKylFcXilxNURERNaFwcfCuKrt4ayyB8BeHyIioqbG4GNhZDIZx/kQERE1EwYfC2S4SzuDDxERUZNi8LFAXL2ZiIioeTD4WCCu5UNERNQ8GHwskH6MzxUGHyIioibF4GOBAtjjQ0RE1CwYfCyQfoxPRn4pqnRC4mqIiIisB4OPBfJ2UcNeLkOlTiBTy7V8iIiImgqDjwWyk8vgq+GUdiIioqbG4GOhuIghERFR02PwsVAMPkRERE2PwcdCcRFDIiKipsfgY6G4iCEREVHTY/CxULzURURE1PQYfCyUvscnLbcEQnAtHyIioqbA4GOh9D0+ReVV0JZUSlwNERGRdWDwsVAOSjt4OCkB8HIXERFRU2kxweftt99Gr1694OjoCDc3t1rbpKam4uGHH4ajoyO8vb0xffp0VFa23N4SjvMhIiJqWi0m+JSXl+OJJ57Aiy++WOvzVVVVePjhh1FeXo5Dhw5h48aN2LBhA+bMmWPmSpuOvxtXbyYiImpKLSb4zJ8/H1OnTkXHjh1rfX7nzp04c+YMNm3ahC5duuChhx7CwoULsWrVKpSXl5u52qYR4OYIgD0+RERETaXFBJ87iY+PR8eOHeHj42PYFxsbC61Wi9OnT9d5XFlZGbRardFmKfQ9Pgw+RERETcNqgk9GRoZR6AFgeJyRkVHncYsWLYJGozFsQUFBzVpnQwRy9WYiIqImJWnwmTlzJmQyWb3buXPnmrWGWbNmIT8/37Bdvny5WV+vIbh6MxERUdOyl/LFp02bhnHjxtXbJjw83KRz+fr64siRI0b7MjMzDc/VRaVSQaVSmfQa5qaf1ZVVUIayyiqo7O0kroiIiKhlkzT4eHl5wcvLq0nO1bNnT7z99tvIysqCt7c3ACAuLg6urq6IiopqktcwNw8nJdQKOUordMjIL0VIKyepSyIiImrRWswYn9TUVCQkJCA1NRVVVVVISEhAQkICCgsLAQCDBw9GVFQUxowZg99//x07duzAW2+9hSlTplhsj86dyGQyo1tXEBER0d2RtMenIebMmYONGzcaHkdHRwMAdu/ejf79+8POzg7bt2/Hiy++iJ49e8LJyQljx47FggULpCq5SQS4OeDitSLO7CIiImoCLSb4bNiwARs2bKi3TUhICH744QfzFGQmXL2ZiIio6bSYS122ijO7iIiImg6Dj4Vjjw8REVHTYfCxcDd7fEolroSIiKjlY/CxcIbVm/NKoNMJiashIiJq2Rh8LJyPqxoyGVBeqUN2Ucu82SoREZGlYPCxcEp7OXxceLNSIiKipsDg0wIY7tLORQyJiIjuCoNPCxDg7giAU9qJiIjuVoMWMNTpdNi7dy/279+PlJQUFBcXw8vLC9HR0Rg0aBCCgoKaq06bxintRERETcOkHp+SkhL83//9H4KCgjB06FD8+OOPyMvLg52dHZKSkjB37lyEhYVh6NCh+PXXX5u7ZpsT4MYxPkRERE3BpB6fNm3aoGfPnvj444/x4IMPQqFQ1GiTkpKC//73vxg9ejRmz56NiRMnNnmxtirAnTcqJSIiagomBZ+dO3eiXbt29bYJCQnBrFmz8PrrryM1NbVJiqNqhkUM8xl8iIiI7oZJl7ratWuHU6dOmXRChUKBiIiIuyqKjOnH+OQVV6CorFLiaoiIiFouk2d1derUCT169MDHH3+MgoKC5qyJbuOiVsBFXd05x5ldREREjWdy8Nm7dy/at2+PadOmwc/PD2PHjsX+/fubsza6hb7X5wqDDxERUaOZHHz69OmDdevWIT09He+//z6Sk5PRr18/tGnTBkuWLEFGRkZz1mnzAgw3K2XwISIiaqwGL2Do5OSE5557Dnv37sX58+fxxBNPYNWqVQgODsYjjzzSHDUSOLOLiIioKdzVys2tW7fGm2++ibfeegsuLi74/vvvm6ouuo0/e3yIiIjuWoNWbr7Vvn37sG7dOnz99deQy+V48sknMWHChKasjW7B1ZuJiIjuXoOCz9WrV7FhwwZs2LABSUlJ6NWrF/71r3/hySefhJOTU3PVSLi1x6dU4kqIiIhaLpODz0MPPYRdu3bB09MTzz77LMaPH4/IyMjmrI1uEXhjjE+GthSVVTrY2/H+skRERA1lcvBRKBT46quvMGzYMNjZ2TVnTVQLL2cVFHYyVFQJZBaUGS59ERERkelMDj7fffed0eOkpCRcuHABffv2hYODA4QQkMlkTV4gVZPLZfDTOCA1pxhpuSUMPkRERI3Q4Osl2dnZGDhwINq0aYOhQ4ciPT0dADBhwgRMmzatyQukm/xv3KWdM7uIiIgap8HBZ+rUqVAoFEhNTYWjo6Nh/1NPPYWffvqpSYsjYwFu1T9vzuwiIiJqnAZPZ9+5cyd27NiBwMBAo/333HMPUlJSmqwwqingRo8Pgw8REVHjNLjHp6ioyKinRy8nJwcqlapJiqLacfVmIiKiu9Pg4NOnTx98+umnhscymQw6nQ5Lly7FgAEDmrQ4MsbVm4mIiO5Ogy91LV26FAMHDsSxY8dQXl6OGTNm4PTp08jJycHBgwebo0a64dbVmzmLjoiIqOEa3OPToUMHnD9/Hr1798ajjz6KoqIiPP744zhx4gQiIiKao0a6Qd/jU1xehfySComrISIiankada8ujUaD2bNnN3UtdAdqhR08nZW4XliOK7klcHNUSl0SERFRi2JSj09qamqDTpqWltaoYujOOM6HiIio8UwKPt26dcPkyZNx9OjROtvk5+fj448/RocOHfD11183WYFkjHdpJyIiajyTLnWdOXMGb7/9Nh588EGo1WrExMTA398farUaubm5OHPmDE6fPo2uXbti6dKlGDp0aHPXbbP0PT6c0k5ERNRwJvX4tGrVCitWrEB6ejo++OAD3HPPPbh+/Tr+/PNPAMBf//pX/Pbbb4iPj2foaWb6Hp+r+Qw+REREDdWgwc0ODg4YNWoURo0a1Vz10B1wEUMiIqLGa/B0dpLWzTE+pRJXQkRE1PIw+LQw+uBzvbAMpRVVEldDRETUsjD4tDBujgo4KOwAAOn57PUhIiJqCAafFkYmk3GcDxERUSM16u7sUnj77bfRq1cvODo6ws3NrdY2MpmsxrZ582bzFmoGXMSQiIiocRocfHx8fDB+/HgcOHCgOeqpU3l5OZ544gm8+OKL9bZbv3490tPTDduIESPMU6AZ6cf5XGHwISIiapAGB59NmzYhJycHDzzwANq0aYPFixfj6tWrzVGbkfnz52Pq1Kno2LFjve3c3Nzg6+tr2NRqdbPXZm4BbtXviT0+REREDdPg4DNixAhs27YNaWlpeOGFF/Df//4XISEhGDZsGLZu3YrKysrmqNNkU6ZMgaenJ7p3745169ZBCFFv+7KyMmi1WqPN0nGMDxERUeM0enCzl5cXXnvtNfzxxx9YsWIFdu3ahVGjRsHf3x9z5sxBcXFxU9ZpkgULFuCLL75AXFwcRo4cib/97W94//336z1m0aJF0Gg0hi0oKMhM1Taev4arNxMRETVGo4NPZmYmli5diqioKMycOROjRo3Czz//jHfeeQdbt241aWzNzJkzax2QfOt27tw5k2v6xz/+gfvvvx/R0dF44403MGPGDCxbtqzeY2bNmoX8/HzDdvnyZZNfTyr6Hp/0vFLodPX3aBEREdFNDbplBQBs3boV69evx44dOxAVFYW//e1veOaZZ4xmWvXq1Qvt2rW747mmTZuGcePG1dsmPDy8oSUa9OjRAwsXLkRZWRlUKlWtbVQqVZ3PWSofVzXkMqC8SofrhWXwdrW+cUxERETNocHB57nnnsPo0aNx8OBBdOvWrdY2/v7+mD179h3P5eXlBS8vr4aWYLKEhAS4u7u3uGBzJwo7OXxd1biaX4oreSUMPkRERCZqcPBJT0+Ho6NjvW0cHBwwd+7cRhdVm9TUVOTk5CA1NRVVVVVISEgAALRu3RrOzs743//+h8zMTNx3331Qq9WIi4vDP//5T7z++utNWoel8HdzwNX8UlzNK0HXYHepyyEiImoRGhx8Kisra535JJPJoFKpoFQqm6Sw282ZMwcbN240PI6OjgYA7N69G/3794dCocCqVaswdepUCCHQunVrrFixAhMnTmyWeqQW4O6AYym5nNlFRETUAA0OPm5ubpDJZHU+HxgYiHHjxmHu3LmQy5vujhgbNmzAhg0b6nx+yJAhGDJkSJO9nqXj6s1EREQN1+Dgs2HDBsyePRvjxo1D9+7dAQBHjhzBxo0b8dZbb+HatWtYvnw5VCoV3nzzzSYvmKrpV29OY/AhIiIyWYODz8aNG/HOO+/gySefNOwbPnw4OnbsiDVr1uDnn39GcHAw3n77bQafZnQz+PAO7URERKZq8LWoQ4cOGcbX3Co6Ohrx8fEAgN69eyM1NfXuq6M63Vy92fwLRRIREbVUDQ4+QUFBWLt2bY39a9euNax6nJ2dDXd3zjRqTvoxPtrSShSUVkhcDRERUcvQ4Etdy5cvxxNPPIEff/zRsI7PsWPHcO7cOXz11VcAgKNHj+Kpp55q2krJiLPKHhoHBfJLKnA1rxSRvgqpSyIiIrJ4DQ4+jzzyCBITE7FmzRokJiYCAB566CFs27YNoaGhAIAXX3yxSYuk2vm7OdwIPiWI9HWRuhwiIiKL16DgU1FRgSFDhmD16tVYtGhRc9VEJgpwc8DZdC2ucGYXERGRSRo0xkehUOCPP/5orlqogQLcqm9VwbV8iIiITNPgwc3PPPNMrYObyfxuzuxi8CEiIjJFo25ZsW7dOuzatQsxMTFwcnIyen7FihVNVhzVz5+LGBIRETVIg4PPqVOn0LVrVwDA+fPnjZ6r71YW1PQCeNsKIiKiBmlw8Nm9e3dz1EGNoL/UlaktRUWVDgq7prs3GhERkTVq9DdlUlISduzYgZKS6t4GIUSTFUWm8XRSQWknh04AGfm8dQUREdGdNDj4ZGdnY+DAgWjTpg2GDh2K9PR0AMCECRMwbdq0Ji+Q6iaXy+B/Y2YXx/kQERHdWYODz9SpU6FQKJCamgpHR0fD/qeeego//fRTkxZHd+bPcT5EREQma/AYn507d2LHjh0IDAw02n/PPfcgJSWlyQoj0xju0s4p7URERHfU4B6foqIio54evZycHKhUqiYpikxn6PHJZ/AhIiK6kwYHnz59+uDTTz81PJbJZNDpdFi6dCkGDBjQpMXRnelndl1hjw8REdEdNfhS19KlSzFw4EAcO3YM5eXlmDFjBk6fPo2cnBwcPHiwOWqkenAtHyIiItM1uMenQ4cOOH/+PHr37o1HH30URUVFePzxx3HixAlEREQ0R41Uj4BbVm/mkgJERET1a3CPDwBoNBrMnj27qWuhRvDVVE9nL63QIbe4Ah5OSokrIiIislyNCj55eXk4cuQIsrKyoNPpjJ579tlnm6QwMo1aYQcvFxWuFZQhLbeEwYeIiKgeDQ4+//vf//DXv/4VhYWFcHV1Nbo/l0wmY/CRgL+bQ3XwyStBx0CN1OUQERFZrAaP8Zk2bRrGjx+PwsJC5OXlITc317Dl5OQ0R410B4G8SzsREZFJGhx80tLS8PLLL9e6lg9JQ3/bCs7sIiIiql+Dg09sbCyOHTvWHLVQI3H1ZiIiItM0eIzPww8/jOnTp+PMmTPo2LEjFAqF0fOPPPJIkxVHpuHqzURERKZpcPCZOHEiAGDBggU1npPJZKiqqrr7qqhB9Ks3s8eHiIiofg0OPrdPXyfp6S91ZReVo7SiCmqFncQVERERWaYGj/Ehy6NxUMBJWR12OLOLiIiobiYHn6FDhyI/P9/wePHixcjLyzM8zs7ORlRUVJMWR6aRyWQ3x/kw+BAREdXJ5OCzY8cOlJWVGR7/85//NFq3p7KyEomJiU1bHZmM43yIiIjuzOTgc/sNMHlDTMvCHh8iIqI74xgfK6Ef4HyFwYeIiKhOJgcfmUxmdF8u/T6yDAHs8SEiIrojk6ezCyEwbtw4qFQqAEBpaSleeOEFODk5AYDR+B8yP8MYHwYfIiKiOpkcfMaOHWv0+JlnnqnRhndml45+jE96XimqdAJ2cvbGERER3c7k4LN+/frmrIPuko+LCnZyGSp1AtcKyuCrUUtdEhERkcXh4GYrYW8nh69rddhJyyuWuBoiIiLLxOBjRQx3ac8rlbgSIiIiy9Qigk9ycjImTJiAsLAwODg4ICIiAnPnzkV5eblRuz/++AN9+vSBWq1GUFAQli5dKlHF0uAihkRERPVr8E1KpXDu3DnodDqsWbMGrVu3xqlTpzBx4kQUFRVh+fLlAACtVovBgwdj0KBBWL16NU6ePInx48fDzc0NkyZNkvgdmIe/W/WlLk5pJyIiql2LCD5DhgzBkCFDDI/Dw8ORmJiIDz/80BB8PvvsM5SXl2PdunVQKpVo3749EhISsGLFCpsJPgFujgA4pZ2IiKguLeJSV23y8/Ph4eFheBwfH4++fftCqVQa9sXGxiIxMRG5ubl1nqesrAxardZoa6nY40NERFS/Fhl8kpKS8P7772Py5MmGfRkZGfDx8TFqp3+ckZFR57kWLVoEjUZj2IKCgpqnaDMI5BgfIiKiekkafGbOnGm4FUZd27lz54yOSUtLw5AhQ/DEE09g4sSJd13DrFmzkJ+fb9guX7581+eUin4Rw4KySmhLKySuhoiIyPJIOsZn2rRpGDduXL1twsPDDX+/evUqBgwYgF69euGjjz4yaufr64vMzEyjffrHvr6+dZ5fpVIZbsPR0jkq7eHuqEBucQXSckvg6qeQuiQiIiKLImnw8fLygpeXl0lt09LSMGDAAMTExGD9+vWQy407q3r27InZs2ejoqICCkX1F35cXBwiIyPh7u7e5LVbKn83B+QWV+BqXgna+blKXQ4REZFFaRFjfNLS0tC/f38EBwdj+fLluHbtGjIyMozG7vzlL3+BUqnEhAkTcPr0aWzZsgXvvfceXnvtNQkrN7+bixhynA8REdHtWsR09ri4OCQlJSEpKQmBgYFGzwkhAAAajQY7d+7ElClTEBMTA09PT8yZM8dmprLr+TP4EBER1alFBJ9x48bdcSwQAHTq1An79+9v/oIsGGd2ERER1a1FXOoi0+l7fLiWDxERUU0MPlaGY3yIiIjqxuBjZfQ9PlkFZSiv1ElcDRERkWVh8LEyns5KKO3lEALIyC+VuhwiIiKLwuBjZWQyGS93ERER1YHBxwox+BAREdWOwccK8S7tREREtWPwsUIBbo4AuJYPERHR7Rh8rJChxyefwYeIiOhWDD5WKICrNxMREdWKwccK3Tq4WX8vMyIiImLwsUq+GjVkMqCsUofsonKpyyEiIrIYDD5WSGVvBy9nFQBe7iIiIroVg4+V0o/z4ZR2IiKimxh8rBQXMSQiIqqJwcdKMfgQERHVxOBjpTilnYiIqCYGHyvlr7kxxoeLGBIRERkw+Fgp9vgQERHVxOBjpfxvjPHJLa5AcXmlxNUQERFZBgYfK6VxUMBFZQ+AU9qJiIj0GHysmL9hZlepxJUQERFZBgYfK8ZxPkRERMYYfKyYv5saAC91ERER6TH4WLEAN0cAXMSQiIhIj8HHiul7fBh8iIiIqjH4WLFAjvEhIiIywuBjxfSzujK0pajSCYmrISIikh6DjxXzdlHDXi5DlU4gnbeuICIiYvCxZnZyGdr5uQIAvk24KnE1RERE0mPwsXLP3R8KANhwKBlllVXSFkNERCQxBh8rN6yTP3xd1bhWUIZvT7DXh4iIbBuDj5VT2ssxvncoAOCj/Reh4yBnIiKyYQw+NuDp7sFwUdkjKasQuxOzpC6HiIhIMgw+NsBFrcBfegQDANbsuyhxNURERNJh8LERz90fBoWdDEcu5SDhcp7U5RAREUmCwcdG+GrUeKRzAADgo30XJK6GiIhIGgw+NmRS33AAwE+nMpCSXSRxNURERObH4GNDIn1d0D/SCzoBfLL/ktTlEBERmR2Dj43R9/p8+dtl5BSVS1wNERGRebWI4JOcnIwJEyYgLCwMDg4OiIiIwNy5c1FeXm7URiaT1dh+/fVXCSu3PD3DW6FDgCtKK3T4T3yK1OUQERGZVYsIPufOnYNOp8OaNWtw+vRprFy5EqtXr8abb75Zo+2uXbuQnp5u2GJiYiSo2HLJZDJM6hsBAPg0PhmlFbyNBRER2Q57qQswxZAhQzBkyBDD4/DwcCQmJuLDDz/E8uXLjdq2atUKvr6+5i6xRRnawRdL3R1wJbcEX/12Bc/cFyJ1SURERGbRInp8apOfnw8PD48a+x955BF4e3ujd+/e+O677+54nrKyMmi1WqPN2tnbyTGhdxgA4JP9F1HF21gQEZGNaJHBJykpCe+//z4mT55s2Ofs7Ix33nkHX375Jb7//nv07t0bI0aMuGP4WbRoETQajWELCgpq7vItwpP3BkHjoEBydjHizmRIXQ4REZFZyIQQkv3v/syZM7FkyZJ625w9exZt27Y1PE5LS0O/fv3Qv39/fPLJJ/Ue++yzz+LSpUvYv39/nW3KyspQVlZmeKzVahEUFIT8/Hy4urqa+E5apuU7EvHB7iREB7th64u9IJPJpC6JiIioUbRaLTQazR2/vyUd4zNt2jSMGzeu3jbh4eGGv1+9ehUDBgxAr1698NFHH93x/D169EBcXFy9bVQqFVQqlUn1WpuxvULx0f6LOJGah2MpuegWWvPSIRERkTWRNPh4eXnBy8vLpLZpaWkYMGAAYmJisH79esjld75Kl5CQAD8/v7st02p5uagwsmsAPj9yGWv2XmTwISIiq9ciZnWlpaWhf//+CAkJwfLly3Ht2jXDc/oZXBs3boRSqUR0dDQAYOvWrVi3bt0dL4fZuuf7hOPzI5ex62wmkrIK0drbWeqSiIiImk2LCD5xcXFISkpCUlISAgMDjZ67dYjSwoULkZKSAnt7e7Rt2xZbtmzBqFGjzF1uixLh5YxB7Xyw62wmPtl/EYtHdpK6JCIiomYj6eBmS2Tq4ChrcjQ5B0+sjofSTo4DMwfA20UtdUlEREQN0iIGN5NluDfEHdHBbjiRmodPD6Xg9dhIqUsiIqIWqLSiCtcKypCpLUXWjT8ztWXIKihF1o0/tSWViJ/1gGQziRl8CDKZDJP7huOFTcfxn19T8GL/CDip+J8GERFVqy3QGP68EWgytWXIL6kw6XwFZZVwVSuauera8duNAAAPRvkitJUjkrOL8cWxy3ju/jCpSyIiomamDzT64HJroLl2S4+NqYEGAFT2cni7quDjooa3qwreLmr4uKrh7aKq/tNVBUeFXTO+q/ox+BAAwE4uw/N9wvHWtlNYe+ASxtwXAnu7FrmwNxERAaio0iFTW4qreaW4mleCtLwSXL2xpeeXIkNbirzipg00Pi5quDrYW/SCuAw+ZDAqJhAr487jSm4JfjiVgUc6+0tdEhER1UIIgbziClzNLzEEG+NwU4qsglKYcitGpb0cPrcFGn2IaUmBxlQMPmSgVtjh2Z6hWLnrPD7adwHDO/lZxX/kREQtTWlFFTLybwSa/NqDTUlF1R3Po7STw89NDX+NA/zdHBDgpoa/mwP83Bzgp1FbVaAxFYMPGRnTMwQf7k3CqTQt4i9ko1drT6lLIiKyOoVllUjNLkZqThGu5FZferoZbkpxvbDszicB4Omsgv8twcbfTY0AN/3fHdDKSQm53HZCjSkYfMiIh5MST8QE4T+/pmDNvosMPkREjaDTCWQVlCE1pxgp2UW4nFOMlJxipGQX43JOMbKLyu94DgeFXXWocXOoNdj4atRQSzhIuKVi8KEanu8Ths8Op2Dv+WtIzChApK+L1CUREVmc0ooqXMktvhFuqv9Mza4OOJdzilFWqav3eA8nJYI8HBHkfiPUaNSGnpoANwe4OSps6hKUuTD4UA0hrZwwpIMvfjiZgY/2XcQ7T3aWuiQiIrMTQiC3uAIp2UWGUJOaczPYZGhLUd+9D+zkMgS4OSDYwxHBrRwR7OGIEA9HBN14LNU6NraOwYdqNalvBH44mYHvfk/D9NhI+Gp4Gwsisj5CCGRoS3EhqwgpOUU3w82NS1IFZZX1Hu+ktENwKyeEeDgipFV1qAm5EXL83Ryg4LIgFofBh2rVJcgN3cM8cORSDtYfvIRZQ9tJXRIRUaOVVlQhObsIF68V4UJWIS5cK8SFa0W4eK0QReX1z47ydVUbem1Cbum9CfZwhIeTkpejWhgGH6rT5L7hOHIpB/89nIqXHmgNF3bLEpEFE0Igu6i8OtxcKzQKOJdzi+u8LGUnlyHEwxGhnk6GQBPSqnoLdHfkAGIrw+BDdRoQ6Y3W3s5IyirE50dSMalvhNQlERGhokqHyznFuFBLwKnv1gouantEeDlXb95Ohr8HezhCac9LUraCwYfqJJfLMKlPOGZ8/QfWHUjGuF5h/OVARGaTX1KBizcCza0BJyW7GJV1LEkskwGB7g4I96wZcDydeVmKGHzoDh6N9seynYnI0Jbif79fxciYQKlLIiIrk19SgcSMAiRmaHEuowBJWYW4eL0I1wrqXsTPQWGHCG+nGgEnzNOJl6aoXgw+VC+VvR3G9QrFsh2J+Hj/RTzeNYD/x0REjVJRpcPFa0U4dyPgJGYU4Fy6FlfzS+s8xtdVjXAvfa+NEyK8q4OOr6uaKxJTozD40B090yMEq3Yn4VxGAfaev4b+kd5Sl0REFkwIgUxtmVHAOZuuxYVrhaioqv0SVYCbAyJ9XRDp64I2Ps6G3htOqqCmxuBDd6RxVGB0t2CsO3gJH++/yOBDRAZFZZU4n1lgFHASMwuQV1z7IGNnlT0ifV3QVr/5uaKNjws0Dgw4ZB4MPmSS8b1DsTE+GQeTsnEqLR8dAjRSl0REZlSlE0jJLsK5jOqQc+5GwEnJLq61vZ1chjBPp5sBx9cVkb4uCHR34OVykhSDD5kk0N0Rwzr54duEq/ho30X86+loqUsiomaSX1yB01fzcfaWgHM+swClFbXfe8rLRWUIOJG+rmjr64LW3s4cZEwWicGHTDapbzi+TbiK70+mY8aQSAS6O0pdEhHdJW1pBU6l5ePklXycTKve6urFUSvkaONzM+C0uzEmp5WzysxVEzUegw+ZrL2/Br1be+JA0nWsPXAJc4e3l7okImoAbWkFTqdpcTItDyfTtDh5JQ/JdYScIA8HRPm5GgWckFZOsONMKmrhGHyoQSb1DceBpOvYcvQyXh3YBhpHDkgkskQFpRU4fVVr1JNz6XpRrW0D3R3QMUCDjoEadAzQoIO/Bu5OSjNXTGQeDD7UIH3u8URbXxecyyjApsMpmDKgtdQlEdm8wrJKnE67GXBOpuXj4rXaQ06A220hJ0ADD4YcsiEMPtQgMpkMk/qG47Uvfsf6g8mY0DuMAxiJzKiorBKnr2rxx5W86rE5afm4eL2o1htwBrg5oEOA642g44aODDlEDD7UcMM7+2PZjkSk55di24k0jO4eLHVJRFaptKIKp9Ly8cctl6suXCusNeT4adTVAeeW3hwOOiaqicGHGkxhJ8f4+8Pw9g9n8dH+i3jy3iAuHU90l4QQSMsrwfHUPBxPycWJ1Fycvqqt9Wacvq5qQ7jRX67ycmHIITIFgw81yujuQfjXz3/i4rUi/HIuC4OifKQuiahFKa2owumr+Tiekofjqbk4npqLTG3Nm3J6OqvQJUiDjgFu6Bjoig4BGni7qCWomMg6MPhQo7ioFfjLfcFYs/ciPtp3kcGH6A7S80uMQs7pNC3Kq4wXBLSTyxDl54quwW7oGuKOrsHuXOmYqIkx+FCjjb8/DOsOXMKR5BycSM1FdLC71CURWYTySl11b05qddA5kZJb6x3IWzkpDQGna7AbOgW6wUHJyQJEzYnBhxrNx1WNR7sE4KvfruCjfRfx4TMxUpdEJIksbemNnpzq8Tl/pOWjvNK4N0cuA9r5uVaHnBA3dA12R7CHI3tziMyMwYfuyqS+4fjqtyv46XQGkq8XIdTTSeqSiJpVRZUOZ9O1OJ5SHXR+S8lFWl5JjXbujoobIccd0cFu6BzoBicVf+USSY3/CumutPFxwYBIL+xOvIZPDlzE/43oKHVJRE1KW1qBY8k5OHKpemzOH1fyatysUy6r/regv2wVE+KO0FbszSGyRAw+dNcm9g3H7sRr+PLYFUwd1IZrh1CLll1YhqPJOTh8KQdHLuXgTLq2xro5GgcFooPdDCGnc5AbnNmbQ9Qi8F8q3bWe4a3QMUCDk2n5+DQ+BVMfbCN1SUQmS88vwZFLN4NOUlZhjTZhnk7oFuqOe0M90DXYHeGeTly7iqiFYvChu6a/jcXfPz+BT+OT8UK/CM5MIYskhEBqTrEh5By5lIPUnJp3J4/0cUGPcA90D/NA91APeLty3Rwia8HgQ03ioQ6+CHR3wJXcEnx1/ArG3BcidUlEEEIgKasQhw09Otk1FgmUy4AOARp0D60OOt1CPXhnciIrxuBDTcLeTo7ne4dh3v/O4JP9F/GX7sGw46UAMrMqncDZdK0h5BxNzkVOUblRG4WdDJ0D3ap7c8I8EBPiDhe1QqKKicjcGHyoyTzZLQjv/vwnUrKLsfN0Bh7q6Cd1SWTlyit1OJmWf+OyVTaOJeeioKzSqI1aIUfXYHd0D/NAj7BWiA52g1rBS7FEtqrFBJ9HHnkECQkJyMrKgru7OwYNGoQlS5bA39/f0OaPP/7AlClTcPToUXh5eeHvf/87ZsyYIWHVtsVRaY8x94Xg/V+SsGbfRQzp4MvpvNSkSiuqcCI1D4cvZePIpRwcT82tMbXcWWWPe0Pd0SOsFbqHeaBjgAZKe7lEFRORpWkxwWfAgAF488034efnh7S0NLz++usYNWoUDh06BADQarUYPHgwBg0ahNWrV+PkyZMYP3483NzcMGnSJImrtx3P9gzFmn0XkXA5D0eTc9E9zEPqkqgFK6usQkJqHuIvZiP+QjZOXM6rsSKyu6PixmWrVugR5oF2fq68zEpEdZIJcfsKFS3Dd999hxEjRqCsrAwKhQIffvghZs+ejYyMDCiV1QMTZ86ciW3btuHcuXMmn1er1UKj0SA/Px+urq7NVb5Vm7X1JD4/kopB7bzxydhuUpdDLUh5pQ5/XMlD/IVsxF/Mxm8puSi7Leh4u6jQI7zVjUtXHmjt5cyp5URk8vd3i+nxuVVOTg4+++wz9OrVCwpF9aDE+Ph49O3b1xB6ACA2NhZLlixBbm4u3N1rv4FmWVkZyspuzvLQarXNW7wNeL5PGDYfTcWus1n45w9n8dIDreHKwaNUi8qq6jE6+h6dY8m5KKmoMmrj6azEfeGtcF94K/SMaIVwTydeQiWiRmtRweeNN97ABx98gOLiYtx3333Yvn274bmMjAyEhYUZtffx8TE8V1fwWbRoEebPn998RdugCC9nPHtfCDbGp+CjfRex9fgVTBsciSfvDeIlCBtXpRM4fTXf0KNz9FIOisqNg467o8IQcnqGt0Jrb2cGHSJqMpJe6po5cyaWLFlSb5uzZ8+ibdu2AIDr168jJycHKSkpmD9/PjQaDbZv3w6ZTIbBgwcjLCwMa9asMRx75swZtG/fHmfOnEG7du1qPX9tPT5BQUG81HWXhBDYk3gNC78/g4vXigAAUX6umDM8CveFt5K4OjIXnU7gTLoWv17Mxq8Xs3H4Ug4KSo1nXWkcFOgR5lEddCJaoY23Cy9dEVGDmXqpS9Lgc+3aNWRnZ9fbJjw83Ojyld6VK1cQFBSEQ4cOoWfPnnj22Weh1Wqxbds2Q5vdu3fjgQceQE5OTp09PrfjGJ+mVVGlw6fxKXhv13lob3zhDWnvizeHtkNwK0eJq6OmptMJnM8qqO7RuVAddPJLKozauKjt0SPMw9Cr087XlUGHiO5aixjj4+XlBS8vr0Ydq9NVD3jU99b07NkTs2fPRkVFhWHcT1xcHCIjI00OPdT0FHZyTOgdhseiA7AiLhH/PZyKn05n4JdzWZjQJwxTBrTmzR1bMP3KyPoxOocv5dRYMNBZZY9uoe7oGVE9Tqe9v4aXPIlIMi1iVtfhw4dx9OhR9O7dG+7u7rhw4QL+8Y9/IDMzE6dPn4ZKpUJ+fj4iIyMxePBgvPHGGzh16hTGjx+PlStXNmg6O3t8mte5DC0Wbj+Dg0nVPX1eLipMj43EqK6B/L/+FkCnE/gzqxBHknNw+GI2fr2Yg+uFxreAcFTa4d5QD9wX7mG4ga29HdfRIaLm1SIudZnq5MmTeOWVV/D777+jqKgIfn5+GDJkCN566y0EBAQY2t26gKGnpyf+/ve/44033mjQazH4ND8hBHadzcLb359Bcnb1DSI7BmgwZ3gUuoVy3R9LUlmlw+mr2uqVkZNzcDQ5B3nFxpeuVPZy3Bvqjp43Ll11CnSDgkGHiMzMqoKPOTH4mE9ZZRU2HkrG+z8nGW4zMKyTH2YNbYcANweJq7NNpRVV1YtP3gg6v6Xkovi2WVeOSjt0DXZHt9DqAcmdgzRQ2fMWEEQkLQafRmLwMb/rhWV4Z2ciNh+9DCGqexAm9w3HC/0j4Kjk+J/mVFBagd9ScnHkUnVvzu+X81FeZbxgoMZBgW6h7obVkdv7u7JHh4gsDoNPIzH4SOf01Xws+N8ZHL6UAwDwdVXjjYci8WjnAI7/aSLZhWU4mpx749JVNs5c1UJ3228AbxeV4c7l3cM8OL2ciFoEBp9GYvCRlhACO05n4O0fzuJyTgkAoEuQG+YMj0LXYM7Oa6ireSWG8TlHLuUgKauwRptgD8fqkBNaHXRCWjlywUAianEYfBqJwccylFZUYd3BS1j1S5JhZd8RXfzxxkNt4afh+J/aCCFw6XpRddC5EXau5JbUaNfGx9lw2ap7qAd8NWoJqiUialoMPo3E4GNZsrSlWLYjEV8dvwIhAAeFHV7oF4FJfcPhoLTtAbXa0gqczyjAybR8HL3Ro3O90HgNHTu5DB38XdHtRm9Ot1APuDvVXBCUiKilY/BpJAYfy3TySj7m/+80jqXkAgD8NWrMHNoOwzv5Wf1lmbLKKlzIKkJiphbnMgpwPqMAiRkFuJpfWqOt0l6OLkFuhstWXUPcuUAkEdkEBp9GYvCxXEIIfH8yHYt+OIe0vOpLOPeGuGPO8Ch0CnSTtrgmoNMJXM4tNoSbc5nVAefS9SJU3T4C+QY/jRqRvi7oFlrdm9MpUAO1wrZ7wojINjH4NBKDj+UrrajCx/su4t97LqCkonr8z6iYQEyPjYSPa8sYr3KtoAyJGQVIzCxAYoYWiRkFOJ9ZaHg/t3NR26OtrwsifV0Q6euKSB8XRPq4QOOoMHPlRESWicGnkRh8Wo6M/FIs/ekctp5IA1C9sN6UAa0xoXeYxfR6FJVVIjHzRg/OjUtU5zMLkH3b/az0lPZytPZyNoScNr4uaOvrAl9XtdVf0iMiuhsMPo3E4NPynEjNxYLtZ3AiNQ9A9YJ77o4KKO3lUNjJDX+q9I/t5FDYV/+ptJdVP76lndLw3K37ZTePr6NdWWUVzmcWGnpwEjMLDFPybyeTASEejsY9OL4uCG3lyPtaERE1AoNPIzH4tExCCHz3+1Us/vEc0msZ9CslLxdVdQ+Oz80enHu8XWx+VhoRUVMy9fub0z3IKshkMjzaJQCx7X3xZ2YhyiqrUF6pQ3mVDuWVOlRUCZRXVaGiUqCsSoeKG8/p/7zZ7pb2lbXt16FMf+xtbeUyoLW3840enBt/+rrAg9PHiYgsBoMPWRW1wg4dAzVSl0FERBaKgwmIiIjIZjD4EBERkc1g8CEiIiKbweBDRERENoPBh4iIiGwGgw8RERHZDAYfIiIishkMPkRERGQzGHyIiIjIZjD4EBERkc1g8CEiIiKbweBDRERENoPBh4iIiGwGgw8RERHZDHupC7A0QggAgFarlbgSIiIiMpX+e1v/PV4XBp/bFBQUAACCgoIkroSIiIgaqqCgABqNps7nZeJO0cjG6HQ6XL16FS4uLpDJZFKX0yJptVoEBQXh8uXLcHV1lbocm8XPwTLwc7AM/BwsQ3N+DkIIFBQUwN/fH3J53SN52ONzG7lcjsDAQKnLsAqurq78BWMB+DlYBn4OloGfg2Vors+hvp4ePQ5uJiIiIpvB4ENEREQ2g8GHmpxKpcLcuXOhUqmkLsWm8XOwDPwcLAM/B8tgCZ8DBzcTERGRzWCPDxEREdkMBh8iIiKyGQw+REREZDMYfIiIiMhmMPhQo6xatQqhoaFQq9Xo0aMHjhw5UmfbDRs2QCaTGW1qtdqM1Vqfffv2Yfjw4fD394dMJsO2bdvueMyePXvQtWtXqFQqtG7dGhs2bGj2Oq1dQz+HPXv21Pi3IJPJkJGRYZ6CrdSiRYvQrVs3uLi4wNvbGyNGjEBiYuIdj/vyyy/Rtm1bqNVqdOzYET/88IMZqrVejfkcpPh+YPChBtuyZQtee+01zJ07F8ePH0fnzp0RGxuLrKysOo9xdXVFenq6YUtJSTFjxdanqKgInTt3xqpVq0xqf+nSJTz88MMYMGAAEhIS8Oqrr+L555/Hjh07mrlS69bQz0EvMTHR6N+Dt7d3M1VoG/bu3YspU6bg119/RVxcHCoqKjB48GAUFRXVecyhQ4fw9NNPY8KECThx4gRGjBiBESNG4NSpU2as3Lo05nMAJPh+EEQN1L17dzFlyhTD46qqKuHv7y8WLVpUa/v169cLjUZjpupsDwDxzTff1NtmxowZon379kb7nnrqKREbG9uMldkWUz6H3bt3CwAiNzfXLDXZqqysLAFA7N27t842Tz75pHj44YeN9vXo0UNMnjy5ucuzGaZ8DlJ8P7DHhxqkvLwcv/32GwYNGmTYJ5fLMWjQIMTHx9d5XGFhIUJCQhAUFIRHH30Up0+fNke5dEN8fLzRZwYAsbGx9X5m1Hy6dOkCPz8/PPjggzh48KDU5Vid/Px8AICHh0edbfhvovmZ8jkA5v9+YPChBrl+/Tqqqqrg4+NjtN/Hx6fOcQqRkZFYt24dvv32W2zatAk6nQ69evXClStXzFEyAcjIyKj1M9NqtSgpKZGoKtvj5+eH1atX4+uvv8bXX3+NoKAg9O/fH8ePH5e6NKuh0+nw6quv4v7770eHDh3qbFfXvwmOt2oapn4OUnw/8O7s1Ox69uyJnj17Gh736tUL7dq1w5o1a7Bw4UIJKyMyr8jISERGRhoe9+rVCxcuXMDKlSvxn//8R8LKrMeUKVNw6tQpHDhwQOpSbJqpn4MU3w/s8aEG8fT0hJ2dHTIzM432Z2ZmwtfX16RzKBQKREdHIykpqTlKpFr4+vrW+pm5urrCwcFBoqoIALp3785/C03kpZdewvbt27F7924EBgbW27aufxOm/h6jujXkc7idOb4fGHyoQZRKJWJiYvDzzz8b9ul0Ovz8889Gqb0+VVVVOHnyJPz8/JqrTLpNz549jT4zAIiLizP5M6Pmk5CQwH8Ld0kIgZdeegnffPMNfvnlF4SFhd3xGP6baHqN+RxuZ5bvB7MOpSarsHnzZqFSqcSGDRvEmTNnxKRJk4Sbm5vIyMgQQggxZswYMXPmTEP7+fPnix07dogLFy6I3377TYwePVqo1Wpx+vRpqd5Ci1dQUCBOnDghTpw4IQCIFStWiBMnToiUlBQhhBAzZ84UY8aMMbS/ePGicHR0FNOnTxdnz54Vq1atEnZ2duKnn36S6i1YhYZ+DitXrhTbtm0Tf/75pzh58qR45ZVXhFwuF7t27ZLqLViFF198UWg0GrFnzx6Rnp5u2IqLiw1tbv+9dPDgQWFvby+WL18uzp49K+bOnSsUCoU4efKkFG/BKjTmc5Di+4HBhxrl/fffF8HBwUKpVIru3buLX3/91fBcv379xNixYw2PX331VUNbHx8fMXToUHH8+HEJqrYe+mnRt2/6n/vYsWNFv379ahzTpUsXoVQqRXh4uFi/fr3Z67Y2Df0clixZIiIiIoRarRYeHh6if//+4pdffpGmeCtS22cAwOi/8dt/LwkhxBdffCHatGkjlEqlaN++vfj+++/NW7iVacznIMX3g+xGsURERERWj2N8iIiIyGYw+BAREZHNYPAhIiIim8HgQ0RERDaDwYeIiIhsBoMPERER2QwGHyIiIrIZDD5EZJL+/fvj1VdfNTwODQ3Fu+++K1k9za1v377473//K3UZFmn06NF45513pC6DqFEYfIgs2Lhx4zBixIga+/fs2QOZTIa8vLw6j92wYQPc3NxqfU4mk2Hbtm13VdvRo0cxadIkk9q2tJD03XffITMzE6NHj67x3KJFi2BnZ4dly5ZJUFnNACqFt956C2+//Tby8/MlrYOoMRh8iKhRvLy84OjoKHUZzeJf//oXnnvuOcjlNX9Frlu3DjNmzMC6deskqMw0QghUVlY22/k7dOiAiIgIbNq0qdleg6i5MPgQUQ1FRUV49tln4ezsDD8/v1ova9zaiyOEwLx58xAcHAyVSgV/f3+8/PLLAKp7KFJSUjB16lTIZDLIZDIAQHZ2Np5++mkEBATA0dERHTt2xOeff270Gv3798fLL7+MGTNmwMPDA76+vpg3b55Rm7y8PEyePBk+Pj5Qq9Xo0KEDtm/fbnj+wIED6NOnDxwcHBAUFISXX34ZRUVFdb73a9eu4ZdffsHw4cNrPLd3716UlJRgwYIF0Gq1OHTokNHz8+bNQ5cuXfCf//wHoaGh0Gg0GD16NAoKCgxtCgoK8Ne//hVOTk7w8/PDypUra/Ti/Pvf/8Y999wDtVoNHx8fjBo1CkB1D+DevXvx3nvvGX6WycnJhh7AH3/8ETExMVCpVDhw4ADKysrw8ssvw9vbG2q1Gr1798bRo0cNr6M/bseOHYiOjoaDgwMeeOABZGVl4ccff0S7du3g6uqKv/zlLyguLjZ6r8OHD8fmzZvr/DkSWSoGHyKqYfr06di7dy++/fZb7Ny5E3v27MHx48frbP/1119j5cqVWLNmDf78809s27YNHTt2BABs3boVgYGBWLBgAdLT05Geng4AKC0tRUxMDL7//nucOnUKkyZNwpgxY3DkyBGjc2/cuBFOTk44fPgwli5digULFiAuLg4AoNPp8NBDD+HgwYPYtGkTzpw5g8WLF8POzg4AcOHCBQwZMgQjR47EH3/8gS1btuDAgQN46aWX6nwvBw4cgKOjI9q1a1fjubVr1+Lpp5+GQqHA008/jbVr19Zoc+HCBWzbtg3bt2/H9u3bsXfvXixevNjw/GuvvYaDBw/iu+++Q1xcHPbv32/0sz127BhefvllLFiwAImJifjpp5/Qt29fAMB7772Hnj17YuLEiYafZVBQkOHYmTNnYvHixTh79iw6deqEGTNm4Ouvv8bGjRtx/PhxtG7dGrGxscjJyTGqed68efjggw9w6NAhXL58GU8++STeffdd/Pe//8X333+PnTt34v333zc6pnv37jhy5AjKysrq/FkSWaRmvQUqEd2VsWPHCjs7O+Hk5GS0qdVqAUDk5ubWeez69esFgBrHOjk5CQDim2++qfW4goICoVQqxRdffGHYl52dLRwcHMQrr7xi2BcSEiJWrlwphBDinXfeEW3atBHl5eW1nvPWtvV5+OGHxbRp0wyP+/XrJ3r37m3Uplu3buKNN94QQgixY8cOIZfLRWJiYq3nmzBhgpg0aZLRvv379wu5XC5KSkpqPWblypUiPDy8xv78/Hzh4OAgEhIShBBCnDhxQjg7O4uCggJDm7lz5wpHR0eh1WoN+6ZPny569OghhBBCq9UKhUIhvvzyS8PzeXl5wtHR0fCz/frrr4Wrq6vROW7Vr18/o89BiJt3id+2bZthX2FhoVAoFOKzzz4z7CsvLxf+/v5i6dKlRsft2rXL0GbRokUCgLhw4YJh3+TJk0VsbKzRa/7+++8CgEhOTq61TiJLxR4fIgs3YMAAJCQkGG2ffPKJURtnZ2fD9sILLxj2u7i41Dg2ISGh3te7cOECysvL0aNHD8M+Dw8PREZG1nnME088gZKSEoSHh2PixIn45ptv7jjGpKqqCgsXLkTHjh3h4eEBZ2dn7NixA6mpqUbtOnXqZPTYz88PWVlZAICEhAQEBgaiTZs2tb7G77//jg0bNhj9fGJjY6HT6XDp0qVajykpKYFara6x//PPP0dERAQ6d+4MAOjSpQtCQkKwZcsWo3ahoaFwcXGptd6LFy+ioqIC3bt3Nzyv0WiMfrYPPvggQkJCEB4ejjFjxuCzzz6rcZmpLvfee6/h7xcuXEBFRQXuv/9+wz6FQoHu3bvj7NmzRsfd+jP28fGBo6MjwsPDjfbp34Oeg4MDAJhcG5GlsJe6ACKqn5OTE1q3bm2078qVK0aPbw0zrq6uhr/L5fIaxzaHoKAgJCYmYteuXYiLi8Pf/vY3LFu2DHv37oVCoaj1mGXLluG9997Du+++i44dO8LJyQmvvvoqysvLjdrdfrxMJoNOpwNw88u3LoWFhZg8ebJhvNGtgoODaz3G09MTubm5NfavXbsWp0+fhr39zV+bOp0O69atw4QJE0yq1xQuLi44fvw49uzZg507d2LOnDmYN28ejh49WucsPT0nJyeTX+dWt9Ysk8lMeg/6y2VeXl6Nek0iqbDHh8gKtG7d2rB5e3vf1bkiIiKgUChw+PBhw77c3FycP3++3uMcHBwwfPhw/Otf/8KePXsQHx+PkydPAgCUSiWqqqqM2h88eBCPPvoonnnmGXTu3Bnh4eF3fI3bderUCVeuXKnzuK5du+LMmTNGPx/9plQqaz0mOjoaGRkZRuHn5MmTOHbsGPbs2WPUc6Z/n+fOnTOp3vDwcCgUCqMBxvn5+TXqt7e3x6BBg7B06VL88ccfSE5Oxi+//AKg9p9lbSIiIqBUKnHw4EHDvoqKChw9ehRRUVEm1VufU6dOITAwEJ6ennd9LiJzYo8PERlxdnbGhAkTMH36dLRq1Qre3t6YPXt2rVO79TZs2ICqqir06NEDjo6O2LRpExwcHBASEgKg+vLPvn37MHr0aKhUKnh6euKee+7BV199hUOHDsHd3R0rVqxAZmZmg76U+/Xrh759+2LkyJFYsWIFWrdujXPnzkEmk2HIkCF44403cN999+Gll17C888/DycnJ5w5cwZxcXH44IMPaj1ndHQ0PD09cfDgQQwbNgxAdW9P9+7dDYOMb9WtWzesXbvWpHV9XFxcMHbsWEyfPh0eHh7w9vbG3LlzIZfLDbPdtm/fjosXL6Jv375wd3fHDz/8AJ1OZ7gcFhoaisOHDyM5ORnOzs7w8PCo9bWcnJzw4osvGl4rODgYS5cuRXFxsVEPVWPt378fgwcPvuvzEJkbe3yIqIZly5ahT58+GD58OAYNGoTevXsjJiamzvZubm74+OOPcf/996NTp07YtWsX/ve//6FVq1YAgAULFiA5ORkRERGGSyNvvfUWunbtitjYWPTv3x++vr61LtZ4J19//TW6deuGp59+GlFRUZgxY4ahR6RTp07Yu3cvzp8/jz59+iA6Ohpz5syBv79/neezs7PDc889h88++wwAUF5ejk2bNmHkyJG1th85ciQ+/fRTVFRUmFTvihUr0LNnTwwbNgyDBg3C/fffj3bt2hnGFbm5uWHr1q144IEH0K5dO6xevRqff/452rdvDwB4/fXXYWdnh6ioKHh5edUYE3WrxYsXY+TIkRgzZgy6du2KpKQk7NixA+7u7ibVWpfS0lJs27YNEydOvKvzEElBJoQQUhdBRGRJMjIy0L59exw/ftzQa9VcioqKEBAQgHfeeadJemLM4cMPP8Q333yDnTt3Sl0KUYOxx4eI6Da+vr5Yu3Ztvb0pjXXixAl8/vnnuHDhAo4fP46//vWvAIBHH320yV+ruSgUihrr+hC1FOzxISIyoxMnTuD5559HYmIilEolYmJisGLFCsOCj0TUvBh8iIiIyGbwUhcRERHZDAYfIiIishkMPkRERGQzGHyIiIjIZjD4EBERkc1g8CEiIiKbweBDRERENoPBh4iIiGwGgw8RERHZjP8HljiKO+uuTBMAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.lineplot(x=r, y=energy)\n",
    "plt.xlabel(\"H-H distance (Angstrom)\")\n",
    "plt.ylabel(\"Energy (eV)\")\n",
    "plt.title(\"H2 Dissociation Curve\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4e4e8fc",
   "metadata": {},
   "source": [
    "## Molecular Orbitals of Benzene"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "efa3bc23",
   "metadata": {},
   "outputs": [],
   "source": [
    "def benzene_geom():\n",
    "    # Structure of benzene ring calculated with DFT using B3LYP functional\n",
    "    # and 6-31G* basis set <https://cccbdb.nist.gov/>\n",
    "    return [\n",
    "        [\"C\", ( 0.0000000,  1.3966140, 0.0000000)],\n",
    "        [\"C\", ( 1.2095030,  0.6983070, 0.0000000)],\n",
    "        [\"C\", ( 1.2095030, -0.6983070, 0.0000000)],\n",
    "        [\"C\", ( 0.0000000, -1.3966140, 0.0000000)],\n",
    "        [\"C\", (-1.2095030, -0.6983070, 0.0000000)],\n",
    "        [\"C\", (-1.2095030,  0.6983070, 0.0000000)],\n",
    "        [\"H\", ( 0.0000000,  2.4836130, 0.0000000)],\n",
    "        [\"H\", ( 2.1508720,  1.2418060, 0.0000000)],\n",
    "        [\"H\", ( 2.1508720, -1.2418060, 0.0000000)],\n",
    "        [\"H\", ( 0.0000000, -2.4836130, 0.0000000)],\n",
    "        [\"H\", (-2.1508720, -1.2418060, 0.0000000)],\n",
    "        [\"H\", (-2.1508720,  1.2418060, 0.0000000)],\n",
    "    ]\n",
    "\n",
    "opts, _ = nanoDFT_options(backend=\"ipu\", float32=True, basis=\"sto-3g\")\n",
    "mol = build_mol(benzene_geom(), opts.basis)\n",
    "energies, _dbginfo = nanoDFT(mol, opts)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "25acb1b7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlQAAAHHCAYAAAB5gsZZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAByVElEQVR4nO3deVhUZf8G8HsGmGHfBEQUEVdEUxGXcFdQSEtNfdPCcguzRHMpl7efu2W5lKaZmqZZ9mZmlsv7orjmrqlYKpB7iiIqwrAIDMzz+wPn5Mg2AzOMg/fnuubCOec553xnYOD2Oc95jkwIIUBERERE5SY3dwFERERElo6BioiIiKiCGKiIiIiIKoiBioiIiKiCGKiIiIiIKoiBioiIiKiCGKiIiIiIKoiBioiIiKiCGKiIiIiIKoiBip5ZQ4cORZ06dXSWZWZm4s0334S3tzdkMhnGjRtnltoqYt26dZDJZLh27Zq5SzGa4r4v165dg0wmw7p16wze3/79+yGTybB///4y23bp0gVdunQx+BhE9GxhoCKL8eeff2LAgAHw8/ODra0tatasie7du2Pp0qVGO8ZHH32EdevW4e2338a3336L119/3Wj7fpplZ2dj5syZegUMY1i+fLlBQehZ/b4QkeWQ8V5+ZAmOHDmCrl27onbt2hgyZAi8vb1x48YNHDt2DJcvX8alS5cM3ufQoUOxf/9+nZ6c559/HtbW1jh06JARq69cBQUFUKvVUCqVkMlkem1z7949eHp6YsaMGZg5c6ZpCwTQtGlTeHh46B3givu+CCGQm5sLGxsbWFlZGXT8/fv3o2vXrti3b1+ZvU/a9ZUVNonIMlmbuwAifXz44YdwcXHByZMn4erqqrMuJSXFaMdJSUlBYGCg0fZnDlZWVgYHjKddcd8XmUwGW1tbM1X0dNNoNMjLy+P7U0FCCOTk5MDOzs7cpZAF4Ck/sgiXL19GkyZNioQpAPDy8iqy7LvvvkNwcDDs7Ozg7u6OQYMG4caNGyXuXzum5urVq9ixYwdkMple45C+++47tGnTBvb29nBzc0OnTp2wa9cunTbLly9HkyZNoFQq4ePjg9GjRyMtLU2nTZcuXdC0aVNcuHABXbt2hb29PWrWrIn58+cXOebSpUvRpEkT6ZitWrXC999/L60vbgzV77//jvDwcHh4eMDOzg7+/v4YPnw4AODatWvw9PQEAMyaNUt67Y/3VCUkJGDAgAFwd3eHra0tWrVqha1bt+rUpT3u4cOHMWHCBHh6esLBwQEvv/wy7t69K7WrU6cOzp8/jwMHDkjHKqmXqLTvS0ljqPSptSSrVq1CvXr1YGdnhzZt2uDgwYN6baelz8+dId/r3NxczJgxA/Xr14dSqYSvry8mTZqE3NxcnXYymQzR0dHYsGGD9LMWExMDAPjjjz/QuXNn2NnZoVatWpg7dy7Wrl2r8zMyZMgQeHh4QK1WF6mhR48eaNSoUZmv/fjx4+jZsyfc3Nzg4OCAZs2aYcmSJTpt9u7di44dO8LBwQGurq7o06cP4uPjddrMnDkTMpkMly5dwtChQ+Hq6goXFxcMGzYM2dnZUrumTZuia9euRerQaDSoWbMmBgwYoLNs8eLFaNKkCWxtbVG9enW89dZbePDggc62derUwYsvvoidO3eiVatWsLOzw8qVKwEA169fR+/eveHg4AAvLy+MHz8eO3fuLHYs3vHjxxEREQEXFxfY29ujc+fOOHz4cLlep5Y+v2v+97//Se+vk5MTevXqhfPnzxfZF5mIILIAPXr0EE5OTuLPP/8ss+3cuXOFTCYTAwcOFMuXLxezZs0SHh4eok6dOuLBgwdSuyFDhgg/Pz8hhBDJycni22+/FR4eHqJFixbi22+/Fd9++63IzMws8TgzZ84UAES7du3EggULxJIlS8Rrr70mJk+eLLWZMWOGACDCwsLE0qVLRXR0tLCyshKtW7cWeXl5UrvOnTsLHx8f4evrK959912xfPly0a1bNwFA/Pe//5XarVq1SgAQAwYMECtXrhRLliwRI0aMEGPHjpXarF27VgAQV69eFUIIcefOHeHm5iYaNmwoFixYIL766ivxwQcfiMaNGwshhMjMzBRffvmlACBefvll6bWfPXtWCCHEuXPnhIuLiwgMDBSffPKJWLZsmejUqZOQyWTi559/LnLcoKAg0a1bN7F06VIxceJEYWVlJV555RWp3ZYtW0StWrVEQECAdKxdu3YV+x6X9n25evWqACDWrl0rtde31n379gkAYt++fdKy1atXS9/Pzz//XIwbN064urqKunXris6dO5f4c6Cl78+dvt/rgoIC0aNHD2Fvby/GjRsnVq5cKaKjo4W1tbXo06ePzrEBiMaNGwtPT08xa9Ys8cUXX4gzZ86ImzdvCnd3d1GtWjUxa9YssXDhQhEQECCaN2+u8zMSGxsrAIht27bp7Pf27dvCyspKzJ49u9TXvmvXLqFQKISfn5+YMWOG+PLLL8XYsWNFWFiY1CY2NlZYW1uLhg0bivnz50vvj5ubm1SHEP98ZoKCgkS/fv3E8uXLxZtvvikAiEmTJkntZs+eLeRyubh9+7ZOLQcOHBAAxKZNm6Rlb775prC2thZRUVFixYoVYvLkycLBwaHI59DPz0/Ur19fuLm5iSlTpogVK1aIffv2iczMTFG3bl1hZ2cnpkyZIhYvXizatGkjvY+P/xzt2bNHKBQKERISIhYtWiQ+++wz0axZM6FQKMTx48cNfp1C6Pe7Zv369UImk4mIiAixdOlS8cknn4g6deoIV1dXnfeXTIeBiizCrl27hJWVlbCyshIhISFi0qRJYufOnTq/DIUQ4tq1a8LKykp8+OGHOsv//PNPYW1trbP88UCl5efnJ3r16lVmPRcvXhRyuVy8/PLLoqCgQGedRqMRQgiRkpIiFAqF6NGjh06bZcuWCQDi66+/lpZ17txZABDr16+XluXm5gpvb2/Rv39/aVmfPn1EkyZNSq3tyUC1ZcsWAUCcPHmyxG3u3r0rAIgZM2YUWRcaGiqee+45kZOTo/Ma27VrJxo0aFDkuGFhYdJ7IIQQ48ePF1ZWViItLU1a1qRJE71CilZx35fiApW+tT4ZqPLy8oSXl5do0aKFyM3NldppA2xZtRryc6fv9/rbb78VcrlcHDx4UGefK1asEADE4cOHpWUAhFwuF+fPn9dpO2bMGCGTycSZM2ekZffv3xfu7u46PyMFBQWiVq1aYuDAgTrbf/rpp0Imk4krV66U+Nrz8/OFv7+/8PPz0wmOQgidn4MWLVoILy8vcf/+fWnZ2bNnhVwuF2+88Ya0TBs0hg8frrOvl19+WVSrVk16npiYKACIpUuX6rR75513hKOjo8jOzhZCCHHw4EEBQGzYsEGnXUxMTJHlfn5+AoCIiYnRabto0SIBQPzyyy/SsocPH4qAgACdnyONRiMaNGggwsPDdV57dna28Pf3F927dzf4derzuyYjI0O4urqKqKgonfXJycnCxcWlyHIyDZ7yI4vQvXt3HD16FL1798bZs2cxf/58hIeHo2bNmjqnc37++WdoNBq88soruHfvnvTw9vZGgwYNsG/fPqPU88svv0Cj0WD69OmQy3U/RtqB4Lt370ZeXh7GjRun0yYqKgrOzs7YsWOHznaOjo4YPHiw9FyhUKBNmza4cuWKtMzV1RU3b97EyZMn9a5Ve5p0+/btxZ7SKU1qair27t2LV155BRkZGdL7ef/+fYSHh+PixYtISkrS2WbkyJE6g+E7duyIgoICXL9+3aBjG6o8tWr9/vvvSElJwahRo6BQKKTlQ4cOhYuLS5nHNvTnTp/v9aZNm9C4cWMEBATo7LNbt24AUGSfnTt3LjLOLCYmBiEhIWjRooW0zN3dHZGRkTrt5HI5IiMjsXXrVmRkZEjLN2zYgHbt2sHf37/E137mzBlcvXoV48aNK3JKXvtzcPv2bcTFxWHo0KFwd3eX1jdr1gzdu3fHf//73yL7HTVqlM7zjh074v79+1CpVACAhg0bokWLFti4caPUpqCgAD/99BNeeukladzTpk2b4OLigu7du+u8j8HBwXB0dCzyPvr7+yM8PFxnWUxMDGrWrInevXtLy2xtbREVFaXTLi4uDhcvXsRrr72G+/fvS8fKyspCaGgofvvtN2g0GoNepz6/a2JjY5GWloZXX31V5zVaWVmhbdu2Rvu9R6VjoHrK7dixA23btoWdnR3c3NzQt2/fEtuq1WpMnjwZzz33HBwcHODj44M33ngDt27d0mnXu3dv1K5dG7a2tqhRowZef/31Im127tyJ559/Hk5OTvD09ET//v0NntforbfeksajeHp6ok+fPkhISDBoH49r3bo1fv75Zzx48AAnTpzA1KlTkZGRgQEDBuDChQsAgIsXL0IIgQYNGsDT01PnER8fb7QB7JcvX4ZcLi91ALs2QDw5/kShUKBu3bpFAkatWrWKXJXn5uamM85j8uTJcHR0RJs2bdCgQQOMHj26yNiMJ3Xu3Bn9+/fHrFmz4OHhgT59+mDt2rVFxuEU59KlSxBCYNq0aUXezxkzZgAoelFA7dq1i7wGAEXGqxhbeWrV0n4vGjRooLPcxsYGdevWLfPYhv7c6fO9vnjxIs6fP19kfw0bNiz2tRQXeq5fv4769esXWV7csjfeeAMPHz7Eli1bAACJiYk4depUmVNUXL58GUDhmKaSlPRZAIDGjRtLoeNx+vwcDRw4EIcPH5aC8v79+5GSkoKBAwdKbS5evIj09HR4eXkVeS8zMzP1fh/r1atX5Hv25Pt48eJFAIVj0p481urVq5Gbm4v09HSDXqc+v2u0x+3WrVuR4+7atcuoF+5QyXiV31Ns8+bNiIqKwkcffYRu3bohPz8f586dK7F9dnY2Tp8+jWnTpqF58+Z48OAB3n33XfTu3Ru///671K5r167497//jRo1aiApKQnvvfceBgwYgCNHjgAArl69ij59+mDChAnYsGED0tPTMX78ePTr1w+nT5/Wu/7g4GBERkaidu3aSE1NxcyZM9GjRw9cvXq1QlehKRQKtG7dGq1bt0bDhg0xbNgwbNq0CTNmzIBGo4FMJsP//ve/Yo/h6OhY7uOaWknviXhsZpPGjRsjMTER27dvR0xMDDZv3ozly5dj+vTpmDVrVrHby2Qy/PTTTzh27Bi2bduGnTt3Yvjw4Vi0aBGOHTtW6nui/d/0e++9V+R/7VpP/lHR53WYQnlqNeaxDfm50+c90mg0eO655/Dpp58W29bX11fneUWvRAsMDERwcDC+++47vPHGG/juu++gUCjwyiuvVGi/5aXPezRw4EBMnToVmzZtwrhx4/Djjz/CxcUFERERUhuNRgMvLy9s2LCh2P1pL8jQqsj7qP0ZXLBggU6v4OPK87Og73G//fZbeHt7F1lvbc0/9ZWB7/JTKj8/H++++y4WLFiAESNGSMtL+1+Ki4sLYmNjdZYtW7YMbdq0wd9//y39T2j8+PHSej8/P0yZMgV9+/aFWq2GjY0NTp06hYKCAsydO1fqYn7vvffQp08fqQ0A/Prrr5g1axYuXLgAHx8fDBkyBB988IH04R05cqR0nDp16mDu3Llo3rw5rl27hnr16lXwHSrUqlUrAIWnFACgXr16EELA399f+p+8KdSrVw8ajQYXLlwo8Renn58fgML/6T/ey5GXl4erV68iLCysXMd2cHDAwIEDMXDgQOTl5aFfv3748MMPMXXq1FIvk3/++efx/PPP48MPP8T333+PyMhI/PDDD3jzzTdLnK9KW7eNjU256y2OvvNjGaIitWq/VxcvXpROqQGFvb5Xr15F8+bNS93eFD939erVw9mzZxEaGlru98vPz6/YOdpKmrftjTfewIQJE3D79m18//336NWrl9RjUlqdAHDu3LkS3/fHPwtPSkhIgIeHBxwcHEo9TnH8/f3Rpk0bbNy4EdHR0fj555/Rt29fKJVKnfp2796N9u3blzss+fn54cKFCxBC6Hwvnnwfte+Fs7Oz0T4v+vyu0R7Xy8vLqJ9TMgxP+T2lTp8+jaSkJMjlcgQFBaFGjRp44YUXSu2hKk56ejpkMlmx0w0AheNOtOMktEEpODgYcrkca9euRUFBAdLT0/Htt98iLCxManPw4EG88cYbePfdd3HhwgWsXLkS69atw4cffljscbKysrB27Vr4+/sX+Z+1Pvbt21fs/9i0Yy+0pxL69esHKysrzJo1q0h7IQTu379v8LGL07dvX8jlcsyePbvImAjtccPCwqBQKPD555/r1LJmzRqkp6ejV69eBh/3yfoVCgUCAwMhhChxfNSDBw+KvBfaX8za03729vYAUGQ6By8vL3Tp0gUrV66UQuvjHp8OwRAODg5FjlVRFam1VatW8PT0xIoVK5CXlyctX7dunV51muLn7pVXXkFSUhK++uqrIusePnxY5BRZccLDw3H06FHExcVJy7Sf+eK8+uqrkMlkePfdd3HlyhWdcV4ladmyJfz9/bF48eIi75X2vahRowZatGiBb775RqfNuXPnsGvXLvTs2bPM45Rk4MCBOHbsGL7++mvcu3dP53QfUPg+FhQUYM6cOUW2zc/P1+v7Gx4ejqSkJJ3xmjk5OUW+N8HBwahXrx4WLlyIzMzMIvspz+dFn9814eHhcHZ2xkcffVTs74Hyfk7JMOyhekppB6fOnDkTn376KerUqYNFixahS5cu+Ouvv3QGdpYkJycHkydPxquvvgpnZ2eddZMnT8ayZcuQnZ2N559/Htu3b5fW+fv7Y9euXXjllVfw1ltvoaCgACEhIToDR2fNmoUpU6ZgyJAhAAp7B+bMmYNJkyZJ41WAwjmYJk2ahKysLDRq1AixsbE6g371NWbMGGRnZ+Pll19GQEAA8vLycOTIEWzcuBF16tTBsGHDABT+T23u3LmYOnUqrl27hr59+8LJyQlXr17Fli1bMHLkSLz33nsGH/9J9evXxwcffIA5c+agY8eO6NevH5RKJU6ePAkfHx/MmzcPnp6emDp1KmbNmoWIiAj07t0biYmJWL58OVq3bq3XH6sn9ejRA97e3mjfvj2qV6+O+Ph4LFu2DL169YKTk1Ox23zzzTdYvnw5Xn75ZdSrVw8ZGRn46quv4OzsLP0hs7OzQ2BgIDZu3IiGDRvC3d0dTZs2RdOmTfHFF1+gQ4cOeO655xAVFYW6devizp07OHr0KG7evImzZ88a/DqCg4Px5ZdfYu7cuahfvz68vLx0eobKq7y12tjYYO7cuXjrrbfQrVs3DBw4EFevXsXatWv1GkNlip+7119/HT/++CNGjRqFffv2oX379igoKEBCQgJ+/PFHaa6k0kyaNAnfffcdunfvjjFjxsDBwQGrV6+WTsM/2fPl6emJiIgIbNq0Ca6urnqFfrlcji+//BIvvfQSWrRogWHDhqFGjRpISEjA+fPnsXPnTgCFp8FeeOEFhISEYMSIEXj48CGWLl0KFxeXCs3O/8orr+C9997De++9B3d39yI9NJ07d8Zbb72FefPmIS4uDj169ICNjQ0uXryITZs2YcmSJTpzVhXnrbfewrJly/Dqq6/i3XffRY0aNbBhwwapR1j7PsrlcqxevRovvPACmjRpgmHDhqFmzZpISkrCvn374OzsjG3bthn0+vT5XePs7Iwvv/wSr7/+Olq2bIlBgwbB09MTf//9N3bs2IH27dtj2bJlBh2XyqESrygkIcTkyZMFgFIf8fHxYsOGDQKAWLlypbRtTk6O8PDwECtWrCjzOHl5eeKll14SQUFBIj09vcj6u3fvisTERLFr1y7Rvn170bNnT+kS3Nu3b4sGDRqI999/X5w+fVocOHBAdO7cWYSGhkptPDw8hK2trXBwcJAetra2AoDIysqSjpOWlib++usvceDAAfHSSy+Jli1biocPHxr8vv3vf/8Tw4cPFwEBAcLR0VEoFApRv359MWbMGHHnzp0i7Tdv3iw6dOgg1RYQECBGjx4tEhMTpTYVmTZB6+uvvxZBQUFCqVQKNzc30blzZxEbG6vTZtmyZSIgIEDY2NiI6tWri7fffrvI5eWdO3cudjqEJ2tcuXKl6NSpk6hWrZpQKpWiXr164v3339f5Hj85bcLp06fFq6++KmrXri2USqXw8vISL774ovj99991jnXkyBERHBwsFApFkSkULl++LN544w3h7e0tbGxsRM2aNcWLL74ofvrppyLHfXJ6huLmfEpOTha9evUSTk5Oek1LoO+0CfrWWlxNQgixfPly4e/vL5RKpWjVqpX47bffROfOnfWe4kGfnzt9v9dCFH6OP/nkE9GkSRPpZyw4OFjMmjVL53sOQIwePbrYms6cOSM6duwolEqlqFWrlpg3b574/PPPBQCRnJxcpP2PP/4oAIiRI0fq9Zq1Dh06JLp37y6cnJyEg4ODaNasWZEpDXbv3i3at28v7OzshLOzs3jppZfEhQsXdNpopxO4e/euzvInf64f1759ewFAvPnmmyXWt2rVKhEcHCzs7OyEk5OTeO6558SkSZPErVu3pDalff6vXLkievXqJezs7ISnp6eYOHGi2Lx5swAgjh07ptP2zJkzol+/ftLn1M/PT7zyyitiz5495X6d+vyu2bdvnwgPDxcuLi7C1tZW1KtXTwwdOrTIZ51Mg/fyq2R3794ts/u/bt26OHz4MLp164aDBw+iQ4cO0rq2bdsiLCysxFNrQOG4j1deeQVXrlzB3r17Ua1atVKPd/PmTfj6+uLIkSMICQnBtGnTEBMTo3NpvrbN0aNH8fzzz8POzg6zZs1Cv379iq3/yct7gcKxQ25ubli9ejVeffXVUmsiItMZN24cVq5ciczMzCKDon/99Vf07dsXv/32Gzp27GimCi3D4sWLMX78eNy8eRM1a9Y0dzlkZjzlV8m0l7KWJTg4GEqlEomJiVKgUqvVuHbtmjTAszjaMHXx4kXs27evzDAF/HOFiHY8TXZ2dpFApP2lq23bsmVLJCYmGnTVlCicSFavy/WJyDgePnyoMxj7/v37+Pbbb9GhQ4dirzD76quvULduXZ3/yFHR9zEnJwcrV65EgwYNGKYIAAPVU8vZ2RmjRo3CjBkz4OvrCz8/PyxYsAAA8K9//UtqFxAQgHnz5uHll1+GWq3GgAEDcPr0aWzfvh0FBQVITk4GUDiZn0KhwPHjx3Hy5El06NABbm5uuHz5MqZNm4Z69eohJCQEANCrVy989tlnmD17Nl599VVkZGTg3//+N/z8/BAUFAQAmD59Ol588UXUrl0bAwYMgFwux9mzZ3Hu3DnMnTsXV65cwcaNG9GjRw94enri5s2b+Pjjj2FnZ1ehAahEZJiQkBB06dIFjRs3xp07d7BmzRqoVCpMmzZNp90PP/yAP/74Azt27MCSJUtMciWmJevXrx9q166NFi1aID09Hd999x0SEhJKHOBPzyDznnGk0uTl5YmJEycKLy8v4eTkJMLCwsS5c+d02uCxMSTaMSXFPbRjRf744w/RtWtX4e7uLpRKpahTp44YNWqUuHnzps5+//Of/4igoCDh4OAgPD09Re/evUV8fLxOm5iYGNGuXTtpPESbNm3EqlWrhBBCJCUliRdeeEF4eXkJGxsbUatWLfHaa6+JhIQE07xZRFSsqVOnigYNGgg7Ozthb28vOnToUGTsjRCFv0scHR3FiBEjhFqtNkOlT7fPPvtMNGnSRBov2rJlS/HDDz+Yuyx6inAMFREREVEFcR4qIiIiogpioCIiIiKqIA5KrwQajQa3bt2Ck5MTB3oSERFZCCEEMjIy4OPjU+x0QE82tgj3798Xr732mnBychIuLi5i+PDhIiMjo9RtHj58KN555x3h7u4uHBwcRL9+/YpMZIdiBnD/5z//0Wmzb98+ERQUJBQKhahXr16RiQTLcuPGjTIn8+SDDz744IMPPp7Ox40bN8r8W28xg9JfeOEF3L59GytXroRarcawYcPQunVrfP/99yVu8/bbb2PHjh1Yt24dXFxcEB0dDblcjsOHD0ttZDIZ1q5dq3N3cldXV+mWAlevXkXTpk0xatQovPnmm9izZw/GjRuHHTt2lHhH+yelp6fD1dUVN27cKHILGCIiIno6qVQq+Pr6Ii0tDS4uLqW2tYhAFR8fj8DAQJw8eVK6d1VMTAx69uyJmzdvwsfHp8g26enp8PT0xPfffy/dpykhIQGNGzeWZvsGCgPVli1b0Ldv32KPPXnyZOzYsUPnpsSDBg1CWloaYmJi9KpfpVLBxcUF6enpDFREREQWwpC/3xYxKP3o0aNwdXXVuRFoWFgY5HI5jh8/Xuw2p06dglqt1rlRZkBAAGrXro2jR4/qtB09ejQ8PDzQpk0bfP311zp3iz969GiRm21q7+BektzcXKhUKp0HERERVV0WMSg9OTkZXl5eOsusra3h7u4uzQRe3DYKhQKurq46y6tXr66zzezZs9GtWzfY29tj165deOedd5CZmYmxY8dK+6levXqRfahUqiK3ItCaN28eZs2aVZ6XSkRERBbIrD1UU6ZMgUwmK/WRkJBg0hqmTZuG9u3bIygoCJMnT8akSZOkW7yU19SpU5Geni49bty4YaRqiYiI6Glk1h6qiRMnYujQoaW2qVu3Lry9vZGSkqKzPD8/H6mpqfD29i52O29vb+Tl5SEtLU2nl+rOnTslbgMAbdu2xZw5c5CbmwulUglvb2/cuXNHp82dO3fg7OxcbO8UACiVSiiVylJfFxEREVUdZg1Unp6e8PT0LLNdSEgI0tLScOrUKQQHBwMA9u7dC41Gg7Zt2xa7TXBwMGxsbLBnzx70798fAJCYmIi///5buglwceLi4uDm5iYFopCQEPz3v//VaRMbG1vqPoiIiOjZYhFjqBo3boyIiAhERUVhxYoVUKvViI6OxqBBg6Qr/JKSkhAaGor169ejTZs2cHFxwYgRIzBhwgS4u7vD2dkZY8aMQUhIiHSF37Zt23Dnzh08//zzsLW1RWxsLD766CO899570rFHjRqFZcuWYdKkSRg+fDj27t2LH3/8ETt27DDLe0FERERPH4sIVACwYcMGREdHIzQ0FHK5HP3798fnn38urVer1UhMTER2dra07LPPPpPa5ubmIjw8HMuXL5fW29jY4IsvvsD48eMhhED9+vXx6aefIioqSmrj7++PHTt2YPz48ViyZAlq1aqF1atX6z0HFREREVV9FjEPlaXjPFRERESWp8rNQ0VERET0NGOgIiIiIqogBioiIiKiCmKgIiIiIqogBioLptEI3MvMxaWUTHOXQkRE9ExjoLJgf6dmo9Xc3Xhp6SFzl0JERPRMY6CyYJ5OhbO5P1QXICs338zVEBERPbsYqCyYg9Ia9gorAMDdjFwzV0NERPTsYqCycB6Ohb1UdzMZqIiIiMyFgcrCaU/7sYeKiIjIfBioLJynIwMVERGRuTFQWTj2UBEREZkfA5WFY6AiIiIyPwYqC6cNVPc4KJ2IiMhsGKgsnCev8iMiIjI7BioL58FTfkRERGbHQGXhHj/lp9EIM1dDRET0bGKgsnAejgoAgLpAIP2h2szVEBERPZsYqCyc0toKLnY2ADiOioiIyFwYqKoATp1ARERkXgxUVYD2Sj9OnUBERGQeDFRVAHuoiIiIzIuBqgpgoCIiIjIvBqoqwIM3SCYiIjIrBqoqQOqh4hgqIiIis2CgqgJ4yo+IiMi8GKiqAE+e8iMiIjIrBqoqQNtDlZqdh/wCjZmrISIievYwUFUB7g4KyGWAEEBqVp65yyEiInrmMFBVAVZyGao9Ou2XwtN+RERElc5iAlVqaioiIyPh7OwMV1dXjBgxApmZmaVuk5OTg9GjR6NatWpwdHRE//79cefOHWn9unXrIJPJin2kpKQAAPbv31/s+uTkZJO+XkNJUyfwSj8iIqJKZzGBKjIyEufPn0dsbCy2b9+O3377DSNHjix1m/Hjx2Pbtm3YtGkTDhw4gFu3bqFfv37S+oEDB+L27ds6j/DwcHTu3BleXl46+0pMTNRp9+R6c+OVfkREROZjbe4C9BEfH4+YmBicPHkSrVq1AgAsXboUPXv2xMKFC+Hj41Nkm/T0dKxZswbff/89unXrBgBYu3YtGjdujGPHjuH555+HnZ0d7OzspG3u3r2LvXv3Ys2aNUX25+XlBVdXV9O8QCPglX5ERETmYxE9VEePHoWrq6sUpgAgLCwMcrkcx48fL3abU6dOQa1WIywsTFoWEBCA2rVr4+jRo8Vus379etjb22PAgAFF1rVo0QI1atRA9+7dcfjw4VLrzc3NhUql0nmYGnuoiIiIzMciAlVycnKRU2zW1tZwd3cvcSxTcnIyFApFkV6l6tWrl7jNmjVr8Nprr+n0WtWoUQMrVqzA5s2bsXnzZvj6+qJLly44ffp0ifXOmzcPLi4u0sPX11fPV1p+2kB1j2OoiIiIKp1ZA9WUKVNKHBSufSQkJFRKLUePHkV8fDxGjBihs7xRo0Z46623EBwcjHbt2uHrr79Gu3bt8Nlnn5W4r6lTpyI9PV163Lhxw9Tls4eKiIjIjMw6hmrixIkYOnRoqW3q1q0Lb29v6ao7rfz8fKSmpsLb27vY7by9vZGXl4e0tDSdXqo7d+4Uu83q1avRokULBAcHl1l3mzZtcOjQoRLXK5VKKJXKMvdjTJ68yo+IiMhszBqoPD094enpWWa7kJAQpKWl4dSpU1Lg2bt3LzQaDdq2bVvsNsHBwbCxscGePXvQv39/AIVX6v39998ICQnRaZuZmYkff/wR8+bN06vuuLg41KhRQ6+2lYU9VEREROZjEVf5NW7cGBEREYiKisKKFSugVqsRHR2NQYMGSVf4JSUlITQ0FOvXr0ebNm3g4uKCESNGYMKECXB3d4ezszPGjBmDkJAQPP/88zr737hxI/Lz8zF48OAix168eDH8/f3RpEkT5OTkYPXq1di7dy927dpVKa9dX9oeqoycfOSoC2BrY2XmioiIiJ4dFhGoAGDDhg2Ijo5GaGgo5HI5+vfvj88//1xar1arkZiYiOzsbGnZZ599JrXNzc1FeHg4li9fXmTfa9asQb9+/YqdFiEvLw8TJ05EUlIS7O3t0axZM+zevRtdu3Y1yessL2c7ayis5Mgr0OBuRi583e3NXRIREdEzQyaEEOYuoqpTqVRwcXFBeno6nJ2dTXac9h/vRVLaQ/z8Tju0rO1msuMQERE9Cwz5+20R0yaQfjy0UydwHBUREVGlYqCqQnilHxERkXkwUFUhvNKPiIjIPBioqhAGKiIiIvNgoKpCPB0VABioiIiIKhsDVRUi9VBxDBUREVGlYqCqQnjKj4iIyDwYqKoQT0dbAMC9zFxwejEiIqLKw0BVhXg4FY6hylFrkJmbb+ZqiIiInh0MVFWIvcIajsrCuwnxtB8REVHlYaCqYjiOioiIqPIxUFUxHtqpE3ilHxERUaVhoKpi2ENFRERU+Rioqhjpfn4MVERERJWGgaqK0fZQ3eMpPyIiokrDQFXF8JQfERFR5WOgqmJ4+xkiIqLKx0BVxWhnS2cPFRERUeVhoKpi/hlDlQeNhrefISIiqgwMVFVMtUfzUBVoBB5k55m5GiIiomcDA1UVY2Mlh5u9DQCOoyIiIqosDFRVkHTaL4M9VERERJWBgaoK+udKvxwzV0JERPRsYKCqgjhbOhERUeVioKqCOLknERFR5WKgqoIYqIiIiCoXA1UV5OHI2dKJiIgqEwNVFcQeKiIiosrFQFUFPT5bOhEREZkeA1UVpL3KLzUrD+oCjZmrISIiqvosJlClpqYiMjISzs7OcHV1xYgRI5CZmVnqNqtWrUKXLl3g7OwMmUyGtLS0cu33jz/+QMeOHWFrawtfX1/Mnz/fmC/N6NzsFbCSywAA99lLRUREZHIWE6giIyNx/vx5xMbGYvv27fjtt98wcuTIUrfJzs5GREQE/v3vf5d7vyqVCj169ICfnx9OnTqFBQsWYObMmVi1apXRXpuxyeUyeDy6px/HUREREZmeTAghzF1EWeLj4xEYGIiTJ0+iVatWAICYmBj07NkTN2/ehI+PT6nb79+/H127dsWDBw/g6upq0H6//PJLfPDBB0hOToZCURhSpkyZgl9++QUJCQl61a9SqeDi4oL09HQ4OzuX4x0w3ItLD+JckgpfD22FbgHVK+WYREREVYkhf78toofq6NGjcHV1lUIPAISFhUEul+P48eMm3e/Ro0fRqVMnKUwBQHh4OBITE/HgwYNyH9vUOFs6ERFR5bE2dwH6SE5OhpeXl84ya2truLu7Izk52aT7TU5Ohr+/v06b6tWrS+vc3NyK7Dc3Nxe5uf8EGZVKVe4ay8uDgYqIiKjSmLWHasqUKZDJZKU+9D2t9jSZN28eXFxcpIevr2+l18CpE4iIiCqPWXuoJk6ciKFDh5bapm7duvD29kZKSorO8vz8fKSmpsLb27vcx9dnv97e3rhz545OG+3zko49depUTJgwQXquUqkqPVRxck8iIqLKY9ZA5enpCU9PzzLbhYSEIC0tDadOnUJwcDAAYO/evdBoNGjbtm25j6/PfkNCQvDBBx9ArVbDxsYGABAbG4tGjRoVe7oPAJRKJZRKZbnrMgYGKiIiospjEYPSGzdujIiICERFReHEiRM4fPgwoqOjMWjQIOkKv6SkJAQEBODEiRPSdsnJyYiLi8OlS5cAAH/++Sfi4uKQmpqq935fe+01KBQKjBgxAufPn8fGjRuxZMkSnR6op5En7+dHRERUaSwiUAHAhg0bEBAQgNDQUPTs2RMdOnTQmQtKrVYjMTER2dnZ0rIVK1YgKCgIUVFRAIBOnTohKCgIW7du1Xu/Li4u2LVrF65evYrg4GBMnDgR06dPL3MOLHNjDxUREVHlsYh5qCydOeahyshR47mZuwAAF2aHw15hERd0EhERPTWq3DxUZDhHpTWU1oXf3nsZvNKPiIjIlBioqiiZTPbPab/MHDNXQ0REVLUxUFVh/4yjYg8VERGRKTFQVWG80o+IiKhyMFBVYbzSj4iIqHIwUFVhDFRERESVg4GqCmOgIiIiqhwMVFWYB8dQERERVQoGqipM20N1jz1UREREJsVAVYU9fpUfJ8QnIiIyHQaqKkzbQ5WXr4EqJ9/M1RAREVVdDFRVmK2NFZxsC+/hx4HpREREpsNAVcXxSj8iIiLTY6Cq4jhbOhERkekxUFVx7KEiIiIyPQaqKk6ai4qBioiIyGQYqKo4aS4qnvIjIiIyGQaqKo6n/IiIiEyPgaqKY6AiIiIyPQaqKo5X+REREZkeA1UV5/Woh+p+Zi4KNLz9DBERkSkwUFVx7g4KyGSARgCpWXnmLoeIiKhKYqCq4qyt5HC3VwDgOCoiIiJTYaB6BnDqBCIiItNioHoG8Eo/IiIi02KgegbwSj8iIiLTYqB6BrCHioiIyLQYqJ4BDFRERESmxUD1DGCgIiIiMi0GqmcAx1ARERGZFgPVM8CD0yYQERGZlMUEqtTUVERGRsLZ2Rmurq4YMWIEMjMzS91m1apV6NKlC5ydnSGTyZCWlqaz/tq1axgxYgT8/f1hZ2eHevXqYcaMGcjLy9NpI5PJijyOHTtmipdpEtoeqrRsNXLzC8xcDRERUdVjbe4C9BUZGYnbt28jNjYWarUaw4YNw8iRI/H999+XuE12djYiIiIQERGBqVOnFlmfkJAAjUaDlStXon79+jh37hyioqKQlZWFhQsX6rTdvXs3mjRpIj2vVq2a8V6cibnY2cDGSgZ1gcD9zDz4uNqZuyQiIqIqxSICVXx8PGJiYnDy5Em0atUKALB06VL07NkTCxcuhI+PT7HbjRs3DgCwf//+Ytdrw5ZW3bp1kZiYiC+//LJIoKpWrRq8vb0r/mLMQC6XwcNRidvpObibkctARUREZGQGB6rc3FwcP34c169fR3Z2Njw9PREUFAR/f39T1AcAOHr0KFxdXaUwBQBhYWGQy+U4fvw4Xn75ZaMdKz09He7u7kWW9+7dGzk5OWjYsCEmTZqE3r17l7iP3Nxc5Ob+M15JpVIZrb7y8nT6J1ARERGRcekdqA4fPowlS5Zg27ZtUKvVcHFxgZ2dHVJTU5Gbm4u6deti5MiRGDVqFJycnIxaZHJyMry8vHQLt7aGu7s7kpOTjXacS5cuYenSpTq9U46Ojli0aBHat28PuVyOzZs3o2/fvvjll19KDFXz5s3DrFmzjFaXMfBKPyIiItPRa1B67969MXDgQNSpUwe7du1CRkYG7t+/j5s3byI7OxsXL17E//3f/2HPnj1o2LAhYmNj9Tr4lClTih3w/fgjISGhQi9QX0lJSYiIiMC//vUvREVFScs9PDwwYcIEtG3bFq1bt8bHH3+MwYMHY8GCBSXua+rUqUhPT5ceN27cqIyXUCrORUVERGQ6evVQ9erVC5s3b4aNjU2x6+vWrYu6detiyJAhuHDhAm7fvq3XwSdOnIihQ4eW2qZu3brw9vZGSkqKzvL8/HykpqYaZVzTrVu30LVrV7Rr1w6rVq0qs33btm1LDY1KpRJKpbLCdRmThyOnTiAiIjIVvQLVW2+9hYIC/S63DwwMRGBgoF5tPT094enpWWa7kJAQpKWl4dSpUwgODgYA7N27FxqNBm3bttXrWCVJSkpC165dERwcjLVr10IuL7vTLi4uDjVq1KjQcSsbe6iIiIhMR+8xVDVr1sTQoUMxfPhwNGzY0JQ1FdG4cWNEREQgKioKK1asgFqtRnR0NAYNGiRd4ZeUlITQ0FCsX78ebdq0AVA49io5ORmXLl0CAPz5559wcnJC7dq14e7ujqSkJHTp0gV+fn5YuHAh7t69Kx1T2/P1zTffQKFQICgoCADw888/4+uvv8bq1asr8y2oMAYqIiIi09F7Ys/Ro0fjp59+QuPGjdGxY0esW7cO2dnZpqxNx4YNGxAQEIDQ0FD07NkTHTp00Dk9p1arkZiYqFPTihUrEBQUJI2J6tSpE4KCgrB161YAQGxsLC5duoQ9e/agVq1aqFGjhvR43Jw5cxAcHIy2bdvi119/xcaNGzFs2LBKeNXGIwUqnvIjIiIyOpkQQhiywf79+7F27Vps3rwZVlZWeOWVV/Dmm29W+NRbVaZSqeDi4oL09HQ4OzubpYZr97LQZeF+2CuscGF2RNkbEBERPeMM+ftt8K1nunTpgm+++QbJyclYtGgR4uPjERISgiZNmuDTTz8td9FkWtoequy8AmTl5pu5GiIioqql3Pfyc3R0xJtvvolDhw5h27ZtSE5Oxvvvv2/M2siIHJTWsFdYAeA4KiIiImMrd6DKzs7GunXr0LlzZ/Tu3RvVqlXDhx9+aMzayMg4dQIREZFpGHzrmSNHjuDrr7/Gpk2bkJ+fjwEDBmDOnDno1KmTKeojI/J0UuLv1Gz2UBERERmZ3oFq/vz5WLt2Lf766y+0atUKCxYswKuvvmr028yQ6fD2M0RERKahd6BasGABBg8ejE2bNqFp06amrIlMhHNRERERmYbegerWrVtFbj2Tk5MDW1tboxdFpsFARUREZBp6D0rXhimNRoM5c+agZs2acHR0xJUrVwAA06ZNw5o1a0xTJRkFAxUREZFpGHyV39y5c7Fu3TrMnz8fCoVCWt60aVOLux3Ls4ZjqIiIiEzD4EC1fv16rFq1CpGRkbCyspKWN2/eHAkJCUYtjoxL20N1jz1URERERmVwoEpKSkL9+vWLLNdoNFCr1UYpikzD47H7+Rl4xyEiIiIqhcGBKjAwEAcPHiyy/KeffkJQUJBRiiLT8HAsPEWrLhBIf8jwS0REZCwGT+w5ffp0DBkyBElJSdBoNPj555+RmJiI9evXY/v27aaokYxEaW0FFzsbpD9U425GLlztFWVvRERERGUyuIeqT58+2LZtG3bv3g0HBwdMnz4d8fHx2LZtG7p3726KGsmIeKUfERGR8RncQwUAHTt2RGxsrLFroUrg6ajEpZRMXulHRERkRHr1UHEAc9XBHioiIiLj0ytQNWnSBD/88APy8vJKbXfx4kW8/fbb+Pjjj41SHBmfpxPnoiIiIjI2vU75LV26FJMnT8Y777yD7t27o1WrVvDx8YGtrS0ePHiACxcu4NChQzh//jyio6Px9ttvm7puKicPR/ZQERERGZtegSo0NBS///47Dh06hI0bN2LDhg24fv06Hj58CA8PDwQFBeGNN95AZGQk3NzcTF0zVQBP+RERERmfQYPSO3TogA4dOpiqFqoEDFRERETGZ/C0CWTZtPfzu8cxVEREREbDQPWM0fZQ3c/KQ36BxszVEBERVQ0MVM8YdwcF5DJACCA1q/SrNomIiEg/DFTPGCu5DNUcOXUCERGRMTFQPYM4dQIREZFxGRyoOnfujPXr1+Phw4emqIcqAa/0IyIiMi6DA1VQUBDee+89eHt7IyoqCseOHTNFXWRCnjzlR0REZFQGB6rFixfj1q1bWLt2LVJSUtCpUycEBgZi4cKFuHPnjilqJCNjDxUREZFxlWsMlbW1Nfr164dff/0VN2/exGuvvYZp06bB19cXffv2xd69e41dJxkRAxUREZFxVWhQ+okTJzBjxgwsWrQIXl5emDp1Kjw8PPDiiy/ivffeM1aNZGQMVERERMZlcKBKSUnBokWL0LRpU3Ts2BF3797Ff/7zH1y7dg2zZs3C6tWrsWvXLqxYscKohaampiIyMhLOzs5wdXXFiBEjkJmZWeo2q1atQpcuXeDs7AyZTIa0tLQiberUqQOZTKbz+Pjjj3Xa/PHHH+jYsSNsbW3h6+uL+fPnG/OlVTrOlk5ERGRcBt3LDwBq1aqFevXqYfjw4Rg6dCg8PT2LtGnWrBlat25tlAK1IiMjcfv2bcTGxkKtVmPYsGEYOXIkvv/++xK3yc7ORkREBCIiIjB16tQS282ePRtRUVHScycnJ+nfKpUKPXr0QFhYGFasWIE///wTw4cPh6urK0aOHGmcF1fJ2ENFRERkXAYHqj179qBjx46ltnF2dsa+ffvKXdST4uPjERMTg5MnT6JVq1YAgKVLl6Jnz55YuHAhfHx8it1u3LhxAID9+/eXun8nJyd4e3sXu27Dhg3Iy8vD119/DYVCgSZNmiAuLg6ffvqp5QaqRz1Uqpx85KgLYGtjZeaKiIiILJvBp/zKClOmcPToUbi6ukphCgDCwsIgl8tx/PjxCu//448/RrVq1RAUFIQFCxYgPz9f59idOnWCQqGQloWHhyMxMREPHjwodn+5ublQqVQ6j6eJs501FFaF33qe9iMiIqo4g3uogoKCIJPJiiyXyWSwtbVF/fr1MXToUHTt2tUoBQJAcnIyvLy8dJZZW1vD3d0dycnJFdr32LFj0bJlS7i7u+PIkSOYOnUqbt++jU8//VQ6tr+/v8421atXl9a5ubkV2ee8efMwa9asCtVlSjKZDJ5OSiSlPcTdjFzUcrM3d0lEREQWzeAeqoiICFy5cgUODg7o2rUrunbtCkdHR1y+fBmtW7fG7du3ERYWhl9//bXMfU2ZMqXIgPAnHwkJCeV6YfqaMGECunTpgmbNmmHUqFFYtGgRli5ditzc8vfcTJ06Fenp6dLjxo0bRqzYODw4joqIiMhoDO6hunfvHiZOnIhp06bpLJ87dy6uX7+OXbt2YcaMGZgzZw769OlT6r4mTpyIoUOHltqmbt268Pb2RkpKis7y/Px8pKamljj2qbzatm2L/Px8XLt2DY0aNYK3t3eRCUu1z0s6tlKphFKpNGpdxsbZ0omIiIzH4ED1448/4tSpU0WWDxo0CMHBwfjqq6/w6quvSqfMSuPp6VnsVYJPCgkJQVpaGk6dOoXg4GAAwN69e6HRaNC2bVtDX0Kp4uLiIJfLpVOMISEh+OCDD6BWq2FjYwMAiI2NRaNGjYo93WcptFf63cvIM3MlREREls/gU362trY4cuRIkeVHjhyBra0tAECj0Uj/NobGjRsjIiICUVFROHHiBA4fPozo6GgMGjRIusIvKSkJAQEBOHHihLRdcnIy4uLicOnSJQDAn3/+ibi4OKSmpgIoHHC+ePFinD17FleuXMGGDRswfvx4DB48WApLr732GhQKBUaMGIHz589j48aNWLJkCSZMmGC012cO0tQJmTlmroSIiMjyGdxDNWbMGIwaNQqnTp2S5po6efIkVq9ejX//+98AgJ07d6JFixZGLXTDhg2Ijo5GaGgo5HI5+vfvj88//1xar1arkZiYiOzsbGnZihUrdAaHd+rUCQCwdu1aDB06FEqlEj/88ANmzpyJ3Nxc+Pv7Y/z48TphycXFBbt27cLo0aMRHBwMDw8PTJ8+3WKnTNDydCy8apFjqIiIiCpOJoQQhm60YcMGLFu2DImJiQCARo0aYcyYMXjttdcAAA8fPpSu+qPCyUFdXFyQnp4OZ2dnc5cDAIg5dxujvjuNlrVd8fM77c1dDhER0VPHkL/fBvVQ5efn46OPPsLw4cMRGRlZYjs7OztDdktm8M8pP/ZQERERVZRBY6isra0xf/58nYkvyTJ5Ohb2Ht7NyEU5OimJiIjoMQYPSg8NDcWBAwdMUQtVIg+nwjFUOWoNMnMZkImIiCrC4EHpL7zwAqZMmYI///wTwcHBcHBw0Fnfu3dvoxVHpmOvsIaj0hqZufm4l5kHJ1sbc5dERERksQwOVO+88w4AFDvPlEwmQ0FBQcWrokrh6aREZm4+7mbkwt/DoewNiIiIqFgGn/LTaDQlPhimLIs0WzqnTiAiIqoQgwPV43JyOCmkJdOOo7qbwe8jERFRRRgcqAoKCjBnzhzUrFkTjo6OuHLlCgBg2rRpWLNmjdELJNPh/fyIiIiMw+BA9eGHH2LdunWYP38+FAqFtLxp06ZYvXq1UYsj05LmouIpPyIiogoxOFCtX78eq1atQmRkJKysrKTlzZs3R0JCglGLI9NioCIiIjIOgwNVUlIS6tevX2S5RqOBWq02SlFUOThbOhERkXEYHKgCAwNx8ODBIst/+uknBAUFGaUoqhza2dLvZeSZuRIiIiLLZvA8VNOnT8eQIUOQlJQEjUaDn3/+GYmJiVi/fj22b99uihrJRLQ9VPcyc6HRCMjlMjNXREREZJkM7qHq06cPtm3bht27d8PBwQHTp09HfHw8tm3bhu7du5uiRjKRao6FFxXkawTSHvJ0LRERUXkZ3EMFAB07dkRsbKyxa6FKZmMlh5u9DR5kq3E3IxfuDoqyNyIiIqIiyhWoACAvLw8pKSnQaDQ6y2vXrl3hoqjyeDoppUDVyNvJ3OUQERFZJIMD1cWLFzF8+HAcOXJEZ7kQgvfys0CeTkr8dScTdzM5WzoREVF5GRyohg4dCmtra2zfvh01atSATMaBzJaM9/MjIiKqOIMDVVxcHE6dOoWAgABT1EOV7J8r/Th1AhERUXmVax6qe/fumaIWMgPOlk5ERFRxBgeqTz75BJMmTcL+/ftx//59qFQqnQdZFg+e8iMiIqowg0/5hYWFAQBCQ0N1lnNQumViDxUREVHFGRyo9u3bZ4o6yEx4Pz8iIqKKMzhQde7c2RR1kJlor/JLzcqDukADGyuDzwITERE988r11/PgwYMYPHgw2rVrh6SkJADAt99+i0OHDhm1ODI9N3sFrB7dw+8+r/QjIiIqF4MD1ebNmxEeHg47OzucPn0aubmFp4rS09Px0UcfGb1AMi25XAaPR/f0u8fTfkREROVicKCaO3cuVqxYga+++go2NjbS8vbt2+P06dNGLY4qBwemExERVYzBgSoxMRGdOnUqstzFxQVpaWnGqIkqGWdLJyIiqhiDA5W3tzcuXbpUZPmhQ4dQt25doxRFlUuai4qn/IiIiMrF4EAVFRWFd999F8ePH4dMJsOtW7ewYcMGvPfee3j77bdNUSOZGE/5ERERVYzBgWrKlCl47bXXEBoaiszMTHTq1Alvvvkm3nrrLYwZM8YUNQIAUlNTERkZCWdnZ7i6umLEiBHIzMwsdZtVq1ahS5cucHZ2hkwmK3JKcv/+/ZDJZMU+Tp48CQC4du1aseuPHTtmqpda6RioiIiIKsbgeahkMhk++OADvP/++7h06RIyMzMRGBgIR0dHU9QniYyMxO3btxEbGwu1Wo1hw4Zh5MiR+P7770vcJjs7GxEREYiIiMDUqVOLrG/Xrh1u376ts2zatGnYs2cPWrVqpbN89+7daNKkifS8WrVqFXxFTw8GKiIioooxOFBpKRQKBAYGGrOWEsXHxyMmJgYnT56Ugs7SpUvRs2dPLFy4ED4+PsVuN27cOACFPVHFUSgU8Pb2lp6r1Wr8+uuvGDNmDGQymU7batWq6bStSrSD0jltAhERUflYxLTYR48ehaurq06vUVhYGORyOY4fP26042zduhX379/HsGHDiqzr3bs3vLy80KFDB2zdutVox3wasIeKiIioYsrdQ1WZkpOT4eXlpbPM2toa7u7uSE5ONtpx1qxZg/DwcNSqVUta5ujoiEWLFqF9+/aQy+XYvHkz+vbti19++QW9e/cudj+5ubnShKcAoFKpjFajKWgDVUZuPh7mFcBOYWXmioiIiCyLWXuopkyZUuKgcO0jISGhUmq5efMmdu7ciREjRugs9/DwwIQJE9C2bVu0bt0aH3/8MQYPHowFCxaUuK958+bBxcVFevj6+pq6/ApxVFpDaV34o8DTfkRERIYzaw/VxIkTMXTo0FLb1K1bF97e3khJSdFZnp+fj9TUVKONa1q7di2qVatWYq/T49q2bYvY2NgS10+dOhUTJkyQnqtUqqc6VMlkMng6KXHzwUOkZOTC193e3CURERFZFL0ClSFjhvQJJFqenp7w9PQss11ISAjS0tJw6tQpBAcHAwD27t0LjUaDtm3b6n28kgghsHbtWrzxxhs6t9MpSVxcHGrUqFHieqVSCaVSWeG6KpM2UHEcFRERkeH0ClR9+/bVa2cymQwFBQUVqadYjRs3RkREBKKiorBixQqo1WpER0dj0KBB0hV+SUlJCA0Nxfr169GmTRsAhWOvkpOTpZnd//zzTzg5OaF27dpwd3eX9r93715cvXoVb775ZpFjf/PNN1AoFAgKCgIA/Pzzz/j666+xevVqo79Oc/LkbOlERETlpleg0mg0pq6jTBs2bEB0dDRCQ0Mhl8vRv39/fP7559J6tVqNxMREZGdnS8tWrFiBWbNmSc+19yBcu3atzqnGNWvWoF27dggICCj22HPmzMH169dhbW2NgIAAbNy4EQMGDDDyKzQv7cD0e+yhIiIiMphMCCHMXURVp1Kp4OLigvT0dDg7O5u7nGIt3v0XFu++iNfa1sZHLz9n7nKIiIjMzpC/3+UalJ6VlYUDBw7g77//Rl5ens66sWPHlmeXZGaci4qIiKj8DA5UZ86cQc+ePZGdnY2srCy4u7vj3r17sLe3h5eXFwOVhZLGUDFQERERGczgeajGjx+Pl156CQ8ePICdnR2OHTuG69evIzg4GAsXLjRFjVQJPNhDRUREVG4GB6q4uDhMnDgRcrkcVlZWyM3Nha+vL+bPn49///vfpqiRKsHjV/lxWB0REZFhDA5UNjY2kMsLN/Py8sLff/8NAHBxccGNGzeMWx1VGu0Yqrx8DVQ5+WauhoiIyLIYPIYqKCgIJ0+eRIMGDdC5c2dMnz4d9+7dw7fffoumTZuaokaqBLY2VnCytUZGTj7uZebCxa7sCU6JiIiokME9VB999JE0S/iHH34INzc3vP3227h79y5Wrlxp9AKp8vBKPyIiovIxuIeqVatW0r+9vLwQExNj1ILIfDwdlbhyN4uBioiIyEAG91B169YNaWlpRZarVCp069bNGDWRmbCHioiIqHwMDlT79+8vMpknAOTk5ODgwYNGKYrMw4P38yMiIioXvU/5/fHHH9K/L1y4gOTkZOl5QUEBYmJiULNmTeNWR5WKPVRERETlo3egatGiBWQyGWQyWbGn9uzs7LB06VKjFkeVi4GKiIiofPQOVFevXoUQAnXr1sWJEyfg6ekprVMoFPDy8oKVlZVJiqTKoQ1U93jKj4iIyCB6Byo/Pz8AgEajMVkxZF68nx8REVH5GDxtAgBcvnwZixcvRnx8PAAgMDAQ7777LurVq2fU4qhyeT3qobqflYcCjYCVXGbmioiIiCyDwVf57dy5E4GBgThx4gSaNWuGZs2a4fjx42jSpAliY2NNUSNVEncHBWQyoEAj8CC76JWcREREVDyDe6imTJmC8ePH4+OPPy6yfPLkyejevbvRiqPKZW0lh7u9Avez8nA3I1eaRoGIiIhKZ3APVXx8PEaMGFFk+fDhw3HhwgWjFEXmwyv9iIiIDGdwoPL09ERcXFyR5XFxcfDy8jJGTWRGDFRERESG0/uU3+zZs/Hee+8hKioKI0eOxJUrV9CuXTsAwOHDh/HJJ59gwoQJJiuUKof2Sj9OnUBERKQ/vQPVrFmzMGrUKEybNg1OTk5YtGgRpk6dCgDw8fHBzJkzMXbsWJMVSpWDPVRERESG0ztQCSEAADKZDOPHj8f48eORkZEBAHBycjJNdVTppEDFHioiIiK9GXSVn0ymOy8Rg1TVwx4qIiIiwxkUqBo2bFgkVD0pNTW1QgWReXG2dCIiIsMZFKhmzZoFFxcXU9VCTwEPnvIjIiIymEGBatCgQZwaoYrT9lClZauRl6+BwtrgmTWIiIieOXr/tSzrVB9VDS52NrCxKvxe389iLxUREZE+9A5U2qv8qGqTy2XSLWc4joqIiEg/ep/y02g0pqyDniKeTkrcTs9hoCIiItITB8hQEbzSj4iIyDAMVFQE56IiIiIyjMUEqtTUVERGRsLZ2Rmurq4YMWIEMjMzS20/ZswYNGrUCHZ2dqhduzbGjh2L9PR0nXZ///03evXqBXt7e3h5eeH9999Hfn6+Tpv9+/ejZcuWUCqVqF+/PtatW2eKl/jUkMZQceoEIiIivVhMoIqMjMT58+cRGxuL7du347fffsPIkSNLbH/r1i3cunULCxcuxLlz57Bu3TrExMRgxIgRUpuCggL06tULeXl5OHLkCL755husW7cO06dPl9pcvXoVvXr1QteuXREXF4dx48bhzTffxM6dO036es2JPVRERESGkQkLuHwvPj4egYGBOHnyJFq1agUAiImJQc+ePXHz5k34+PjotZ9NmzZh8ODByMrKgrW1Nf73v//hxRdfxK1bt1C9enUAwIoVKzB58mTcvXsXCoUCkydPxo4dO3Du3DlpP4MGDUJaWhpiYmL0Oq5KpYKLiwvS09Ph7Oxs4KuvfP/98zbe2XAareu4YdOoduYuh4iIyCwM+fttET1UR48ehaurqxSmACAsLAxyuRzHjx/Xez/aN8Ta2lra73PPPSeFKQAIDw+HSqXC+fPnpTZhYWE6+wkPD8fRo0dLPE5ubi5UKpXOw5Kwh4qIiMgwFhGokpOTi8zQbm1tDXd3dyQnJ+u1j3v37mHOnDk6pwmTk5N1whQA6bl2vyW1UalUePjwYbHHmjdvHlxcXKSHr6+vXjU+LXiVHxERkWHMGqimTJkCmUxW6iMhIaHCx1GpVOjVqxcCAwMxc+bMihdehqlTpyI9PV163Lhxw+THNCZtD1VWXgGycvPLaE1EREQG3cvP2CZOnIihQ4eW2qZu3brw9vZGSkqKzvL8/HykpqbC29u71O0zMjIQEREBJycnbNmyBTY2NtI6b29vnDhxQqf9nTt3pHXar9plj7dxdnaGnZ1dscdUKpVQKpWl1vU0c1Baw15hhey8AtzLzIWD0qw/JkRERE89s/6l9PT0hKenZ5ntQkJCkJaWhlOnTiE4OBgAsHfvXmg0GrRt27bE7VQqFcLDw6FUKrF161bY2toW2e+HH36IlJQU6ZRibGwsnJ2dERgYKLX573//q7NdbGwsQkJCDHqtlsbTSYnr97NxNyMXftUczF0OERHRU80ixlA1btwYERERiIqKwokTJ3D48GFER0dj0KBB0hV+SUlJCAgIkHqcVCoVevTogaysLKxZswYqlQrJyclITk5GQUEBAKBHjx4IDAzE66+/jrNnz2Lnzp34v//7P4wePVrqYRo1ahSuXLmCSZMmISEhAcuXL8ePP/6I8ePHm+fNqCS8nx8REZH+LOZczoYNGxAdHY3Q0FDI5XL0798fn3/+ubRerVYjMTER2dnZAIDTp09LVwDWr19fZ19Xr15FnTp1YGVlhe3bt+Ptt99GSEgIHBwcMGTIEMyePVtq6+/vjx07dmD8+PFYsmQJatWqhdWrVyM8PLwSXrX5aAem3+PknkRERGWyiHmoLJ2lzUMFANN+OYdvj13H2G71MaFHI3OXQ0REVOmq3DxUVPmkuajYQ0VERFQmBioqFif3JCIi0h8DFRWLk3sSERHpj4GKisUeKiIiIv0xUFGxPJy0V/nlgdctEBERlY6Biorl4agAAOQVaKB6yNvPEBERlYaBioqltLaCi13hbXruZuaYuRoiIqKnGwMVlUg7jiqF46iIiIhKxUBFJeKVfkRERPphoKIS8Uo/IiIi/TBQUYk4WzoREZF+GKioRB485UdERKQXBioqUXXnwkB1/X62mSshIiJ6ujFQUYla13EHAMTdSEN6ttrM1RARET29GKioRL7u9mhY3REFGoEDF++auxwiIqKnFgMVlaprgBcAYG/8HTNXQkRE9PRioKJShQZUBwDs/+suCjS8px8REVFxGKioVC1ru8LFzgZp2Wqc+fuBucshIiJ6KjFQUamsreTo3NATALA3IcXM1RARET2dGKioTN2046gYqIiIiIrFQEVl6tzQE3IZkJCcgaS0h+Yuh4iI6KnDQEVlcnNQINjPDQB7qYiIiIrDQEV64fQJREREJWOgIr1op084cvk+HuYVmLkaIiKipwsDFemlYXVH1HS1Q26+Bkcu3zN3OURERE8VBirSi0wm49V+REREJWCgIr09HqiE4KzpREREWgxUpLeQetVgayPH7fQcJCRnmLscIiKipwYDFenN1sYKHep7AOBpPyIioscxUJFBtNMn7OH0CURERBKLCVSpqamIjIyEs7MzXF1dMWLECGRmZpbafsyYMWjUqBHs7OxQu3ZtjB07Funp6VKbs2fP4tVXX4Wvry/s7OzQuHFjLFmyRGc/+/fvh0wmK/JITk422Wt9mmnHUZ25kYbUrDwzV0NERPR0sDZ3AfqKjIzE7du3ERsbC7VajWHDhmHkyJH4/vvvi21/69Yt3Lp1CwsXLkRgYCCuX7+OUaNG4datW/jpp58AAKdOnYKXlxe+++47+Pr64siRIxg5ciSsrKwQHR2ts7/ExEQ4OztLz728vEz3Yp9iNVzs0LiGM+Jvq7A/MQX9WtYyd0lERERmJxMWcLlWfHw8AgMDcfLkSbRq1QoAEBMTg549e+LmzZvw8fHRaz+bNm3C4MGDkZWVBWvr4rPk6NGjER8fj7179wIo7KHq2rUrHjx4AFdX13LVr1Kp4OLigvT0dJ1QZqkW7kzEsn2X8GKzGlj2Wktzl0NERGQShvz9tohTfkePHoWrq6sUpgAgLCwMcrkcx48f13s/2jekpDClbePu7l5keYsWLVCjRg10794dhw8fLvU4ubm5UKlUOo+qRDuO6sBfd6Eu0Ji5GiIiIvOziECVnJxc5BSbtbU13N3d9R7LdO/ePcyZMwcjR44ssc2RI0ewceNGnTY1atTAihUrsHnzZmzevBm+vr7o0qULTp8+XeJ+5s2bBxcXF+nh6+urV42WooWvK9wdFMjIycep6w/MXQ4REZHZmTVQTZkypdgB348/EhISKnwclUqFXr16ITAwEDNnziy2zblz59CnTx/MmDEDPXr0kJY3atQIb731FoKDg9GuXTt8/fXXaNeuHT777LMSjzd16lSkp6dLjxs3blT4NTxNrOQydGnkCYDTJxAREQFmHpQ+ceJEDB06tNQ2devWhbe3N1JSdP9w5+fnIzU1Fd7e3qVun5GRgYiICDg5OWHLli2wsbEp0ubChQsIDQ3FyJEj8X//939l1t2mTRscOnSoxPVKpRJKpbLM/ViybgFe+Pl0EvYmpODfPRubuxwiIiKzMmug8vT0hKenZ5ntQkJCkJaWhlOnTiE4OBgAsHfvXmg0GrRt27bE7VQqFcLDw6FUKrF161bY2toWaXP+/Hl069YNQ4YMwYcffqhX3XFxcahRo4Zebauqjg08YS2X4VJKJq7fz4JfNQdzl0RERGQ2FjGGqnHjxoiIiEBUVBROnDiBw4cPIzo6GoMGDZKu8EtKSkJAQABOnDgBoDBM9ejRA1lZWVizZg1UKhWSk5ORnJyMgoICAIWn+bp27YoePXpgwoQJ0vq7d+9Kx168eDF+/fVXXLp0CefOncO4ceOwd+9ejB49uvLfiKeIi50NWtVxA8DTfkRERBYzD9WGDRsQHR2N0NBQyOVy9O/fH59//rm0Xq1WIzExEdnZ2QCA06dPS1cA1q9fX2dfV69eRZ06dfDTTz/h7t27+O677/Ddd99J6/38/HDt2jUAQF5eHiZOnIikpCTY29ujWbNm2L17N7p27WriV/z0Cw2ojmNXUrE3IQXD2vubuxwiIiKzsYh5qCxdVZuHSutSSibCPj0AhZUcZ6Z3h4PSYvI5ERFRmarcPFT0dKrn6QC/avbIK9Dg0KV75i6HiIjIbBioqNxkMhm6NiqcH2xvPMdRERHRs4uBiioktHFhoNqXmAKNhmePiYjo2cRARRXSxt8dDgorpGTk4vytqnWLHSIiIn0xUFGFKK2t0KGBBwBOn0BERM8uBiqqsNCA6gCAvQl3zFwJERGReTBQUYV1CSic7f7szXSkZOSYuRoiIqLKx0BFFeblZItmtVwAAPsT75bRmoiIqOphoCKj4PQJRET0LGOgIqPQTp9w8OJd5OVrzFwNERFR5WKgIqNo6uMCTyclsvIKcOJqqrnLISIiqlQMVGQUcrkMXRsVDk7n9AlERPSsYaAio+nG6ROIiOgZxUBFRtOhgQdsrGS4dj8bV+5mmrscIiKiSsNARUbjqLTG83WrAeBpPyIierYwUJFRaadP2MPpE4iI6BnCQEVGpZ0+4eS1VKhy1GauhoiIqHIwUJFR+VVzQD1PB+RrBA7+dc/c5RAREVUKBioyum4Bj2ZN5zgqIiJ6RjBQkdFpp0/Yn5gCjUaYuRoiIiLTY6Aio2tVxw1Otta4n5WHszfTzF0OERGRyTFQkdHZWMnRqSFnTSciomcHAxWZRLdGHEdFRETPDgYqMokujTwhkwHnb6mQnJ5j7nKIiIhMioGKTKKaoxJBvq4A2EtFRERVHwMVmQynTyAiomcFAxWZjHb6hMOX7iFHXWDmaoiIiEyHgYpMpnENJ9RwscVDdQGOXblv7nKIiIhMhoGKTEYmk6ErT/sREdEzgIGKTOrx6ROE4KzpRERUNVlMoEpNTUVkZCScnZ3h6uqKESNGIDMzs9T2Y8aMQaNGjWBnZ4fatWtj7NixSE9P12knk8mKPH744QedNvv370fLli2hVCpRv359rFu3zhQvsUpqX98DSms5bj54iIspJX+/iIiILJnFBKrIyEicP38esbGx2L59O3777TeMHDmyxPa3bt3CrVu3sHDhQpw7dw7r1q1DTEwMRowYUaTt2rVrcfv2benRt29fad3Vq1fRq1cvdO3aFXFxcRg3bhzefPNN7Ny50xQvs8qxU1ghpF41ADztR0REVZdMWMB5mPj4eAQGBuLkyZNo1aoVACAmJgY9e/bEzZs34ePjo9d+Nm3ahMGDByMrKwvW1tYACnuotmzZohOiHjd58mTs2LED586dk5YNGjQIaWlpiImJ0eu4KpUKLi4uSE9Ph7Ozs17bVCXfHr2Gab+eR5s67vhxVIi5yyEiItKLIX+/LaKH6ujRo3B1dZXCFACEhYVBLpfj+PHjeu9H+4Zow5TW6NGj4eHhgTZt2uDrr7/WGetz9OhRhIWF6bQPDw/H0aNHy/lqnj3agem/X09FWnaemashIiIyPuuym5hfcnIyvLy8dJZZW1vD3d0dycnJeu3j3r17mDNnTpHThLNnz0a3bt1gb2+PXbt24Z133kFmZibGjh0rHbt69eo621SvXh0qlQoPHz6EnZ1dkWPl5uYiNzdXeq5SqfSqsaqq5WaPRtWdkHgnAwf+uos+LWqauyQiIiKjMmsP1ZQpU4odFP74IyEhocLHUalU6NWrFwIDAzFz5kydddOmTUP79u0RFBSEyZMnY9KkSViwYEGFjjdv3jy4uLhID19f3wrtryro1pjTJxARUdVl1kA1ceJExMfHl/qoW7cuvL29kZKi+4c4Pz8fqamp8Pb2LvUYGRkZiIiIgJOTE7Zs2QIbG5tS27dt2xY3b96Uepi8vb1x584dnTZ37tyBs7Nzsb1TADB16lSkp6dLjxs3bpT1VlR52tvQHPjrLvILNGauhoiIyLjMesrP09MTnp6eZbYLCQlBWloaTp06heDgYADA3r17odFo0LZt2xK3U6lUCA8Ph1KpxNatW2Fra1vmseLi4uDm5galUikd+7///a9Om9jYWISElDy4WqlUSttToSBfV7ja2yAtW40zN9LQuo67uUsiIiIyGosYlN64cWNEREQgKioKJ06cwOHDhxEdHY1BgwZJV/glJSUhICAAJ06cAFAYpnr06IGsrCysWbMGKpUKycnJSE5ORkFB4X3ltm3bhtWrV+PcuXO4dOkSvvzyS3z00UcYM2aMdOxRo0bhypUrmDRpEhISErB8+XL8+OOPGD9+fOW/ERbM2kqOzg0LwzNP+xERUVVjEYEKADZs2ICAgACEhoaiZ8+e6NChA1atWiWtV6vVSExMRHZ2NgDg9OnTOH78OP7880/Ur18fNWrUkB7aU3A2Njb44osvEBISghYtWmDlypX49NNPMWPGDGm//v7+2LFjB2JjY9G8eXMsWrQIq1evRnh4eOW+AVWA9rTf3ngGKiIiqlosYh4qS/esz0OllZadh+C5u1GgETg0uStqudmbuyQiIqISVbl5qKhqcLVXILi2GwBgH0/7ERFRFcJARZVKO33CHgYqIiKqQhioqFJpx1EduXwf2Xn5Zq6GiIjIOBioqFI18HJELTc75OVrcOTSfXOXQ0REZBQMVFSpZDLZP1f7JfK0HxERVQ0MVFTptIFqX0IKeJEpERFVBQxUVOmer1sNdjZWuJ2eg/jbGeYuh4iIqMLMeusZejbZ2lihfX0P7I6/g70JdxDo8+zOzUVEVNmEEMjIzced9BzcTs9BsioHd9JzkPZQDVc7G1RzVKKaowIejkp4OCpQzVEJB4UVZDKZuUt/qjFQkVmENvZ6FKhSEN2tgbnLISMTQiA3X4PcfA3y8jVQFxR+zXvy6xPrch9//ti/cwt0n6sLCk8VW8tlsLaSw8ZKBiu5DDZW8n+WyWWwspLBRi6HtVXhMmu5DNbadlayR88L/21jJX+0j8Jlhv7tkEH/DazkMjjZWsPZ1gaOttawkvMPFRmHRiNwLysXd9JzcTv9Ie6oHgtNj/59Jz0HWXkFBu1XaS3XCVjVHBTwcHr09bEAVs1RAXd7Baytnr0TYAxUZBZdGxWOozpzIw2XUjLham8DuazwT5JcJoNMDt3nssKvclnhwHbtV6ocefkaPMjOQ2pWCY/sPKRm5uFBdh7uZ+XhQVYe8jUcH6cvR6W1FLCcbK3hbGdT7HMnWxs4P/rqYqd9bgNbG7nen4f8Ag1y8jXIURcgV/tVrUFOvu7X3Ceea9vn5hdAaW0FB6U1HG2t4ai0gqPSBg5KKzgpCwOi9t+G1GUphBBQFwjkawqDfYFGIL9AA/Wjr0IAMllhwJbJHv370e+yIsu1/8Zjv9dQuKC45QIC9zPzdHqVbqdrg9JD3FHl4o4qR+/PnrOtNWq42KG6iy28nZVwtVcgPVuN+1m5uJeZh3uZubiXmYscdeF/dpLSHiIp7aFe+3azL+zp0gYwDwcFXOxsoLSxgtJa/uhhBaVN4b8V2uePLVdYyR+t/2ebpzmoMVCRWXi72KKJjzPO31Ih7NMD5d6PXPZP4NL+4ikSxOSyImGs8LluUNN9Lnu0XdEg9/i21lYyKKwKfxnYPPqqeOyrzRPPFdbFtZVBYWX1aLns0S+WwjZyE/4xyteIwpCU+SgQPQpC9x//+mh9Rm7F5gz75z2SSe+BwqrwNWp/mUrvy2PvU4nv3aNfqmqNBvkFAvmP/pjlawTUBcUvK9D884cw/7Gv2j+EBdptNQL5BYaFQUMvrsgrEMjMVSNHrQEAZObmIzM3H7fTcwzaj5a1tsfrUfCyksuRqw1A6gLkPPa1oBKDrlxWGBYdH4Uvh0f/drK1hoNCG8gKHw6PltvZWEEjUBhUHn2PCv9dzPNH36/Snmv/XfD4z4G2bYGQfobUj7X952fon9CkfmzfTzuZDPB0VKKGiy2qO9sWfnWxlZ57O9vC28UW9gr9IkB2Xj7uSwErD/czc3E/q/D5/cy8wgCWUfg1NSsPGgE8yFbjQbYal4x8MbdcBp0gprS2kn5nejopsW5YG+Me0AAMVGQ2w9r7Y9ov55CTX4DyXuynEYBG2vjp/0VnyeQywM1eAXcHBdwcFKj2+Fd7Bao5KqT17g4KONlaS2GoqvVSGEtevgYZOWqocvILvz4s/JqRkw/Vo+Wqh/88f3xdxqN1GqENx4V/wAyh7QGwfdRroM9XhbX8Ud35yHoUBDNyH/1buywvH0IUfj5VOflQ5eQD6SZ6E58Cchmk08wymQxCCAgAQgACovDr4/9GYQgvTzZTWMt1glFxocnTUWnUnhx7hTXs3a3h6172/VcLNAJpj3qq72Xk4l7WowCWmYf0h+rCU/j5BYWn+B/1fOXma8P/o1P/Um9o4Tr1Y//B0QjgoboAD9VFT1l6O9sa7TWXB2+OXAl4c2T9aH/BaIR49Mv4sa949FVT+FVnWUnbCKBAaH+Z/bPvx9v8czyBAk3p6zWaf5YJIVAgCv+3Wtq4H3W+QF5BgTTuJy//iXFCBbrbPj7GSJ2vKTUilvXJFWUETCuZDK6PBSFtQHIv7mFf2F0v51ifp4oQAtl5BToBKyMnH/kaAdtHp0pK+qq0lpvs+6nRCDxUFxSGrceCV6Y2dOUVLs98LIhp12fnFcBKXjgmTjsezlouk8a3Wcn/GQtnbaVt96jN4+Pinnxu9dg+Hxs39/iYusfH4f0zzk4uHUtnjN5j6yr6PgpRNGhpw1fh+n/+42j/DA4O12jEYwGsQDeE5f+zXC6ToVNDT6Me25C/3+yhoqeGTCaDlQywMmBwL9GzTCaTweHR6bIaLuau5h9y+T91Vef/IcskezSE4NEzc5byVJLLZbCVW8HWxgqAjbnLKdHTO7qLiIiIyEIwUBERERFVEAMVERERUQUxUBERERFVEAMVERERUQUxUBERERFVEAMVERERUQUxUBERERFVEAMVERERUQUxUBERERFVEAMVERERUQUxUBERERFVEAMVERERUQUxUBERERFVkLW5C3gWCCEAACqVysyVEBERkb60f7e1f8dLw0BVCTIyMgAAvr6+Zq6EiIiIDJWRkQEXF5dS28iEPrGLKkSj0eDWrVtwcnKCTCYz6r5VKhV8fX1x48YNODs7G3XfpmbJtQOWXb8l1w5Ydv2WXDtg2fVbcu2AZddvqbULIZCRkQEfHx/I5aWPkmIPVSWQy+WoVauWSY/h7OxsUT+kj7Pk2gHLrt+Sawcsu35Lrh2w7PotuXbAsuu3xNrL6pnS4qB0IiIiogpioCIiIiKqIAYqC6dUKjFjxgwolUpzl2IwS64dsOz6Lbl2wLLrt+TaAcuu35JrByy7fkuuXV8clE5ERERUQeyhIiIiIqogBioiIiKiCmKgIiIiIqogBioiIiKiCmKgsmBffPEF6tSpA1tbW7Rt2xYnTpwwd0l6mTdvHlq3bg0nJyd4eXmhb9++SExMNHdZ5fLxxx9DJpNh3Lhx5i5Fb0lJSRg8eDCqVasGOzs7PPfcc/j999/NXVaZCgoKMG3aNPj7+8POzg716tXDnDlz9LrHljn89ttveOmll+Dj4wOZTIZffvlFZ70QAtOnT0eNGjVgZ2eHsLAwXLx40TzFPqG02tVqNSZPnoznnnsODg4O8PHxwRtvvIFbt26Zr+AnlPXeP27UqFGQyWRYvHhxpdVXGn1qj4+PR+/eveHi4gIHBwe0bt0af//9d+UXW4yy6s/MzER0dDRq1aoFOzs7BAYGYsWKFeYp1sgYqCzUxo0bMWHCBMyYMQOnT59G8+bNER4ejpSUFHOXVqYDBw5g9OjROHbsGGJjY6FWq9GjRw9kZWWZuzSDnDx5EitXrkSzZs3MXYreHjx4gPbt28PGxgb/+9//cOHCBSxatAhubm7mLq1Mn3zyCb788kssW7YM8fHx+OSTTzB//nwsXbrU3KUVKysrC82bN8cXX3xR7Pr58+fj888/x4oVK3D8+HE4ODggPDwcOTk5lVxpUaXVnp2djdOnT2PatGk4ffo0fv75ZyQmJqJ3795mqLR4Zb33Wlu2bMGxY8fg4+NTSZWVrazaL1++jA4dOiAgIAD79+/HH3/8gWnTpsHW1raSKy1eWfVPmDABMTEx+O677xAfH49x48YhOjoaW7dureRKTUCQRWrTpo0YPXq09LygoED4+PiIefPmmbGq8klJSREAxIEDB8xdit4yMjJEgwYNRGxsrOjcubN49913zV2SXiZPniw6dOhg7jLKpVevXmL48OE6y/r16yciIyPNVJH+AIgtW7ZIzzUajfD29hYLFiyQlqWlpQmlUin+85//mKHCkj1Ze3FOnDghAIjr169XTlEGKKn+mzdvipo1a4pz584JPz8/8dlnn1V6bWUprvaBAweKwYMHm6cgAxVXf5MmTcTs2bN1lrVs2VJ88MEHlViZabCHygLl5eXh1KlTCAsLk5bJ5XKEhYXh6NGjZqysfNLT0wEA7u7uZq5Ef6NHj0avXr10vgeWYOvWrWjVqhX+9a9/wcvLC0FBQfjqq6/MXZZe2rVrhz179uCvv/4CAJw9exaHDh3CCy+8YObKDHf16lUkJyfr/Py4uLigbdu2FvsZlslkcHV1NXcpetFoNHj99dfx/vvvo0mTJuYuR28ajQY7duxAw4YNER4eDi8vL7Rt27bUU5pPm3bt2mHr1q1ISkqCEAL79u3DX3/9hR49epi7tApjoLJA9+7dQ0FBAapXr66zvHr16khOTjZTVeWj0Wgwbtw4tG/fHk2bNjV3OXr54YcfcPr0acybN8/cpRjsypUr+PLLL9GgQQPs3LkTb7/9NsaOHYtvvvnG3KWVacqUKRg0aBACAgJgY2ODoKAgjBs3DpGRkeYuzWDaz2lV+Azn5ORg8uTJePXVVy3mpreffPIJrK2tMXbsWHOXYpCUlBRkZmbi448/RkREBHbt2oWXX34Z/fr1w4EDB8xdnl6WLl2KwMBA1KpVCwqFAhEREfjiiy/QqVMnc5dWYdbmLoCebaNHj8a5c+dw6NAhc5eilxs3buDdd99FbGzsUzNmwRAajQatWrXCRx99BAAICgrCuXPnsGLFCgwZMsTM1ZXuxx9/xIYNG/D999+jSZMmiIuLw7hx4+Dj4/PU115VqdVqvPLKKxBC4MsvvzR3OXo5deoUlixZgtOnT0Mmk5m7HINoNBoAQJ8+fTB+/HgAQIsWLXDkyBGsWLECnTt3Nmd5elm6dCmOHTuGrVu3ws/PD7/99htGjx4NHx8fi+vxfxJ7qCyQh4cHrKyscOfOHZ3ld+7cgbe3t5mqMlx0dDS2b9+Offv2oVatWuYuRy+nTp1CSkoKWrZsCWtra1hbW+PAgQP4/PPPYW1tjYKCAnOXWKoaNWogMDBQZ1njxo2fmiuESvP+++9LvVTPPfccXn/9dYwfP94iewq1n1NL/gxrw9T169cRGxtrMb1TBw8eREpKCmrXri19hq9fv46JEyeiTp065i6vVB4eHrC2trbYz/DDhw/x73//G59++ileeuklNGvWDNHR0Rg4cCAWLlxo7vIqjIHKAikUCgQHB2PPnj3SMo1Ggz179iAkJMSMlelHCIHo6Ghs2bIFe/fuhb+/v7lL0ltoaCj+/PNPxMXFSY9WrVohMjIScXFxsLKyMneJpWrfvn2RKSr++usv+Pn5maki/WVnZ0Mu1/2VZWVlJf2v3ZL4+/vD29tb5zOsUqlw/Phxi/gMa8PUxYsXsXv3blSrVs3cJent9ddfxx9//KHzGfbx8cH777+PnTt3mru8UikUCrRu3dpiP8NqtRpqtbrKfI6fxFN+FmrChAkYMmQIWrVqhTZt2mDx4sXIysrCsGHDzF1amUaPHo3vv/8ev/76K5ycnKQxIy4uLrCzszNzdaVzcnIqMtbLwcEB1apVs4gxYOPHj0e7du3w0Ucf4ZVXXsGJEyewatUqrFq1ytyllemll17Chx9+iNq1a6NJkyY4c+YMPv30UwwfPtzcpRUrMzMTly5dkp5fvXoVcXFxcHd3R+3atTFu3DjMnTsXDRo0gL+/P6ZNmwYfHx/07dvXfEU/UlrtNWrUwIABA3D69Gls374dBQUF0mfY3d0dCoXCXGVLynrvnwyANjY28Pb2RqNGjSq71CLKqv3999/HwIED0alTJ3Tt2hUxMTHYtm0b9u/fb76iH1NW/Z07d8b7778POzs7+Pn54cCBA1i/fj0+/fRTM1ZtJGa+ypAqYOnSpaJ27dpCoVCINm3aiGPHjpm7JL0AKPaxdu1ac5dWLpY0bYIQQmzbtk00bdpUKJVKERAQIFatWmXukvSiUqnEu+++K2rXri1sbW1F3bp1xQcffCByc3PNXVqx9u3bV+zP+ZAhQ4QQhVMnTJs2TVSvXl0olUoRGhoqEhMTzVv0I6XVfvXq1RI/w/v27TN36UKIst/7Jz1N0yboU/uaNWtE/fr1ha2trWjevLn45ZdfzFfwE8qq//bt22Lo0KHCx8dH2NraikaNGolFixYJjUZj3sKNQCbEUzrNMBEREZGF4BgqIiIiogpioCIiIiKqIAYqIiIiogpioCIiIiKqIAYqIiIiogpioCIiIiKqIAYqIiIiogpioCIisgAymQy//PKLucsgohIwUBHRU+fu3bt4++23Ubt2bSiVSnh7eyM8PByHDx/WaXfmzBn861//QvXq1WFra4sGDRogKioKf/31FwDg2rVrkMlkRR6DBw8u8dhdunTBuHHjpOd16tTB4sWLTfEyizVz5ky0aNGiyPLbt2/jhRdeqLQ6iMgwvJcfET11+vfvj7y8PHzzzTeoW7cu7ty5gz179uD+/ftSm+3bt6N///4IDw/Hhg0bUK9ePaSkpGDTpk2YNm0aNm7cKLXdvXs3mjRpIj03xz0j8/LyKnSfO29vbyNWQ0RGZ+573xARPe7BgwcCgNi/f3+JbbKysoSHh4fo27dvifsQQkj3nTtz5ozex3/83oydO3cuck8yrYMHD4oOHToIW1tbUatWLTFmzBiRmZkprffz8xOzZ88Wr7/+unBycpLuZTZp0iTRoEEDYWdnJ/z9/cX//d//iby8PCGEEGvXri3xHpcAxJYtW6T9//HHH6Jr167C1tZWuLu7i6ioKJGRkSGtHzJkiOjTp49YsGCB8Pb2Fu7u7uKdd96RjkVExsVTfkT0VHF0dISjoyN++eUX5ObmFttm586duHfvHiZNmlTseldXV6PU8vPPP6NWrVqYPXs2bt++jdu3bwMALl++jIiICPTv3x9//PEHNm7ciEOHDiE6Olpn+4ULF6J58+Y4c+YMpk2bBgBwcnLCunXrcOHCBSxZsgRfffUVPvvsMwDAwIEDMXHiRDRp0kQ63sCBA4vUlZWVhfDwcLi5ueHkyZPYtGkTdu/eXeT4+/btw+XLl7Fv3z588803WLduHdatW2eU94aInmDuREdE9KSffvpJuLm5CVtbW9GuXTsxdepUcfbsWWn9J598IgCI1NTUUvej7aGys7MTDg4O0uP06dMlbvN4D5UQhT1Nn332mU6bESNGiJEjR+osO3jwoJDL5eLhw4fSdiX1oD1uwYIFIjg4WHo+Y8YM0bx58yLt8FgP1apVq4Sbm5tOj9iOHTuEXC4XycnJQojCHio/Pz+Rn58vtfnXv/4lBg4cWGZNRGQ49lAR0VOnf//+uHXrFrZu3YqIiAjs378fLVu2lHpXhBAG7W/jxo2Ii4uTHoGBgRWq7+zZs1i3bp3Um+bo6Ijw8HBoNBpcvXpVateqVatia2nfvj28vb3h6OiI//u//8Pff/9t0PHj4+PRvHlzODg4SMvat28PjUaDxMREaVmTJk1gZWUlPa9RowZSUlIMOhYR6YeBioieSra2tujevTumTZuGI0eOYOjQoZgxYwYAoGHDhgCAhIQEvfbl6+uL+vXrSw+lUlmh2jIzM/HWW2/phLSzZ8/i4sWLqFevntTu8cADAEePHkVkZCR69uyJ7du348yZM/jggw+Ql5dXoXpKYmNjo/NcJpNBo9GY5FhEzzpe5UdEFiEwMFCah6lHjx7w8PDA/PnzsWXLliJt09LSjDaOSqFQoKCgQGdZy5YtceHCBdSvX9+gfR05cgR+fn744IMPpGXXr18v83hPaty4MdatW4esrCwptB0+fBhyuRyNGjUyqCYiMg72UBHRU+X+/fvo1q0bvvvuO/zxxx+4evUqNm3ahPnz56NPnz4ACnt+Vq9ejR07dqB3797YvXs3rl27ht9//x2TJk3CqFGjjFZPnTp18NtvvyEpKQn37t0DAEyePBlHjhxBdHQ04uLicPHiRfz6669FBoU/qUGDBvj777/xww8/4PLly/j888+LBMI6derg6tWriIuLw71794odmB8ZGQlbW1sMGTIE586dw759+zBmzBi8/vrrqF69utFeOxHpj4GKiJ4qjo6OaNu2LT777DN06tQJTZs2xbRp0xAVFYVly5ZJ7fr06YMjR47AxsYGr732GgICAvDqq68iPT0dc+fONVo9s2fPxrVr11CvXj14enoCAJo1a4YDBw7gr7/+QseOHREUFITp06fDx8en1H317t0b48ePR3R0NFq0aIEjR45IV/9p9e/fHxEREejatSs8PT3xn//8p8h+7O3tsXPnTqSmpqJ169YYMGAAQkNDdd4fIqpcMmHo6E4iIiIi0sEeKiIiIqIKYqAiIiIiqiAGKiIiIqIKYqAiIiIiqiAGKiIiIqIKYqAiIiIiqiAGKiIiIqIKYqAiIiIiqiAGKiIiIqIKYqAiIiIiqiAGKiIiIqIKYqAiIiIiqqD/B2W4cvqW0KhKAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(energies[:, 0])\n",
    "plt.xticks(np.arange(0, energies.shape[0], 2))\n",
    "plt.xlabel(\"SCF Iteration\")\n",
    "plt.ylabel(\"Total energy (eV)\")\n",
    "plt.title(\"Self consistent field energy convergence\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8060c09f",
   "metadata": {},
   "source": [
    "Use PySCF to evaluate the Gaussian atomic orbitals on a uniform grid.\n",
    "Multiplying by the molecular orbital coefficients calculated previously with `nanoDFT`\n",
    "gives use the molecular orbitals as a function of space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "fc79c21d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(125000, 36)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L = 15.0\n",
    "n = 50\n",
    "axes = [np.linspace(-L, L, n) for _ in range(3)]\n",
    "grid = np.stack(np.meshgrid(*axes, indexing=\"ij\"), axis=-1)\n",
    "grid = grid.reshape(-1, 3)\n",
    "atomic_orbitals = mol.eval_gto('GTOval_sph', grid)\n",
    "\n",
    "C = _dbginfo[3]\n",
    "molecular_orbitals = atomic_orbitals @ C\n",
    "molecular_orbitals.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85db5266",
   "metadata": {},
   "source": [
    "Below we define functions for visualising the molecular orbitals of benzene using the [py3Dmol](https://3dmol.org/) package"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "63fbdbd2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def cube_data(axes, value) -> str:\n",
    "    fmt = \"cube format\\n\\n\"\n",
    "    x, y, z = axes\n",
    "    nx, ny, nz = [ax.shape[0] for ax in axes]\n",
    "    fmt += \"0 \" + \" \".join([f\"{v:12.6f}\" for v in [x[0], y[0], z[0]]]) + \"\\n\"\n",
    "    fmt += f\"{nx} \" + \" \".join([f\"{v:12.6f}\" for v in [x[1] - x[0], 0.0, 0.0]]) + \"\\n\"\n",
    "    fmt += f\"{ny} \" + \" \".join([f\"{v:12.6f}\" for v in [0.0, y[1] - y[0], 0.0]]) + \"\\n\"\n",
    "    fmt += f\"{nz} \" + \" \".join([f\"{v:12.6f}\" for v in [0.0, 0.0, z[1] - z[0]]]) + \"\\n\"\n",
    "\n",
    "    line = \"\"\n",
    "    for i in range(len(value)):\n",
    "        line += f\"{value[i]:12.6f}\"\n",
    "\n",
    "        if i % 6 == 0:\n",
    "            fmt += line + \"\\n\"\n",
    "            line = \"\"\n",
    "\n",
    "    return fmt\n",
    "\n",
    "\n",
    "def build_transferfn(value) -> dict:\n",
    "    v = np.percentile(value, [99.9, 75])\n",
    "    a = [0.02, 0.0005]\n",
    "    return {\n",
    "        \"transferfn\": [\n",
    "            {\"color\": \"blue\", \"opacity\": a[0], \"value\": -v[0]},\n",
    "            {\"color\": \"blue\", \"opacity\": a[1], \"value\": -v[1]},\n",
    "            {\"color\": \"white\", \"opacity\": 0.0, \"value\": 0.0},\n",
    "            {\"color\": \"red\", \"opacity\": a[1], \"value\": v[1]},\n",
    "            {\"color\": \"red\", \"opacity\": a[0], \"value\": v[0]},\n",
    "        ]\n",
    "    }\n",
    "\n",
    "\n",
    "def plot_orbital(orbital, mol):\n",
    "    xyzfmt = f\"{len(mol.atom)}\\n\\n\" + mol.tostring()\n",
    "    v = py3Dmol.view(data=xyzfmt, style={\"stick\": {\"radius\": 0.06}, \"sphere\": {\"radius\": 0.2}})\n",
    "    v.addVolumetricData(cube_data(axes, orbital), \"cube\", build_transferfn(orbital))\n",
    "    return v"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8efd393a",
   "metadata": {},
   "source": [
    "Try changing the `mo_index` variable to select the different molecular orbitals benzene."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "127815f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "mo_index = 5\n",
    "\n",
    "orbital = molecular_orbitals[:, mo_index]\n",
    "mol_view = plot_orbital(orbital, mol)\n",
    "mol_view.spin()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}