checks.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2310a17c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3cda1dd0",
   "metadata": {},
   "source": [
    "# Checks for the different functions of the code"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2cb53ce1",
   "metadata": {},
   "source": [
    "## Harmonic Oscillator eigenfunctions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "14b61966",
   "metadata": {},
   "outputs": [],
   "source": [
    "from basis_set import *\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "51ca5c14",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Norm of 0th eigenfunction: 0.99998\n",
      "Norm of 1th eigenfunction: 0.99998\n",
      "Norm of 2th eigenfunction: 0.99998\n",
      "Norm of 3th eigenfunction: 0.99998\n",
      "Norm of 4th eigenfunction: 0.99998\n",
      "Norm of 5th eigenfunction: 0.99998\n",
      "Norm of 6th eigenfunction: 0.99998\n",
      "Norm of 7th eigenfunction: 0.99998\n",
      "Norm of 8th eigenfunction: 0.99998\n",
      "Norm of 9th eigenfunction: 0.99998\n"
     ]
    }
   ],
   "source": [
    "# Normalization check (numerical integration from -L to L)\n",
    "L = 100\n",
    "N = 40000\n",
    "\n",
    "x = np.linspace(-L, L, N)\n",
    "Ax = 2*L/N\n",
    "\n",
    "for n in range(10):\n",
    "    norm = np.sum(HO_wf(x,n)**2)*Ax\n",
    "    print(\"Norm of {}th eigenfunction: {:0.5f}\".format(n, norm))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "6b7b89bb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7kAAAEWCAYAAACjclDSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABsaElEQVR4nO3dd3hsV3ku8PebXjTSqLfTi3vn2IAx4IAJ4BswoRNCSYGQhJIeuCGBS8gNJJB2QyBOqKGaEMCAKaYaAzYc93Ls03yayjnq0vS27h9779FImrL39Nnz/p5Hj3WmLsmf9qzyrW+JUgpEREREREREduBodQOIiIiIiIiI6oWDXCIiIiIiIrINDnKJiIiIiIjINjjIJSIiIiIiItvgIJeIiIiIiIhsg4NcIiIiIiIisg0OcptARD4hIu+t02t9U0ReV4/XqvA+oyJyh4isicgHG/1+m947IiJ7mvmeVBsRebeIfLpOr/UREfnLerxWhffxi8jXRGRFRL7Y6Pfb9N6PiMj1zXxPqg1j3PJ7M8Y7COPb8nszvjsM++KW37vj++Ic5FZBRAZE5MsiEhWRkyLyawX3vV5E7mzUeyulnq+U+mSjXr/AGwHMA+hVSv1xo95ERH4oIr9deJtSqkcpdbxR70mVicibReSgiCRF5BOb7rteRM406r2VUm9SSv11o16/wEsBjAIYVEq9rFFvUuyDVSl1sVLqh416TypPRLwi8lH9+r0mIveLyPML7meMW8AYbz8i8mkRmRGRVRE5XPg5y/i2hvHd3kRkv4gkCidp2Be3xq59cVerG9ChPgQgBe3iegWAb4jIA0qpR1raqvraCeBRpZRqdUOoJaYBvBfAcwH4W9yWRtkJ4LBSKtPqhlDTuQCcBvBMAKcA3AjgFhG5VCl1opUNqzPGePf6WwC/pZRKisgFAH4oIvcppe5pdcPqiPFNgNYn/0WrG9Eg7IvXQinFLwtfAILQBrjnFdz2XwDeB+BCAAkAWQARAMv6/Z+A9kf4DQBrAO4GsLfMezwFwE8BLAN4AMD1Bff9EMBv6987AXwQ2izPEwDeDEABcOn39wH4KIAZAFPQBi1O/b7XA7gTwAcALOnPf35Be9P6zxkBcIN+23sL2nE9gDMF/z4B4E8APAhgBcAXAPgK7r8JwP0AVgEcA/A8AH+j/64S+vv8q/5YBWBfwc/wKQBzAE4CeCcAR6WfgV91i/f3AvjEpviPA8jp/88iACYAvBvALfr/qzUAjwA4UOZ1LwBwO4BFAI8DeHnBfZtj7c/0GJ4G8Nub4sOr//8/BeAsgI8A8BfGKIA/BnBOf43f0O/7P3p8p/Wf4bf0n+HTBe+7Cxv/nn4I4K8B/ET/Gb8DYKjg8ddh/e/2tB6fb8TGv6WvFfy93FDwM/yT/vNN6997K/0M/Kp7rD8I4CWMcca43b4AnK//Xl/O+GZ82+kLwCv1uM3/vwf74uyLGz9vqxvQaV8ArgQQ23Tbn2D9wvd6AHduuv8TABYAXANtBeEzAD5f4vUn9cfeCC2d/Dn6v4f1+wv/sN4E4FEA2wD0A/jupj+sLwP4d2gfaiMAfg7gdwramQbwBv0P9HehXZyloM3v3fQzVPrD+jm0D8sBAIcAvEm/7xr9j+05+s80CeCCzT9PwWsV/mF9CsBXAYSgfWAdhjY7XfFn4Fdd4n3DILfY/3v9tndDu0DeqP+/+FsAd5V4zSC0DsRv6H8PV0L7cLhoc6xBuwDPArgYQADApzfFxz8CuFWPuRCArwH424J2ZgC8B4Bbb1sMQH9Bmz+96Weo1EE6BuA8aKvbPwTwPv2+ndA+NF+lv9cggCuK/e0U/L0YHaT3ALgL2t/oMLQP1b828zPwq25xPqrHr3FdYowzxjv+C8C/6b9LBeBeAD2Mb8a3Xb4A9ELrE24r8v/+9WBfvOv74tyTa10PtBmQQivQ/seX82Wl1M+VllbzGWhpzsX8OoDblFK3KaVySqnbARyE9oe22csB/LNS6oxSagnaajIAbbO6/pw/UEpFlVLnoH2YvLLg+SeVUv+hlMoC+CSAcWidvWr9i1JqWim1CO2DyvgZfwvAx5RSt+s/05RS6rFKLyYiTr2971BKrSktjfCDAF7TwJ+BqnenHrdZaNkNl5d43K8AOKGU+rhSKqOUug/AlwAU21P1cgAfV0o9opSKQfsgAwCIiECbZf9DpdSiUmoNwP/FxhhPA3iPUiqtlLoN2izl+TX8jB9XSh1WSsWhzR5fod/+awC+q5T6nP5eC0qp+02+5qv1Np5TSs1BW6EojPF6/wxUQETc0K7JnzRxXWKMM8Y7hlLq96D1TZ4O4H8AJCs8hfHN+O4kfw3go0opK/vL2Rfvor449+RaF4E2e1SoF9oMYDmzBd/HoA2Wi9kJ4GUi8oKC29wAflDksRPQZlMNhd/v1J83o32OANBmbgofk2+TUiqmP65Uu8zY/DNO6N9vB3BbFa83BO1nOFlw20los09b3rNOPwNVb/P/f5+IuNTW/VI7ATxZRJYLbnNB61RtNgHtg8VQGL/D0FYG7imIcYE2k2hY2PT+5f72zCj1d7wd2gpBNSawNcYnCv5d75+BdCLigBZ3KWgpZpUwxqvDGG8RvdN5p4j8OrQVln8p83DGd3UY300mIldAS9+90uJT2Re3rmP74hzkWncYgEtE9iuljui3XQ5t/wqgLe/X4jSA/1JKvcHEY2egpUcYtm96nSS0/Sb1KMoQhfZhZBiz8NzTAPaWuK/c72se2gzoTmipIACwA9qeBmqdesT4j5RSzzHx2HIxPg9tb9nFSql6xEStMX5Nifsq/b6mocW4cQ3Zod9GDaSvIn0U2mzzjUqpdMHdjPGtGOOdzYX1z2HG91aM785yPbS02VMFAyqniFyklLoK7IsX03V9caYrW6SUikJL+3mPiARF5GnQNnIbM5hnAWwTEU+Vb/FpAC8QkeeKiFNEfHq5/21FHnsLgLeJyKSIhAH8eUE7Z6AVVfigiPSKiENE9orIM6ts1/0AbtSPTxoD8AcWnvtRAL8hIs/W2zGpV3sEtN9X0XO49BnoWwD8jYiERGQngD+C9juiBhIRl4j4oM2oG3FoTIqdBTAoIn1VvvzXAZwnIq8REbf+dbWIXFjksbdAi50LRSQAIH/2olIqB+A/APyjiIzo7Z4UkedW2a77ATxDRHboP9s7LDz3MwBuEJGX67+7QX2mGSgT47rPAXiniAyLyBCAvwJjvBk+DK1AyQv01MVCjPGtGOMdQkRGROSVItKj9yOeC22v6ff0hzC+t2J8d5aboQ3YrtC/PgKtoJQRO+yLb9V1fXEOcqvze9CKFpyDdnH7XbV+fND3oc3mzYrIvNUXVkqdhjZo/t/QqpidBvCnKP7/6j+g/fE8COA+aGkIGWhV0gDgtQA80GZelgD8N7Q8+Wr8F7Tqcif09/yC2ScqpX4OrUDFP0Lbv/wjaDNCAPDPAF4qIksiUiyN6i3QZq6OQ6ve9lkAH6vuRyAL3glthv3t0PamxPXboO/h+ByA4yKyLCITJV+lCH3f1S9D2+MxDS3N5f3QKlRufuw3oaXX/QDAUWjFPYD1vWV/btwuIqvQCj5UtddJ33PzBWh/T/dA68iZfa5xDM0fQ6s2ej/W97N9FMBF+u/qK0We/l5o6XwPAngIWoGYuhxYT8XpH9K/A61zNCvaofcREXk1wBgv8VzGeOdQ0FKTz0D77P8AtD2BtwKM7xLPZXx3EKVUTCk1a3xB20qY0PdEA+yLb9GNfXGjehfZgIg8H8BHlFI7Kz6YqAPpKwUPQzuegWcjku0wxsnOGN9kd+yLtw+u5HYwEfGLyI16as0kgHdBK1VOZBsi8qsi4hWRfmirBV9j54jshDFOdsb4JjtjX7x9cZDb2QRaqfolaCkSh6DtBSGyk9+BtjXgGLT0n99tbXOI6o4xTnbG+CY7Y1+8TTFdmYiIiIiIiGyDK7lERERERERkG7Y7J3doaEjt2rWr1c0gm7vnnnvmlVLDrXhvxjg1A2Oc7I4xTnbXqhhnfFMzVIpv2w1yd+3ahYMHD7a6GWRzInKyVe/NGKdmYIyT3THGye5aFeOMb2qGSvHNdGUiIiIiIiKyDQ5yiYiIiIiIyDY4yCUiIiIiIiLb4CCXiIiIiIiIbIODXCIiIiIiIrINDnIJ8VS21U0gIiIiIiKqCw5yu5hSCn/4hftx4V99C//yvSOtbg4REREREVHNOMjtYt89dA5fvm8KwyEv/um7h3FsLtLqJhEREREREdWEg9wu9tm7T2K8z4evv+U6iAhuOXi61U0iIiIiIiKqSUsHuSLyMRE5JyIPl7hfRORfROSoiDwoIlc1u412FUlm8OMj83jh5RMY7fXh6fuH8K2HZ1vdLFthfJPdMcbJ7hjjZHeMcbKrVq/kfgLA88rc/3wA+/WvNwL4cBPa1BV+8cQiMjmFZ5w3DAB4+v5hnFyIYWo53uKW2conwPgme/sEGONkb58AY5zs7RNgjJMNtXSQq5S6A8BimYfcBOBTSnMXgLCIjDendfb202Pz8DgdeNLOfgDA0/YNarcfnW9ls2yF8U12xxgnu2OMk90xxsmuWr2SW8kkgMKNomf026hGD5xewSWTvfC5nQCA80dDCPlceODMcmsb1l0Y32R3jHGyO8Y42R1jnDpSuw9yTRGRN4rIQRE5ODc31+rmtL1cTuHRmVVcPNGXv01EcPFELx6eWm1hy6gUxjjZHWOc7I4xTnbG+KZ20+6D3CkA2wv+vU2/bQOl1M1KqQNKqQPDw8NNa1ynOr0UQySZwcUTvRtuv3iiD4/NriKbUy1qWdcxFd8AY5w6FmOc7I4xTnbHvjh1pHYf5N4K4LV6ZbenAFhRSs20ulGd7tFpbbX2oi2D3F4k0jkc53m5zcL4JrtjjJPdMcbJ7hjj1JFcrXxzEfkcgOsBDInIGQDvAuAGAKXURwDcBuBGAEcBxAD8Rmtaai+PTK/C6RCcNxracLsx6H10ZhX7N91H1jG+ye4Y42R3jHGyO8Y42VVLB7lKqVdVuF8B+P0mNadrHD67ht1DwXzRKcOuwSAcAhyfi7aoZfbC+Ca7Y4yT3THGye4Y42RX7Z6uTA3wxHwUu4eCW273uZ3Y1h/AMaYrExERERFRh+Igt8vkcgonF2NFB7kAsGc4yJVcIiIiIiLqWBzkdpnplThSmVzpQe5QD56YjyLHCsvUBf75u0fwls/dh7VEutVNISIiIqI6aemeXGq+J+a1Vdpdg8UHuXtHgoins5hZTWAy7G9m04ia6r5TS/jH7x4GAOwZCuIPn3Nei1tERERERPXAldwuc0If5O4ZLr2SC4DHCJHtfeuRWbgcgqt2hPHl+6ag1dYgIiIiok7HQW6XeWI+hoDHiZGQt+j9OwcDAIBTi7FmNouo6X58eB5X7xrAr161DacWYzi5wJgnIiIisgMOcrvME/MR7BwMQkSK3j/a64PbKTi9GG9yy4iaJ53N4ci5NVy+PYwDO/sBAPedXmpxq4jWNSqzgBkLRETUDTjI7TInF2LYPRQoeb/TIZgM+3GaK7lkY8fnokhnFS4cD+G80RCCHifuO7Xc6mYRAdDOMr/6b76LP/j8fXV93ZVYGr/0gR/iFf/+M6Szubq+NhERUTvhILeL5HIKZ5bi2D5QepALANsHAji9xEEu2ddjs6sAgPPHQnA6BBeO9+Kx2bUWt4pIc/MdxzEfSeEr90/j8Nn6xeVnfn4SJxZiuPuJRdz+6Nm6vS4REVG74SC3i8xHkkhlc9hWoWry9oEAV3LJ1s4saen4Owe0Amz7R3tw5OwaUzmp5XI5hdsfPYun7RsEANxxeK5ur/2jx+dwwVgIA0EPB7nUFY6ei2CVR8QRdSUOcrvImWWtYz/ZX2GQ2x/AUiyNSDLTjGYRNd2ZpRiGejzwe5wAgP0jISzF0liIplrcMup2x+ejWImncdPlk9g+4Mf9p5fr8rq5nMKDZ1bw1L2DeMqeARw8uViX1yVqV/efXsZz/vFHePlHfoZsjhOYRN2Gg9wuMq0PcicqruRq93M1l+zqzFIck/3rafv7R7Wjs+qZGkpUjftOaQXQrtwRxvmjIRw5W5/j3KaW44inszh/NITLt4VxejGOhUiyLq9N1I4+c9dJKAU8NruGgyc4qUONEUtlcNfxhZZngs1HknjwzHJL29BuOMjtIlN6iuZkpUFuP48RIns7sxTHtoKMhr3D2iD3xDxjnlrrwTMr6PG6sHe4B/tGQjg+H0GmDkWijpzTJnD2j/bgsm1hAMBDUys1vy5Ru7rn1BKetm8QbqfgB4/XL+2fqNDbPn8/XnnzXfjM3ada1galFF79H3fjhf/6E/z06HzL2tFuOMjtIlPLcfT6XAj53GUfZxSm4kou2VEupzC1FN+wN32s1weP08GJHWq5Y3MR7BvpgcMh2DfSg3RW4fRS7Ue6HZ+LAtAmdM7TMxeO6bcR2c1aIo3jc1E8efcgLhzv5QoXNcRCJInvHtLqG/zPvWda1o4Hz6zgcT0T7Yv3tK4d7YaD3C4ytSlFs5T+gBtBjzNfnIfIThaiKaSyuQ1p+w6HYNsAj86i1js+F8WeIa0gmpF1Y2w1qcX0cgI9XhfCAQ8Ggh6EA24cm6tPKjRRu3lcr5Z/8UQvLp3sw0NnVpDjvlyqs58cW4BSwLV7B/HQ1AoS6WxL2mHUbrh27yDuODzX8tTpdsFBbheZWo5jMuyr+DgRwWS/vy4dK6J2M6/vQxwOeTfcvmMggJOLXNmi1okmM5hdTWDPsDbIndCv1/W4Fs+uxjHWp72eiGDPUBDHOcglm8pX0B8M4OKJPqwlM5hin4bq7JHpFXicDrzymh1IZ1XdaihYdWhmFf0BN5578RgWoinMriZa0o52w0FuF9EGueX34xomwn5Mr/ADgezHGOQO9Wwd5J5a4Eoutc4T89okyx59j7gxKJ1err3DMrOSwHjf+iTn3uEeHD3HSR2yp6mCQpvGpJHx90VUL49Or2L/aA/OHw0BAJ5YaE2MHZpdw4Xjvbhksg8A8NAZ1lsAOMjtGquJNNYSmYrHBxnG+/x16VgRtZv1Qa5nw+07BgJYTWSwEuOZitQax/ODXK1T7nU5MdTjwUwdJhxnVxIY7V0f5O4aCmI+kkQsxaPiyH6mluPoD7gR8Ljy6f8c5FI9KaXw6PQqLhrvxc7BAETQsuyYUwtR7B4K4sJxbbD92CxPigA4yO0a65WVK+/J1R7nw2I0hXiqNfsLiBplfk07C3doU7qyUXCNxaeoVYw94TsG1q/TY30+zKzUNuGYyeZwbi25YSXXqC4+xdoLZEPTy/H8pP5wyIser4vp+VRXK/E0FqIpnDcags/txEjI25JtfvFUFkuxNCbCfgQ8Loz3+XCiRSvK7YaD3C6RH+SaXMk1ivIwZZnsZj6ShMflQMjr2nA7j86iVptejiOsrz4Zhnq8WIymanrd+UgK2ZzasJK7TY93FhgkO5pejmOiT+vHiAh2DwXzmRJE9WBMPhr95dFeH2ZXm3/2uJHpY0xi7hwM4ARjHQAHuV3DGKxOmCg8pT2uflU9idrJXCSJ4R4vRGTD7cZ+9XqkhhJVY3YlgfG+jRORg0EvFiK1dZyMFP2RguyF7fqE55klTuqQ/cysJDZU0N8xEOCEDtVVfnCp96tHQj6ca0HBJ2OwbXx27B4K4iTriwDgILdrTC3F4XE5MBT0Vn4wCjr83JdLNjMfSW3ZjwsAvX4XAh4n96JTy0yvJDDRt3EicrDHg4VoqqYjIYyV4IHgetwP9XjhcTnY8SfbSWVyWEtkMFgQ75P9fkwtx3m0CtWN0VcwMgbG+rwtqWq8Psg1VnKDWIimsJpgfREOcrvEGb2yssMhlR8MLe1CBCy5T7Yzv5bcUlkZ0FLaxvt8XMmllplZiedXBQyDQQ+SmRyiNdRHWIppg9z+gk6/wyHYFvZzkEtt5xsPzuCmf72z6kyy5SLxPhn2I5XJYT5SW+o/kWF2JQGnQ/LHEY6GfFiOpZt+Vm4+U6dXa8dOo74IV3M5yO0WU0vmjw8CAI/L0bJN9ESNNB9JYrDISi5gHJ3FlVxqvngqi+VYeku6srH6WkvKcn4lN7Ax7if7/UxXprbzz987jAfOrOCWg6erev5CkcwFI3WZE/dUL9MrcYyGvHDqi0dGzYO5tebuy12MpuB1OeB3OwGsx3qtBQvtgIPcLjG1HDe9H9cw3sezcslelFJYjqc3zPAXGu/zcWKHWmJ6U/EQg5F1sFBD8amlaAoOAXr97g23b+sP4DRXcqmNpLM5HJ/Tiubcd2q5qtdYKjLInWSdEaqzmeUExgsWj8IB7fpqZM40y2I0hcGgJ19nxMgGYlYaB7ldIZnJYm4tafr4IMNkmGflkr0k0jmkMjn0bersGybCfsxHkkhlck1uGXU7o/7BlsJTPcZKbvUdp8VYCuGAJ7/iYNjW78diNNX09DqiUk4uRJHJaftmqz0GZTFWZJDLI7OozmZW4hgrmJQ0Js+XYs3dC7sUTW2YuB8KeuF2Cvvv4CC3KxidJ7PHBxkmwj4WaiBbWYlrHz5hf4l05T4/lALOtqB4BHW3UhXwjY76YrT6FLilaBr9ga0TO8Z7cXWL2oWxinvdviFMLcWRzVnvfxjp+f0F6fm9Phd6vC6mK1NdKKW0Ct6Fg1z9GrtcxUquUqrqQlEL0dSGCR2HQ7RMTMY6B7ndwLioW01XntALNdSSJkfUTpbjWiyXWskdZ6efWsSYjBzblK4c1jvqyzWsDixu6gQZjKqgnPGndnFW38949a4BZHKqqpRLY5AbLpjYERFMhv0c5FJdLMXSSGZyGzJvjGv1UhV95vd96zFc8zffxfG5SBVt2Xp9ZxFNDQe5XcC4qG+zmK48wWOEyGZW9IFCuMiqFrCeKsq96NRsMytxDPV44HU5N9we9DjhdEhNx0EsxVIbVrUMPA+d2s3cagIOAS6Z7AUAnKuiiM9SNIU+vxtu58Yu7kSYNReoPow4KqyhEPYbe3KtX6s/+dMTSKRz+OHjc5afuxjZen2f4HZDABzkdoXp5ThEtq4QVDLJaoRkM8t6unLpPbnGSi4/HKi5ZlYSW/bjAtoKVK/PlU+1r0aplVweFUft5uyqdsSbcSxLNXvRN6dvGsbDflacpbqYNc6mLSg85XI60OtzWU5XjqUySKS1OiCHz65Zem4qk8NaMrMl3ifCPpxdTVSV7m8nHOR2genlOEZCXnhc1v53c5af7GalwiA34HGhz+9mmg813cxKfEtlZUOf342VeKaq11VKYUkvPLUZj4qjdnNuLYGRXi8GjariVRydpWUuFNmD3udjoTWqC6OPMLHpmt0f9FheyS0shmb13PJiZ0IDWlZaJqeafpxRu+EgtwtMLyfyA1Yr+gNueF0OdoDINiqlKwPa5A5T9KnZZpYTFQa51a3kJtI5pLOqbEVxpudTuzi3lsRoyIdB43zoKvY3rsTTRSd1jEwJruZSraZXEnA5JD8ZYwgHPJaPEDqtn1U+GfbnvzfL2May+fqez0rr8ms7B7ldQDsj1/og1yjU0O1/JGQfK/E0nA5Bj9dV8jETfT6mb1JTrSXSWEtmNqS+Fer1u7Fa5SB3Te8EhXzFY557t6idLES0VGOf24kerwvzVazkriUyRePdmERipg7VanYlgdFe35Zj2foDbsuDXGPS5cCufsysJCydaLKa0DJ8Nse7MaEz2+UTOhzk2pxSClPL8fz+Wqsmwn5MsQNENrEc1wqSGIemFzMe9nGmn5rKiLdSK7m1DHJLdYIM2nnoPCqO2oO2CqutSg32ePKVkq0oOchlMU2qk+nleNETS3p9bqwlrG0tMaox7xvuQSqTQyxlPp3eeK/eLYNcY0Knu2Odg1ybW4imkMrktuwbMIvVCMlOlmPpkmmbhvE+P1biacRS1e2BJLJqOn/MW/HJyFrSlY10tl5fiXTlPh+SmVxVgwmiekplcoins/lYHQh6LBeeUkphLZFGj3drvHMll+plZiWBsSKFAkM+FyIWB7mL0TR6vC6M6vFp5VpsTH6GNl3f+/xu+NwOzHZ5rHOQa3OVOk+VTIT9mFtLIplhoQbqfCvxyoPcyTDPDqXmqriS63NjNZGuarU1P9PvL76SO854pzaRLwyor+Rqq2LWJneSGW0PerGVXJ/biYGgB9NdvrpFtVFKYXYlUXTxqMfnsr6SG0uhP+jGgHHOroV057USmToigvE+VhPnINfmjEHuZH/1g1wAOLvS3RXayB4KU+FKYVVxarYZ/Zi30d7ShafSWYV4FVVh1/fkFo97HhVH7WJz9ftevzufbm/WeuZCiUmdPh9mGOtUg4VoCqlsruikZK/PjVQ2Z2lhaDGawkDAk6+QbGUlt9z1fazXxz25rW4ANZaxn7baPbnsAJGdmFnJNT64OMilZpleSWAk5IXbWfwj2YjZalKWS830GzipQ+3CiO9ePd5DPpflldxIPt6LX+e5ukW1MvZ0F0tXNopaWlnN1VZyPfmzbq2u5DoECHqcW+4b72N9EQ5ybW5qKY6Ax1mxY18KO0BkJ8uxNMIV/hbG+nwQAVPaqGlmVxL5apjFGKnGq1WclVtpJbc/oO3d4jWeWm1180quz43VeMZSmn7lSR3WGaHaGCeOFCs8ZcSdlX25i9EU+gOe/NnOi1HzEztriTRCvuLFNMf6fDi7mkAu171FBVs6yBWR54nI4yJyVETeXuT+14vInIjcr3/9diva2cmm9eODylWTLcdY1eJKbnUY4+0jl1NYTVReyXU7HRgJedkRMokxXrvpleKVOg21rOSuxkvP9APa3i2elVsa47t5Nqcrh3wuPfUzZ/o1jEFuqWPixvp8WE1kEE2ysKCBMW7NbL6GQn1Wcpdj2jaqXp8bDgGWLa7klprQGe/zIZNTmI9273bD0odFNpiIOAF8CMBzAJwB8AsRuVUp9eimh35BKfXmpjfQJrTOU3WpyoBWqGGox4upJXaArGKMt5dIKgOlSq9oFZoI+1mB0wTGeO2UUphZTuD680ZKPsboOEWS1aQrl57pN0zyqLiiGN/NVWxPLqDts/W5i0/SbFYpc2FCH5jMrCSwb6SnpvbaAWPcuumVODxOBwb19OJCRtytmbxW53IKkWQGIZ8bDocg6HUhYmECZjWRKRnrYwVn5Y6EqjthpdO1ciX3GgBHlVLHlVIpAJ8HcFML22NL08txTJZZITBjW78fp5didWpRV2GMtxFj5r6nxKxnoYk+P6vNmsMYr9FKPI14Olt2JXd9kFtN4anSM/2GiT4/i/EUx/huoi2DXJ/1NP1K6co8RmgLxrhFM8sJjPZ54XBsnTg04s7sSm5MLyYY0q/xIa+1I4i0ScxKsd69fZlWDnInAZwu+PcZ/bbNXiIiD4rIf4vI9mIvJCJvFJGDInJwbm6uEW3tSIl0FvORVNVFpwzbBwI4w5XcajDG24gxyA2WSGMrZOzbqubIli7DGK+RMZlSbk9uTxX7vAzlZvoN42EfzvGouGLqFt9A98a4WSvxNAIeZ74Am3FerpXiUxXPhdb7QzOcxDTwGm5RuRoKVvfkGo8zrvE9PmsruWuJTMlK4mP6IPfsavfGersXnvoagF1KqcsA3A7gk8UepJS6WSl1QCl1YHh4uKkNbGfG7E0t6cqAtpI7vRxHtos3rzcQY7xJjFWwHm/ltLfxPj+SmZylUv5UEmO8DGNFabzMSq4xMVPNPsLVMjP9Bh4VVxNT8Q10b4ybtbqp+r0Rt1aOEVpLlM/YGe01Cgty4t4CXsMLTK/Ei56RCxTuyTU3MWNsQTGe12M5XTldchJzIOCB2ylcyW2RKQCFs0Hb9NvylFILSinjU/c/ATypSW2zBWMfba2D3O39AWRyCrNdPBtUJcZ4G8nPmHrN7ckFujvNxyTGeI2MKt4TZVZygx6941TFIFeb6S8f8zwqriTGdxNtPuItvyfXQsG1SDKDoMcJZ5FUUgDwuBwY6vFyJXcdY9yCXE7h7Gqi6PFBQEHWjclr9eZCaUGvy1LRqnLbURwOwWiXn5XbykHuLwDsF5HdIuIB8EoAtxY+QETGC/75QgCHmti+jje1rO2jrTVdeVu/9vwzi9yXaxFjvI1E8unKlVdy2ek3jTFeo+nlOFwOwXDIW/IxTocg4HFWtZK7lkiXTGcz8Ki4khjfTRRNZTZURV5fybV2pEqlugsTfT6u5K5jjFswH00inVUlayh4XU54XA7TA9XIplohIYvpyrFUpuwWLO2s3O6N9ZZVV1ZKZUTkzQC+DcAJ4GNKqUdE5D0ADiqlbgXwVhF5IYAMgEUAr29VezvR6cU4nA7Jbz6v1vaBAADgzFIcT65Hw7oEY7y95AtPmdiTa6SOshhPeYzx2s0sxzHW5yu58mTo8bqqHORWLjxlfEZwkLsR47u5osnshlg1rtUxCwXX1kzsQR/r8+HYXLS6RtoMY9yaGRM1FHp9LtNZN5FNK7k9FgpPpTI5pLOq5PFwgFZh+aEzy6Zez45aNsgFAKXUbQBu23TbXxV8/w4A72h2u+zi1GIME2EfXM7aFuwnwtoeFlZYto4x3j4iFga5g0EPPC5HPpWUSmOM12Z6OWFqS0mP13zHyaCUdjxFpZUt7ag4D1e3imB8N088lcVY7/qkfEBP04+m6pO+aRjv8+POI/NQSpU9WqtbMMbNy9dQKLN41GMh5XgtuXmQ6zY9mRlPaZM/fk/5ldzvPJLo2lhv98JTVINTizHs0Fdha+F1OTEa8rHCMnW0iIXqyiKipbRxZYsabGq5dBGTQsEqVnKTmRyyOWWyojjPyqXWiqYyCBSsSjkdAp/bgVjK/Eru5pTnYibCPkRT2ar2uFN3W6+GX/qaHfC4ELO4khsqrK6cyiBnotBrLJ3R36/MSm6vD8lMDssx62es2wEHuTZ2uk6DXEA/K5d7cqmDRZMZuBwCr8vcZW8i7OcglxoqqxcxMbuSa/UIIWNwECwz028Y56QOtVg8lUVgU82EoMfa5E4smS3b6QfWU01ZfIqsml1NwONyYCDoKfmYoNdpemJm8+R7yOuCUuvn55ZjvEe5eO/2s3I5yLWpaDKDhWgqv5+2VjwrlzqdkbZpNmVnvM/ftR8M1Bxza0lkcsrUIDdo8WgJYH0feqVOP6BN6szwbGhqoWgqs2VCJmBhwFDqNTYzigYxPZ+sml6OY7zPV7YfEfC4EDOZYh9JZuBzO/JnQ1s5Ez2fruwutydXi/XZ1e6MdUuDXBEJikjlT0tqOWP/bD1XcmdXE8hkc3V5PaJm046WMF+GYDLsw1nGPDWQUb3bTAX8kM9laW8iULCSayJdeTLsRzSVxWqcKZzUfNmcQiKdg99T20puPJXd8hqbcSWXqjW9HC973BugreRGTU7MrCUyG441NFLtjfNzyzFzfc/HepdO2Jcd5IqIQ0R+TUS+ISLnADwGYEZEHhWRvxeRfc1pJll1aqG+g9zt/QFkc6pr/1Co80WTlfdqFRoP+5FTwNm1ZOUHE1UhX8SkxHEUhYJep+V0ZWNQbHYlF+CxWdQa8XTx1PqAp4qV3ArX+ZGQFw5BVx+tQtWZXk5gsr/8INfvduVXWSuJJjMbjjVcH+RWfr5xfS83qTMc8sLpkK49K7fSSu4PAOyFVlVtTCm1XSk1AuA6AHcBeL+I/HqD20hVOKXvn93eX7+VXADcl0sdK5rMmjoj15Dv9DNNnxrE2ANrbk+uG1ELR6kA60evmC08VdgmomYyCvVs2ZPrNZ/BYKwGV5rUcTkdGAn58kWEiMxIZ3M4u1a5hoK2kmuyQnI6uyHd2Kd/bybdOW5iT67TIRgJebt2garSJ98NSqkta+ZKqUUAXwLwJREpfyAZtcSZpThCXhfCgfr879kxqA2WTy7GcG1dXpGoudaSGfT5zf89GFkQpxZjuGb3QKOaRV1sejmBkNeF3grnegJAj9eJVDaHZCYLr8vcZI21lVzuU6TWiZbosAc8TpxbNZdNY6wGm4n38bCPK7lkyexKAkppW5nK0aorm5uQ3Jxeb8SumZXgfOEpd/mh3Fifjyu5xRgDXBG5YfN9IvK6wsdQezm1GMP2gUDdzsWa6PPD43LgiXkeoE6dSUtXNr+SOxn2wyHAyQXGPDXG9HLcVKoysL4aa2U111gNMLMXfSjohcfpYLoytUQsPyGzMVaDHvMrufnVYBPxPsHCgmSR2cybgEebkEybqOcRT2+sBm58byZFP24iXRnQKix364SO2cJTfyUiH9YLT42KyNcAvKCRDaPanJiPYudgfVKVAcDhEOweDOL4HDv81JmiFgtPeVwOTIT9OLnAFH1qjDNLcVNFp4D1vVpWivAYA+LNKaDFOByC8TBTOKk1Sh2HYqW6cjRfiMfESq7e8Wc1cTJrysIgFzA7UN2YruyvYiW3UryP9vows5Loylg3O8h9JoBjAO4HcCeAzyqlXtqoRlFtUpkcTi7GsHe4p66vu3soiOPzkbq+JlGzRBKZfHl+s3YNBnGS+9CpAZRSOLUYw87BoKnHG4PcNQvFp6ys5AI8K5daZ32QW2Ql1+TEjhHv/grpm4BWWDCRzmE5xmREMmfaZDV8I+vG1L7adBb+gpg34j9u4pxcY1LHV2H7ynifD7FUFmsWj6CzA7OD3H4A10Ab6CYB7JR65cFS3Z1ajCGbU9gzbK7zZNae4SBOLcR4pAp1HKUUoilr1ZUBbS/6KaYrUwMsxdKIJDOmzzI3JmisHCNkrOSWO0ex0ETYz0EutYSRarx5VSrgcSGZyZnqd5hd2QKAiT7uQSdrppYTGAx68sWhSjFWcs1sLdFWcteHYn63tXRlv9sJh6P8cGxMP0aoG/flmh3k3gXgW0qp5wG4GsAEgJ80rFVUk2Nz2mprI1ZyMzmFM6w2Sx0mns4ip8xVmS20cyCApVgaK3HO9lN9GXu9d5oc5Bqxa+UYoVgqg4CncifIMBn282xoaoloiSI6xoA1ZmZly8Ke3PEwz8ola6aX46Yq4edXY80MVLdUV3ZAZH2/bTmxVNZckTV9Qqcb96CbHeTeoJT6GAAopeJKqbcCeHvjmkW1MPbNNmIlFwBTlqnjGAMDy4NcfV/7Ke7LpTozjnnbYbJ2QshIV7ayJzeVNdXhN0zwbGhqEaNTv3n/uBG/ZqrVxi3uyQV4Vi6Zpw1yTZxpbqzkmjwGqDBdWUTgd5vbh765MnMpRqx343GIZQe5IrILAJRSpzbfp5S6QzTbGtQ2qtLxuQhGQl6ETBxLYcWeoR799Zm+SZ0log8MQlbTlQe0iZ2Ti4x5qi9j4sTsWeb5fV4WBrmxZKaqs6GZskzNli8atXlPrtf8gKHUanAxQz1euByC6S5c3SLrlFJ6ocDK1+uAyT25mWwOqWxuy3aSgMdpKnMhlsqaqrcw3ueHx+noyn5Mpd/O34uIA8BXAdwDYA6AD8A+AL8E4NkA3gXgTCMbSdYcm4vUfRUXAPqDHoQDbhznMULUYYy9MdWu5LLCMtXbycUYRkJeUzPxwHrnP9LAlVzj/EcOcqnZYskMRLR0zUJWVnJjJVaDi3E6RKs6y1gnE+bWkoins9g1ZGKQa3JPbiKjbQvxezbGvN/jNJXqHE1lTH1+OB2CbQP+rsxIK/vpp5R6mYhcBODVAH4TwDiAOIBDAL4B4G+UUpwGayNKKRybi+JXLhtvyOvvGQri6DmmK1NniZQoalJJ0OvCUI8XJzixQ3WmVVY2f8yb0XE3e5yK9thMPnXOjHG9QAnPyqVmi6WyCLid2FzT1ErqZ6zEanApOwYCONGFHX+yzogTM9XwAyaPAcpXA98Ur363uUFu3OSeXEA7KaIbY73ilUAp9SiAv2hCW6gO5iJJrMTTdS86ZTh/rBffeHAaSqktH0ZE7Wo9Xdl6Cv/e4SCOznFih+rr5EIU1+0bNv14t9MBj8th+ZzcXr/5mA96XQgH3Di9yEEuNVd0095Eg9nUT6D0anApe4aD+NoD7M9QZSf0QoG7TExMGpMslSZmEil9JXdTurLf4zKdrhwOmLu+7xwM4O7jC10X62YLT0FErhWRXxOR1xpfjWwYVefQzBoA4MLx3oa8/kXjIawmMl1ZpY06V7TKlVwAOG80hKNnI115kDo1xmoijbOrSewdsbatpMfrsnSEkNWVXECryn+MkzrUZPFU8f3jQQvHsURLrAaXsme4B6uJDBaiKWuNpa5zciEKl0MqnpELmM+6iaWNc5037cl1O01VV958xm45OwcCiKaymI90V6ybGuSKyH8B+ACA66AdIXQ1gAMNbBdV6dDMKgDgogYNco3Bs/E+RJ3AWMm1ek4uAJw32oO1JCd2qH6O6Vs+9o+ELD0v4HGa6uwboklre3IBLd65JYWardT+cWMl10wGQyyVzT/ejPyJESymSRWcWIhhW78fLmflYZPH6YDTIRWzD4yU5M0pxwGPuerK0aT5ScydQ3oRzYXuinWzK7kHADxNKfV7Sqm36F9vbWTDqDqPTq9iMuxHn8kUBqvOH9M6ZRzkUidZ35NrfZC7f1SL+cNn1+raJupeR/RB5L4Ra9tKerwuS4WnYiVWx8rZNxLCYjSFhQiPEaLmMc503szoxJvp9Jd6jVL25k+M4KQOlXdyIWpqPy6gHQNkZkIyrqck+7akK5vfk2u2cOGuQWOQ2137cs0Och8GMNbIhlB9HJpZxYXj1lYHrAj53Ng+4M+nRRN1gqi+V8tKB8hwnj7IPXKWHSGqj2PnIvC4HNjeXzn1rVDQ6zK1N9FgtboyAOzXB95HuJpLTRQrUUTH6MTHTexRtJq5MNnvh8fl4IkRVJZSCk/MRbF7yPz2kqCn8rW61EqumXNylVKIpc0XntrW74fLITg+313XdbOD3CEAj4rIt0XkVuOrkQ0j6xLpLI7PRxu2H9dw4VgvDs1yJZc6RySZQY/HVVXBhYGgB0M9XjzOlVyqkyPnItgzFDSV+lYo4HEiYjJdOZ3NIZXJWd6Tu3+Ug1xqvliy+JmfRuqnmXTleNraHnSnQ7BrMMCVXCrrzFIc0VQ2P+FtRsDrzJ/bXIoxcbN5NVZLVy4f76lsDtmcMj2p43Y6sHsoiMdnuyvWzU55vbuRjaD6ODSzimxONWw/ruGiiV7cfugs1hJphHyNSYsmqqdIIlNVqrLh/LEePMaJHaqTw2fXcMX2sOXn9XhdmDW5N9xYCbCyRxEAxnp9CHldOMpJHWqiaIlUYxFBwMTKFmC9mjigbRl4eIrXdirN2Kp0/pj57SUBEynHxv3Fqisn0rmqnlvOeWMhPHRmxfTj7cDUNLJS6kcAHgMQ0r8O6bdRG7n31DIA4Kqd/Q19n6t29EMp4IHT3fXHQp0rmsqgx1f9IPeybWE8NrOGhImUOaJylqIpnFmK4+KJPsvPDXhcpo8QMlYCrK7kigj2j/ZwSwo1VTyVzVel3czsHsVYKoOAhU4/oBXpPLUYw2oibel51D2MLC5LK7kmrtXlVnJT2Rwy2dIDXWOV2ErNhQtGQzi1GLO05aXTma2u/HIAPwfwMgAvB3C3iLy0kQ0j6+49tYTJsB+jvb6Gvs8VO8IQAe45udTQ9yGql0gyW9NK7uXbwsjkFB5lwTWq0cPT2uTgpZPWB7k9JlLgDEbRE6sruYA2qfPw9AqyOR6bRc2hreQWj9Wg1/y5oaUGyqUYk02Hpnltp+Ien13DZNhvKXPRTIXkUquxRkZDuZg3jhgye4QQoK3kAsDhLqovYnZD0F8AuFop9Tql1GsBXAPgLxvXLKrGfSeXcOWOcMPfp9fnxvmjIdxzioNc6gzRZAY9VZyRazBSSx84vVyfBlHXemiq+kFu0KutDpg5s7nalVyjbbFUluflUlNkcwqJdK5kER2/yXNDY6ni+3rLuXhC2971CAe5VMLjs2s4b9RaJXxThadKVFc2/l0ueyG/HcVC5sL5+kr041209crsINehlDpX8O8FC8+lJphZiWN6JYGrdjQ2Vdlw1c5+3HdyCTnO9FMHiCQyljs/hcb6fBjt9XKQSzV7eGoFOwYCVR3zFvS6kMkppMqksRnyK7lVxP3l27UB+INdtn+LWsPo7Je6Rps9HzqatHaEEACM9Pow1OPlIJeKiqeyOHIugosmrNW6MZNiH09l4XVphdUKBUwcmxUrUZm5nB0DAYS8rq66rpsdqH5Lr6z8ehF5PYBvALitcc0iq358ZB4A8NS9g015v6t39WMtmeEHA3WESLK2PbmAlrJ8Hwe5VAOlFO4/tYxLt1lfxQXWV2XNdPjzK7lVZDDsGepBj9eFB88sW34ukVWxpJF6WXpPbqV05WxOIZnJVTWpc+lkL2OdijK2bVy53doCUsBjrrpysZhfH+SWXgnOpzpbGOQ6HILLt4dxn16/pxuYLTz1pwBuBnCZ/nWzUurPG9kwsubHR+YxHPLigrHGnZFb6On7hwEAP3z8XIVHErVeNJVBTw17cgHgmt0DOLkQw5ml7jpMnernzJKWcfPk3QNVPd/YX2um+FQ0Vf1KrsMhuGJ7GD9/YtHyc4msilUoohP0uCqmK9cyqXNg1wCOnItgIZK0/Fyyt/v0bXlXWNwKaHYlt1h1ZGOfbblCl1E93q1e36/aEcZjs6umCxh2OtMpx0qpLyml/kj/+nIjG0XW5HIKdx6Zw9P3D1V1Dmg1hnq8uHxbH37AQS61OaVUzUcIAesTO3fqWRNEVv3s+AIA4Cl7qsu4MSZqomb2JyaNTlB1e9Gv3TeIx2bXMLfGjj81VqUOu5kiPrEqVrYMT9mjTTpxUoc2u/fkMrYP+DHU47X0vKDHVbFCcqzESq4/vye3zHOrSFcGgCt39COnumcrStlBrojcqf93TURWC77WRIR5qm3i3lNLWIql8czzhpv6vtefP4L7Ty9jKZpq6vsSWZHM5JDJqZpXcs8b7cFIyIsfH+Ugl6pz17EFDAQ92D9irYiJIZBPVza/klvtXvTr9g0BAH56jPFOjVWpw+63MMitJt4vnQzD53bgbg5yqUA2p/Cz4wt48m7rk5JmKiQnSqzkWklXtj7I1U5HufuJBUvP61RlB7lKqev0/4aUUr0FXyGllLVd2NQwtz4wDa/LgWdfONrU973hwlHkFPCtR2ab+r5EVhgDgloHuSKCZ5w3jDsOzyGZ4Xm5ZE0mm8MPHj+HZ9SQcZNfyTWzJ7fCPsdKLp7oQ5/fna/3QNQosQqp9dpKbvmJnWgNmQselwNX7xrAHUfmTFUup+7w4JllrMTTeEYVC0jGdTdW5lodS2WLxmu+unKZAXKlv5lSwgEPLt8Wxg8en7P0vE5l9pzc/zJzGzVfJpvDbQ/N4NkXjtTcibfqksle7BvpwZfuOdPU9yWyIpI09mrV/vfxvy4bx1oigx91yQcE1c/PTyxiKZbGcy8eq/o1ghb35HqcDnhc1R2E4HQInn3BCL7zyCwndaihYsny+2n9HhcS6VzZ0xyq7fQbfvmiURyfi+LIOR6bRZo7Ds9DZD2rxQojo6Dsamw6u+X4IGB9oqbcnt54KgMRwOe2fn1/9gUjeOD0cldsRTH727m48B8i4gLwpPo3h6z69iNnMR9J4Vev3Nb09xYRvPiqSRw8ucTzFKltRfIrudWfk2u4bt8QBoIefPWB6Zpfi7rL1x+cgdflqGpVwGB0nCpV7QS0zlWgxph/wRUTWE1kcMdhruZS46yf+Vl8gGpUFS+/sqWv5FYZ88+9eAwiwDcfYmYaabU8vvbgNA7s7MdA0GP5+X4TxwAl0iUKT5lYyY3qqc7VZAX90gUjAIDvHTpr+bmdptKe3HeIyBqAywr34wI4C+CrTWkhlfXRO49jx0AAz9KDttle9qTt8Loc+PcfHWvJ+xNVYqR29nitn0u6mdvpwAsvn8B3HpnF7Eqi5tej7hBJZvDV+6bwK5dN1JRRYKx0mVrJTWZrOhsaWJ/U+eLB0zW9DlE5lQaoVs4NrTbmR3p9uHrXAL5835myK8bUHR6ZXsXRc5GqF5DMxmyxdGW/qUmdbNVZCxdP9GLvcBBf7IIszEp7cv9WKRUC8Peb9uMOKqXe0aQ2Ugk/ePwc7j21jN+6bveWw6SbZTjkxSuu3o7/uXcKJxeiLWkDUTnRCqlwVv3WdbuRU8B//vh4XV6P7O8LvziNaCqLX3vyjppexxggR0wMcmOpTNWVlQ1upwOvumY7bj90FseZrUMNEq1YeKpy6mcte3INr3nKTpxYiOH7j/HUiG738Z+cgM/twI2XVre9JGAyXblYzQSvywGRyunK1ca6iOCVV+/APSeXcGjG3jWEzaYr/1xE8qfXi0hYRF7UmCaRGbFUBn/9tUexeyiIV11TW8epVr93/T743E688ysPs2gDtZ21OhWeMmwfCOCmKybwqZ+dZJo+VbSaSONDPziKa/cO4iqLZy1u5nU54HRIxSI8gDZwCNQh5l9/7W54XQ78/bcfr/m1iIqJpbLa/kJX9Su5xqpXLYPc510yhsmwH/9w+2FkuZrbtU4uRPGV+6fwyqt3IBywnqoMmNtXm0gV35MrIvC7y5+zW2oV2KyXHdiGkNeFf7j9cNWv0QnMDnLfpZTKH6qklFoG8K6GtIgqyuUU3v6lh/DEQhR/86JLqi4sUi9jfT782fPOx4+PzONfv3+0pW0h2ixfXdlXv8Jsb3/+BfC5HfjjWx4oe2A7dTelFP7yKw9jJZ7GO55/Yc3nmIsIAh6n6erKwRpXcgEtW+ctz9qPbz48i1u5F50aIJbMIOB2wlEiI83MINf4m6g2hRPQMhfe/vwL8OjMKjN1ulQup/C/v/wQ/G4n3vTMvVW/TqWYVUohli49UA14nGWPHyq1CmxWOODB7zxzD25/9Cy+beMTUsyOjoo9ruYeo4g8T0QeF5GjIvL2Ivd7ReQL+v13i8iuWt+z060l0njbF+7HrQ9M409++XxcW0XVt0Z4zVN24sVXTuKDtx/G33/7MaTLHIDdTRjjrRetY3Vlw0jIh/e/5DI8cGYZb/jUQazE03V77U7DGC8uk83hvd84hK/eP423Pms/Lt3WV/lJJvR4XaarK9fS4S/0hqfvwYGd/fiTLz6A7z5q/2IlhRjfjRdNZfMpycUYcdyoarOFfuWycTz/kjG8/1uP4X/utf+eRYAxbshkc3jnVx/GT44u4J3/60KM9fmqfi1/hbNu01mFbE4VLTwFaMcIJcpO6tS+HeWNz9iLSyf78Ce3PICDJ+x5RrTZq8FBEfkHEdmrf/0DgHtqeWMRcQL4EIDnA7gIwKtE5KJND/stAEtKqX0A/hHA+2t5z06VzSk8PLWCf/neETzrgz/CNx6cxp8973z8/i/ta3XT8kQE73vJZXj5gW340A+O4bn/dAc+/pMncGwu0rVFHBjj7WEtoQ9y69ThNzz/0nG878WX4mfHFvDsD/4I//bDozh6rrvinTG+1UIkia/cN4UX/dtP8NE7n8DrnroTb312/a7VQa8LURPpyrFUpm770D0uB25+7QGcN9qD3/7UQbzlc/fhp0fnbZ/FwPhujniFWDU68+XiPpbKIlBltdlCIoIPvOxyPHn3IP7olgfwu5++Bz89Zt9YZ4xrW0pue2gGL/7wT/HZu0/hTc/ci1fWuA1w/Qih4nFjpNcXS1cGjLOhG1N4yuBxOfDvr3kSBns8eMXNd+Hdtz6CR6ZXbJWqb/Y39BYAfwngC/q/bwfw+zW+9zUAjiqljgOAiHwewE0AHi14zE0A3q1//98A/lVERFWx8fOB08v4yI+OQSlAQXu69j3y36PgduNfxlupErdvfLxa/77YbetvseX2Yu+TzSnMR5I4t5bMB90zzhvGH96wH1fu6Lf6K2g4j8uB97/kMvzyRWP4x+8exv/5mva/0ud2YCTkw2CPB16XA26n9uVyCMx+HgnMf3DV+Bm3wXtuugTDIW+1T29qjH/u56dwx+G5/M8vECD/PfIf/tr3xW9HwXM3PqbI7QVPkPW30r8vdrv2XcDjxHifD9sGAnjSzn70+mqvelxONJmB3+1sSHG2V1y9AxeO9+J933wMf/etx/F333ocXpcDk2E/Qj4XenwuuBwOOARwiMDhEDhk4++xlUZCXvyfmy6p5SWaGuPv/MpDmF9Llb2Gl7uurj924/V787V7y+uXuW4br5RI5zC7EsdSTFvV3zUYwP971ZV4weUTVn/MsoIm05Wjyfqt5ALAQNCD/37Ttfh/3z+Cj//kBL72wDREgIk+P/qDbgQ9Lvg9TjhE8tcMKfy+RTF//fnDeMXVVXdYmxrfZ1cTePetj0B7r8rX8MJra7Hrd8lrcpFrd9FrvP7aA0EPJsJ+XDzRiwvGQjUPJDerlHXgN7G/sdJqsBVBrwuf/M1r8OEfHsPNdxzDNx+ehUOAyX4/+gMehHwuOB0OOEU7T9oh2lerrum//0v7cMlk1ZkiTY3x7x06iy8e1FbIjetpsWt44bW6VP97cz97y7+NJxR9L4VMTmF2JYG5SBJKAZNhP/75lVfgpismrf5YWxgrtCUHuRXOdfa7nWWrK8fLpDpbMRH246u/fx3+722H8Om7TuITPz2BgMeJsV4fBoIeeN0OuBwOuJ0Cp0Ms9cVrVfj39G+vvqqq646pK4JSKgpgSwpDjSYBFJ5LcAbAk0s9RimVEZEVAIMANhzaJyJvBPBGANixo/iHWTSVyReJ2fxhW7TzX+KCj2KP3fQaxsOk4NNJjMc5jNul4MOmeAfA6RCcNxrCeJ8P+0Z68OQ9Axjv8xf9+dqFiOCGi0Zxw0WjOHougntOLuLw2QjmI0ksRlNIZnKIJDPIZJXplGYrl9H1Lmd91Jh23dQYn1tL4qh+kP3mDwjjxmK3b+nUF/wKCwcJRR9bcDuK3q42PUYbdGb0SRuHAM+6YARvffZ+XLYtXPTnqlU0lanrftzNLtsWxmff8BScXozhp8fmceRsBDOrCawlMogk0siqLHI5hZxSyCntd5JrkwJtdVidaGqMn1yI4exqouI1vNz1e+M1ev35mzv6GwYXJq7bAwEHrtoRxvaBAK7ZPYDLt4UbMrESNJmuHEvVZ09uIZ/biT997gX4/V/ah58cXcDDUys4uRDFaiKDSDKDxWgKSgE5td6hVPnvWxPzNXT+gTrGN1A5xlOZHI7NRTb87oDNnfXi1+8tE+xFJ+hLDypKvXYup/LF+wBtIPCm6/fi1dfsKLmH1qp4hSI6pgpP1aGaeCGPy4G33bAfb3zGHtxxZA4PT63g1GIMK/E0VuPr1/Vs/treumt6ud+LCU29hq/E0zg+H9kw6aI/V/sv1q/hWybLi0zybJ342dTvdhj/dGx5XadDcP5oKH/NftLOfrid9alx43AIfG5HyXRlYwDr9xR/P7+nsYWnCvUF3Hj/Sy/DHz/3PNx5ZB4PTa3g3GoSC9EkEukcMtlMPr26Wer1eWGq5yci5wH4EwC7Cp+jlHpWXVpRI6XUzQBuBoADBw4U/c1cu3cI3/nDZza1Xd1u30gP9o30tLoZtmAmxt/67P1467P3N7Vd1cjpGQpH5yL48ZF53PKL03jRh36CP3/eBfidGgo9lLKWyNStsnI52wcCeMVAayuddzIzMf5fv7W579V9Ah4XFqOxso/J5ZTWCWpQ3Ac8LjznolE856LRhry+XVWK8e0DgbbspyTSWUwvx3Hw5BK+ePA0/vIrD+N7h87iw69+Uk3FbwzRVPlrtJnjWOrZ6S/k9zjx3IvH8NyLqztKppuYuYa/+KptePFV1Z0922mCHlfFldxSe3L9bifmI6mSrx1LZuB31/f6PhLy2e7/j9nf0BcBfATAfwKo18aEKQDbC/69Tb+t2GPOiIgLQB+AhTq9P1GjMcaLcDgEI70+jPT6cO3eIfzu9Xvxji89hL/95mNwiOANz9hT1/eLJuu3N5G2YIw3WY+3/F4tYH2VoN4ruV2I8Q1tBX/PcA/2DPfgZU/ahk/ffQp/9dWH8bbP34d/f82Tak5fjiWzGCmzLcjMcSyNGuR2AcZ4g5RbjY2nM/pjig/DAh4X4ul40fsqVWamdWbX5TNKqQ8rpX6ulLrH+KrxvX8BYL+I7BYRD4BXArh102NuBfA6/fuXAvh+NXsAiFqEMW5Cr8+Nf3nVlbjx0jG871uP4d5TS3V9/Wgy25SV3C7FGG+ygIl0ZaNAT6NWcrsI43sTEcFrnrITf3HjhfjOo2fxxXtqr0AcS2fKFgbU6ngIomUHuZm67kHvIozxBilXPCqe0rbClauuXGqAnMzkoBQQ4OR9RWYHuV8Tkd8TkXERGTC+anljpVQGwJsBfBvAIQC3KKUeEZH3iMgL9Yd9FMCgiBwF8Eeo/75gooZhjJvndAje/5LLMNzjxbu++khdKxRHks1JV+5GjPHm6zFRXTmW5EpuPTC+S/vNp+3GgZ39+LtvPV7z3v5YMluxw+53OxGvkK5cj9TpbsMYbxy/p/S12ki9L5mu7HGULDxlDJwDJZ5L68z2/IwZnD8tuE0BqCmvUCl1G4DbNt32VwXfJwC8rJb3IGolxrh5IZ8bf/a88/FHtzyA7zw6i+ddMl6X140kM3U9I5c2Yow3V9Dj0ouB5OAqUSQlv5LLla2aMb6LczgEf/zL5+NV/3EXbjl4Gq996q6qXytqYhU2UGZ/I6Cl6HNSpzqM8cYIlk1XNgpPlTpCyFVyD3qM13fTTK3kKqV2F/mq78Y5Iup6N10xiW39fnzsJyfq9ppRDnLJRoz95bEyq2fGYIB70amRnrJnAJdt68Pnfn668oNLyOYUEulcxf2FAa+zbMxHk/U7QoioHsqlKycqDHJ9bicS6VzRrDbjNZm5UJmpQa6IvLbYV6MbR0TdxekQvPapO/HzJxbzxyHVKpLMIMRBLtmEMWFTbl+ucR9n+qmRRAQvuWobDs2s4tDMalWvsV4krdJKrhOxMjFf7yOEiGrlL7saWz7l2IjlRGbrIDn/XMZ7RWb35F5d8PV0aIdCv7DcE4iIqnHTFZMQAb7+4HTNr5XO5pDM5LiSS7ZhdGzKDXK5kkvN8oLLJ+B0CL7x4ExVz8/vTay0kusuna7MarPUjoLlCk9VWMk19uoWS3dmurJ5ZtOV31Lw9QYAVwHgAahEVHejvT5cs2sAX6+y01TIGAhwkEt20ZNfyS2XuqnHPTtB1GADQQ+u3B7GHUfmqnp+vkhapcJTHmfJQjz5arOMd2oj5Y4QSqSyEAG8ruLDMGPwWyzm41zJNc3sSu5mUQC769kQIiLDcy4axdFzEUwtFz8nzqyI3tlnujLZhdGRN7OSy04QNcP15w/jwTMrmI8kLT83mq8yW/4aHSxzPvR6ej7jndpHwKPtIy922lIslYXf7Sx5xnT5lVxe380yuyf3ayJyq/71dQCPA/hyY5tGRN3q6fuHAQB3Vrk6YIhwJZdsJr+SW6bSrDFwYNxTMxjX658dW7D83LjJ1Hq/21VyVYyFeKgdBTwuZHMKyUxuy33xdLbk8UHac0uv5JpN8acKRwiJiFcplQTwgYKbMwBOKqVqPwGciKiI80Z7MBLy4o4j83jF1Tuqfp31dGV+GJA9GLFcdiU3mYWjTCocUT1dNNELn9uBe08t4QWXT1h6bjS/KlW58FSpM0eNgQBXtqid5AeqqSx8mwa08fTW2woZA+Bi2Qsxk38zVPmc3J9B23/720qp1zShPUREEBFct28IPzo8B6VUyZSeSiL6fq8ermiRTeSrK5fo8APrZ0NX+3dDZIXb6cDl28K49+SS5efGTKYalzuOhemb1I6MeIyls+jfdF88Vb5QWrk9uYx38ypN83pE5NcAXCsiL9781YwGElF3unJnPxaiKZxZqn5frrHa1ePjIJfswcwRQrFUhkWnqKmu2tmPR6ZX8+d/mpWvBF4hXv0eJ1KZHLLFzg3lkVnUhoxzm4sdfRVPZ8umG/s9pffkxlPM1DGr0hXhTQBeDSAM4AWb7lMA/qcBbSIiwpXbwwCA+04vY/tAoKrXiCRYZZbsxThXMVKuunIqiwBT9KmJLt8WRian8NjsGq7Qr91m5I9DqRCvxjU8lsog5HNveg2ubFH7CXpKpxwXS2EuVKnwVMDDTB0zyvb8lFJ3ArhTRA4qpT7apDYREeH8sRC8LgfuP7WMF1rc52UwCk8xXZnswuEQLXWzzEpuNMmVXGqui8Z7AQCPzaxaHOSaG6AWrmxtGeRyTy61IX+5QW46i4Ggp/JzSxSeYtEpc8yudX9eRN4pIjcDgIjsF5FfaWC7iKjLuZ0OXDrZh/tPW9/nZeA5uWRHQa+r7J7cWDLLYmvUVNv6/Qh6nHhsds3S84zCUz5X5T25hY8vZEz4+DmxQ23ESJ+Pp4ukK6fKV1c27kuUWMkNcpBritlB7scApABcq/97CsB7G9IiIiLdJZN9eGx2Dbki+7DMiCQz8Lgc8HDvCtlI0ONEtGy6MldyqbkcDsF5YyEcmlm19LxYMoOgxwmHo3zqZb6IT5HJnfV9vez4U/sw4rHYtbrSEUL5dOUSK7ncf26O2Z7fXqXU3wFIA4BSKgaAyeBE1FAXjIUQS2VxeilW1fMjyQxCXMUlmwl6XRUKT2URYNxTk10w1ovHZteglPlJyajJWM2vipVI/QR4bii1l0rFo8rFq8vpgMfpKHmEEFPzzTE7yE2JiB9asSmIyF4AyYa1iogI2r5cAJZT4AyRZIaVlcl2gp7y6cqRZAY9TFemJrtwPISVeBpnV813D7VK4JVjNVBmf2MslYHTIfA4mbFD7SNQUCxts0oruYA2SC5WrdzsxBCZH+S+G8C3AGwXkc8A+B6AP2tUo4iIAOC8UW2Q+3i1g9wE0zbJfoLe8unKsSTT2aj59gz1AACOz0dMPyeazJqKVX+ZdOVoMouA28lqs9RWAiWKRymlKh4hBGgpy0XT85PmJoao8hFCAACl1HdE5B4AT4GWpvw2pdR8Q1tGRF0v6HVh+4C/+kEuV3LJhgJeF6KLxVP4czmFWJqFSaj5dg1pR72dmI/h2r3mnhNLZUwVSVtfFSue+skjs6jdeF0OOEQrBFgomclBqcrp9QGPE/F0bsvtsQqpzrTOVO9PRL4G4LMAblVKRRvbJCKideeP9uLxs9UPcsd6fXVuEVFr9XhK78lNZLJQihXFqfkm+vzwuBw4sWC+mxhNZRH2uys+rmy6ctrcajBRM4kIAh7Xlpg19uhWSlf2uZ2IFy20xgw1s8ymK38AwNMBPCoi/y0iLxUR9hyJqOH2DgdxaiGGbBUVlrmSS3akFZ4qnq5snA3NPVvUbA6HYOdAAE/Mmx/kxpJmV3LLFfHJVBwwELWC3+PccoRQvlBahZjVVnJL7cllvJthapCrlPqRUur3AOwB8O8AXg7gXCMbRkQEALuHgkhlc5hejlt+biSR4YoW2U7Q60Q0lSlaxdZIjWO6MrXCrqEgTlgZ5KbMrcKWS1fW9vUy3qn9FDvuzYjhintyPc4tkzqZbA6pTI4ruSaZLkWnV1d+CYA3AbgawCcb1SgiIsOuoSAAWFodMPAIIbKjoNcFpYqfoWhUXebkDrXCrsEATi7GTJ9tHjVZXdnpEHhcjuKFeNKsNkvtyV8kXTlhciXX53Zuea5RxIqTOuaYGuSKyC0ADgF4FoB/hXZu7lsa2TAiIgDYU+UgN5XJIZnJoYedH7IZY1BQLGU5ml/JZdxT8+0aCiKVyWF6xVzmTSxpfoAa8Gzt9ANaunKA6crUhgJF0pXNruQGihwhZGTqcA+6OWZ/Sx8F8CqlVOkzC4iIGmA45EXQ47Q8yDUK83BPLtmNsUobTWYwHPJuuM9YyeWeLWqF3YPapOSphRi29QfKPjaVySGVzZlOrQ8UWdkCmK5M7SvgcebrJBjiJldj/UXifT1Th/FuhtkjhL4tIpeIyEUAfAW3f6phLSMiglahcNdQ0FLFTmC9AA/TNslujFn8aNEzFLVOETMYqBUm+/0AgDMmaigY+w3NrkoFvK4tq2IATJ05StQKAY8T51aTG24z4t5XIfvAX6TwlHF9Z6E1c8weIfQuANcDuAjAbQCeD+BOABzkElHD7R4K4qGpFUvPMQa53JNLdtOTX8ktvSeXK1vUCmN9PogAU0uVB7lWV6UCRYr4ANqRKox3akcBjwuxLdWVM/n7yvG7txaeirHmgiVmC0+9FMCzAcwqpX4DwOUA+hrWKiKiArsGgzizFEcmu/Vg9FIiTFcmmzIGBcXOyjVu455cagWvy4nhHq+pavixlLnOvqFYpz+TzSGRzrHTT22pWIXkeErrx5g5QiiTU0gX9HtiKRaessLsIDehlMoByIhIL7Tjg7Y3rllEROu29fuRzSnMriZMPyeS0Ae57PyQzeT35BZLV9Y7Qez0U6tM9vsxZWKQmy+SZmEld/OqWDTF9HxqX8WPENJi2Ex1Ze3x689n9XxrzA5yfyEiYQD/AeAeAPcC+FmjGkVEVMgoYHLGRAqcIb+Syw8DspnCwlObRZMZuJ3acStErTAZNjnItbiSG/BuPY6F13lqZ36PC/F0dsORWvnV2AqTO8bfRWGF5XxlZu7JNcXsVaEXwMsA/BDAtwD0KqUebFSjiIgKbTOKmVQzyGW6MtlM+SOEMjxeglpqst+P7zxyFrmcgsMhJR8Xs3jcVaBIujKr6FM7M9KKE5nsesHAZAYelwNuZ/mJSL9Hu78w5mMsqGmJ2anejwIYB/D/AHwfwLtE5G0NaxURUYHxsFbM5MxSzPRzmK5MdlXYWdosmsoy5qmlJsN+pLI5zEeSZR9n9bgrrfDUxphnFX1qZ8aE5OaUYzPXaL/bVeS53JNrhdkjhH4gIncAuBrALwF4E4CLAfxzA9tGRARAK2YyGvKZqthpiLAAD9mUx+WAx+nId3gKsdIstdpkeP0YoZFeX8nH5fePmy08pad+FjImM1lFn9qRX4/tWDIL9Gi3RZNZU/vQjWOxCmM+nsrCIYCX21FMMfVbEpHvAfgJgFcAeBzA1UqpCxrZMCKiQpP9fsvpykGPs2y6HFGnCnq3rmoBWgcqwA4/tZBxVm6lScl8JXALK7np7MZqs1Gu5FIbMyYcCwumaX0TMyu5+iB30ypw0OOCCPs1ZpidCngQQArAJQAuA3CJiPgb1ioiok229ftxZtlaujL3aZFdBb2uooPcSDKDHpODBqJGMFZyKxWfWj8OxeSe3CKpn2ssPEVtrFjMxlIZU5MygSIrudFkxnR6P5kc5Cql/lAp9QwALwawAODjAJYb2C4iog229fsxs5wwfVau1tlnx4fsKehxFT1CaC2RRsjrbkGLiDQhnxshnwszFQa50WQGPrcDTpPZNsZgOFYQ91EOcqmNBQrTlXWRZNbUIHf9CKH1eF9LZBDy8fpulqmrgoi8GcDTATwJwAkAHwPw48Y1i4hoo239AWRyCmfXkvmVgnIiyQx6+GFANqWlK2/dk6t1gtjhp9aa6PNjarn8ueZG6qVZxVbFmK5M7Ww9ZjdOzEyGS+9V3/zcwiOEIkle360w+5vyAfgHAPcopbZOHRMRNVj+GKHFmPlBLtN6yKaCXle+uFohzvRTOxgP+zCzUiFdOZnNF9cxI1+IZ1O6ssfl4LnQ1JaKFY+KmTzmLVDkqLjVRAa9HOSaZjZd+QNKqbs5wCWiVslX7DRZfCqSYLoy2VfQs3VPbjan9AwGxj211kTYj5mV8iu5qxYnZIKeIkeqcFsKtbFiMWt2K5WRnVA4mRlJpNHLSUzTOPVFRB1hQh/kTlfY52XQPkj4YUD2FCiSrmzs0eVMP7XaRJ8Pi9HUhlXXzdYSaUupl/4iqZ+czKR2ZhSJMiYklVKIpswdIeR2OuB3O7GWSOdv43YUa1oyyBWRARG5XUSO6P/tL/G4rIjcr3/d2ux2ElWLMV5/PrcTQz0eTFdIgTOsJdJMV24gxnhr9Xi3Fp5aM84MZSeoZozv2oz3aZOS5VKWI8mMpfNti+3JNVvEh7ZijDdej8cFES1rAQCSmRyyOWW6onjI58pf1wHtGs9JHfNatZL7dgDfU0rtB/A9/d/FxJVSV+hfL2xe84hqxhhvgMlw5WImAJDLKawlM+jzcyW3gRjjLVTsCCFjxp97cuuC8V0DI/OmXMqy1VWp4oPctKWBMm3AGG8wh0MQ8rqwGteuzVargRcOcjPZHOLpLK/vFrRqkHsTgE/q338SwIta1A6iRmGMN8BE2I+ppcpn5UZSGSgF9HKQ20iM8RYKepxIZxVSmfUjtbiSW1eM7xpM6NVjy52Vq1WKNX+NXi88VVip1lzqJxXFGG+CXr+7YJCrTdCYzT4I+dxY1Scvjb25rLlgXqsGuaNKqRn9+1kAoyUe5xORgyJyl4i8qNSLicgb9ccdnJubq3dbiarBGG+AibAf08sJKKXKPs74QGGBhoZijLeQ0UkqXM3lSm5d1TW+ge6K8bE+bZA7UyLzRimlbSmx0GEvVsRnNZHmZGb1eA1vgt6CgarxX7N1EwpXcjmJaV3DflMi8l0AY0Xu+ovCfyillIiU6rHuVEpNicgeAN8XkYeUUsc2P0gpdTOAmwHgwIED5Xu/RHXCGG++ibAf8XQWy7E0+oOeko9bMQa57PzUhDHevowOfzSVyf8tsBNk2Xki8nCR2+se3/rrdE2Me11ODPV4S+7JTWZySGeVtcJTbidENlabXY6luS2lvKbFeDfFtxV9fne+T7JqsW/S63Pni20a13cWFjSvYb8ppdQNpe4TkbMiMq6UmhGRcQDnSrzGlP7f4yLyQwBXAij64UHUbIzx5jOOEZpajpcd5K7G9Q8DPz8MasEYb1/rK7kbz1AEOMi14LBS6kCxOxjftZsI+zBdYk9ufkLGwn5ah0PQ61sfMORyCquJNMIc5JbDGG+xXr8LJ+a1bVZG7JqdmClcyc2nK/PUCNNala58K4DX6d+/DsBXNz9ARPpFxKt/PwTgaQAebVoLiWrDGG+AwkFuOespQfwwaCDGeAsZ+xALKyzn05XZCaoHxneNxvt8JY98qza1PhxwYzmW1l+DtRdqxBhvgsKJmVoGucYqMCcxzWvVIPd9AJ4jIkcA3KD/GyJyQET+U3/MhQAOisgDAH4A4H1KKf5hUadgjDeAUcyk0lm5Vj9IqCqM8RYqvic3A5dD4HO36qPdVhjfNRrv82NmOV60hkLEYpVZQ9jvxvKmAUM4UDqrh8pijDdBn399T671Qa4b8XQW6WwOS7EUAGCgTBYbbdSS6QCl1AKAZxe5/SCA39a//ymAS5vcNKK6YIw3xkDQA5/bgamlCiu53JPbcIzx1srvyS1IV15LpBHyuSAirWqWbTC+azcZ9iOaymI1sfU4t2r3j/cFPFWvitFGjPHm6PW7EUtpA9XVRBpOh+SPw6rE+PuIJDL5QW65rVq0Ead7iahjiIhWYblEMRPDaiIDEWv7vYg6ST5deVMRnn6ualGbGC+TebM+yLU2QO3zu7Gid/aX46n8bUTtyojP1XgaK3GtUJrZiUhjy9VKPI3FaBoepwNBkwNk4iCXiDrMZNiPqRLHUhhW42mEvC44HFzRInvKpysX7MldiqUQDrDDT+1hvE+roVCswvL6ntx6pCsz5ql9GYPc5XgaK/GtWQ3lDPZok5YL0SSW9es7M3XM4yCXiDrKZNhvKl2ZqcpkZ0ZHaSmazt+2FE1zvxa1DaNQ4HSRSclq05XDAa2ITy6n8gWouJJL7WyoxwsAmFvTBqpW4tV47nwkhcVoitd3izjIJaKOMhH2Yz6SRCKdLfmY1USalZXJ1txOB3p9LixGk/nbtJl+doKoPQyHvHA5pGi6slGIx2rhqT6/G0oBa8kM9+RSRxjpXR/kzq0lMRzymn5ufiU3ksJSLMXtKBZxkEtEHWVCXx2YLXH+IqCdk8szcsnuBnu8WIytr+QuxlLoZ+omtQmnQzDa68NMkWv1ckwrkuZyWuuGGgPalVgai9EU/G4nfG7uUaT2NVywkjsfSeZXZ80YDBoruUksRFPoD/L6bgUHuUTUUdZT4EqnLBvFHYjsrD/gzq/kxlNZJNI5Vt6ktjIRLn5WbrWrUka65kI0ibOrCYz2mh8wELVCOOCG2ymYXU1gMZqytJLrcWkZO/ORJM6uJDDa62tgS+2Hg1wi6ijGIPdMmUHuaiJtuWonUacZCHqxENEqzOaPl2A6G7WR8b7i1fC1SuDWr9FGJ//sahLn1pIYCbHTT+1NRDDc48WhmVXkFDDcY+0aPdTjxfG5KKKpbL7/Q+ZwkEtEHWW0zwuR8iu5SzEWaCD7Gwx68oNbDnKpHY2HfZhdSSCXUxtuX46l0FdFrI71aYPa2ZW4tr+RK7nUAUb7fLjr+AIAYNtAwNJzJ/v9uPsJ7blGxXIyh4NcIuooXpcTwz3ekoPcfNomO/tkc/1BDxajKSi1XmmWe3KpnUyG/UhnFeYLCqQB2nEq1cTqQMCjp37q6cpcyaUOsG+4B+msNtGzZyho6bl7C567rZ+DXCs4yCWijjPZ78dUiUHuor6iNcACDWRzg0EP0lmFtWQG59a04j5W9nsRNVr+rNxNxwgtRVMIV1E3weEQjIR8ODYXQSyVzVeuJWpn+0d78t9bTTneM7w+KC58HaqMg1wi6jgTYX/RsxcBYDHCtE3qDsbxEvNryXwFWyOdk6gdjOvxWJh5k8nmsJrIVH3c1XifDz8+MgcA2DVobVWMqBWevn8YAPC0fYOWK4pft28IAPCknf0IeHhqhBX8bRFRx5kM+3H7o2ehlIKIbLjPWMkdtFjcgajTGKtk08sJnF1JIORzsRNEbSVfDb/gGKHVRAaAVnW2GuePhXDw5BIArmxRZ7hwvBeffcOTsW/YerzuGe7B597wlA0rumQOPw2JqONMhv1IZXKYj2wtx78U5UoudQdjf9b0chwzK4n8qhlRuwgH3PC5HZgpWMmttUjaJZN9+e93WiziQ9Qq1+4dqvq5T907WMeWdA+mKxNRx5koc1buYtTYk8tBLtnbWJ8PIsDUclw/M5SDXGovIoKJTccIza1pRaiGeqrbT/uci0Yx2uvFbzxtl+XUTyLqHlzJJaKOMxFe3+d1+fbwhvsWoyk4HYJenpNLNud2OjAa8mFqOY6TizE8/5LxVjeJaIvNNRTO6YPcaotGDfV48ZM/fxYHuERUFq8QRNRxtoW1FLViFZYXYyn0B9xwOGTLfUR2s2sogLuOL2A5lsb+Ee5PpPYz3ufDTMFK7rlVbcA7UkMlcA5wiagSXiWIqOP0+l0IepxFB7lL0RT341LXuGSiD2eWtL+D80ZDLW4N0VbjYT/OrSWRyuQAAHORJDwuB/qqOEKIiMgsDnKJqOOIiJ4Ct3WQuxBNoZ/7calLHNjVn//+yh3h1jWEqITJsA9KAWf1Fdy51SSGe7xbKuMTEdUT9+QSUUcqdVbuudUELt0Wbn6DiFrgWReM4kVXTOCSyT4EvfxIp/ZjHHU1s5LA9oEAZlcTGK1yPy4RkVn8RCSijjTZ78dDUysbblNK4exqEjfUsNeLqJN4XA780yuvbHUziEoqLBQIACcXYhsyEIiIGoHpykTUkSbDfixGU4insvnbVhMZxNNZjPG8UCKitrCtPwCnQ3B8LoJUJoeZlTh2DgZb3SwisjkOcomoI+VXB4pV7eR5oUREbcHndmLfcA8enl7FmaUYcgrYORBodbOIyOY4yCWijjRpHCO0tD7IPbuqnb84ynRlIqK2cfFkLx6eWsHxuSgA7egrIqJG4iCXiDrS5n1eADCrr+SOciWXiKhtXDLRh3NrSXzz4Vk4HYILx3tb3SQisjkOcomoI432+uAQbDgrd1ZPXeYgl4iofTx9/xAA4Ev3nsGF4yEEPKx7SkSNxUEuEXUkt9OB7QOBfPoboFXtHAl54fc4W9gyIiIqtG+kJ3+O84uumGxtY4ioK3AqjYg61gVjIRyaXc3/++RCDDsHudeLiKidiAg+/vqr8YsTS3jWBSOtbg4RdQGu5BJRx7pgrBcn5qP5Y4ROLkZ5NAURURsKBzx4zkWjcDqk1U0hoi7AQS4RdawLx0PIKeDIuTWsJdI4u5rELq7kEhEREXU1pisTUce6YEyr0Pnw1Cpi+mruxRN9rWwSEREREbUYB7lE1LF2DgYwHPLiruMLiKUyAIBLJjnIJSIiIupmHOQSUccSETxt7yB+fGQey/E0tg/4MRzytrpZRERERNRC3JNLRB3tf102gYVoCnccnsNzLxprdXOIiIiIqMW4kktEHe1ZF4zg+ZeM4cRCDG94xp5WN4eIiIiIWoyDXCLqaE6H4MO//qRWN4OIiIiI2gTTlYmIiIiIiMg2OMglIiIiIiIi2+Agl4iIiIiIiGyjJYNcEXmZiDwiIjkROVDmcc8TkcdF5KiIvL2ZbSSqBWOc7I4xTnbG+Ca7Y4yT3bVqJfdhAC8GcEepB4iIE8CHADwfwEUAXiUiFzWneUQ1Y4yT3THGyc4Y32R3jHGytZZUV1ZKHQIAESn3sGsAHFVKHdcf+3kANwF4tOENJKoRY5zsjjFOdsb4JrtjjJPdtfOe3EkApwv+fUa/jcguGONkd4xxsjPGN9kdY5w6VsNWckXkuwDGitz1F0qpr9b5vd4I4I0AsGPHjnq+NFFJjHGyO8Y42dx5IvJwkdvrHt8AY5xaomkxzvimdtOwQa5S6oYaX2IKwPaCf2/Tbyv2XjcDuBkADhw4oGp8XyJTGONkd4xxsrnDSqmSBXdMMB3fAGOcWqJpMc74pnbTzunKvwCwX0R2i4gHwCsB3NriNhHVE2Oc7I4xTnbG+Ca7Y4xTx2rVEUK/KiJnADwVwDdE5Nv67RMichsAKKUyAN4M4NsADgG4RSn1SCvaS2QVY5zsjjFOdsb4JrtjjJPdtaq68pcBfLnI7dMAbiz4920Abmti04jqgjFOdscYJztjfJPdMcbJ7to5XZmIiIiIiIjIEg5yiYiIiIiIyDY4yCUiIiIiIiLb4CCXiIiIiIiIbIODXCIiIiIiIrINDnKJiIiIiIjINjjIJSIiIiIiItvgIJeIiIiIiIhsg4NcIiIiIiIisg1RSrW6DXUlInMATpa4ewjAfBObUw7bUlyntGWnUmq4mY0xMMarwrYUxxivDdtSXKe0hTFeGdtSXKe0pSUxzviuCttSXNXxbbtBbjkiclApdaDV7QDYllLYltq0U5vZluLYltq0U5vZluLYltq0U5vZluLYluq1U3vZluLs0hamKxMREREREZFtcJBLREREREREttFtg9ybW92AAmxLcWxLbdqpzWxLcWxLbdqpzWxLcWxLbdqpzWxLcWxL9dqpvWxLcbZoS1ftySUiIiIiIiJ767aVXCIiIiIiIrIxDnKJiIiIiIjINmw9yBWRd4vIlIjcr3/dWOJxzxORx0XkqIi8vUFt+XsReUxEHhSRL4tIuMTjTojIQ3p7D9a5DWV/ThHxisgX9PvvFpFd9Xz/gvfZLiI/EJFHReQREXlbkcdcLyIrBf/v/qoRbdHfq+zvXDT/ov9eHhSRqxrVFqsY41temzFevD2M8fq0hTG+/j6M8TphjG94XcZ38fYwvuvTFl7D19/H/jGulLLtF4B3A/iTCo9xAjgGYA8AD4AHAFzUgLb8MgCX/v37Aby/xONOABhqwPtX/DkB/B6Aj+jfvxLAFxr0/2UcwFX69yEAh4u05XoAX29SnJT9nQO4EcA3AQiApwC4uxntMtl2xriFn5MxXvJ+xri5tjDG19+HMV6/tjPGTf6MjO+S9zO+zbWF1/D197F9jNt6JdekawAcVUodV0qlAHwewE31fhOl1HeUUhn9n3cB2Fbv96jAzM95E4BP6t//N4Bni4jUuyFKqRml1L3692sADgGYrPf71NFNAD6lNHcBCIvIeKsbZQFjfB1jvDjGuAmM8XWM8abrhhhnfFeP8W0Cr+HruiHGu2GQ+2Z9WftjItJf5P5JAKcL/n0Gjf+f/JvQZiOKUQC+IyL3iMgb6/ieZn7O/GP0i8AKgME6tmELPQ3jSgB3F7n7qSLygIh8U0QubmAzKv3OWxEjVjDGNYzx0hjj9ccY1zHG64Ixzvguh/Fdf7yG6+wa4646Nq4lROS7AMaK3PUXAD4M4K+h/eL+GsAHoQV109uilPqq/pi/AJAB8JkSL3OdUmpKREYA3C4ijyml7mhMi1tLRHoAfAnAHyilVjfdfS+AnUqpiL5/4ysA9jeoKW39O2eMdy7GuDmM8c7FGDeHMd6ZGN/mML47l51jvOMHuUqpG8w8TkT+A8DXi9w1BWB7wb+36bfVvS0i8noAvwLg2UpPMC/yGlP6f8+JyJehpTbU4w/LzM9pPOaMiLgA9AFYqMN7byEibmh/VJ9RSv3P5vsL/9CUUreJyL+JyJBSar7ebTHxO69bjFTZPsa4OYzxEhjj9WsLY3wdY9xS+xjjlTG+S2B8168tvIavs3uM2zpdWTbmav8qgIeLPOwXAPaLyG4R8UDb5H1rA9ryPAB/BuCFSqlYiccERSRkfA9tg3yxNlfDzM95K4DX6d+/FMD3S10AaiEiAuCjAA4ppf6hxGPG9MdBRK6BFqt1/yM3+Tu/FcBrRfMUACtKqZl6t6UajPENGOPF34cxXr+2MMZ1jPH6YYznMb6Lvw/ju35t4TVc1xUxrppQMatVXwD+C8BDAB7Ufznj+u0TAG4reNyN0KqKHYOWztCIthyFlkt+v/71kc1tgVZt7QH965F6t6XYzwngPdD+2AHAB+CLelt/DmBPg34X10FLW3mw4PdxI4A3AXiT/pg367+DB6AVB7i2QW0p+jvf1BYB8CH99/YQgAOtjm3GOGOcMc4YZ4wzxu0W44xvxred45sx3twYF/2JRERERERERB3P1unKRERERERE1F04yCUiIiIiIiLb4CCXiIiIiIiIbIODXCIiIiIiIrINDnKJiIiIiIjINjjIJSIiIiIiItvgIJeIiIiIiIhsg4PcLiYiV4vIgyLiE5GgiDwiIpe0ul1E9cD4JrtjjJPdMcbJ7hjjjSNKqVa3gVpIRN4LwAfAD+CMUupvW9wkorphfJPdMcbJ7hjjZHeM8cbgILfLiYgHwC8AJABcq5TKtrhJRHXD+Ca7Y4yT3THGye4Y443BdGUaBNADIARtFonIThjfZHeMcbI7xjjZHWO8AbiS2+VE5FYAnwewG8C4UurNLW4SUd0wvsnuGONkd4xxsjvGeGO4Wt0Aah0ReS2AtFLqsyLiBPBTEXmWUur7rW4bUa0Y32R3jHGyO8Y42R1jvHG4kktERERERES2wT25REREREREZBsc5BIREREREZFtcJBLREREREREtsFBLhEREREREdkGB7lERERERERkGxzkEhERERERkW1wkEtERERERES28f8Bc14024JMI7sAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1152x288 with 5 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Shape of the eigenfunctions\n",
    "x = np.linspace(-5, 5, 1000)\n",
    "\n",
    "fig, ax = plt.subplots(nrows=1, ncols=5, figsize=[16, 4])\n",
    "ax[0].set_ylabel(\"wavefunction(x)\")\n",
    "for n in range(5):\n",
    "    eigenfunction = HO_wf(x,n)\n",
    "    ax[n].plot(x, eigenfunction)\n",
    "    ax[n].set_title(\"{}th eigenfunction\".format(n))\n",
    "    ax[n].set_ylim(-1.15, 1.15)\n",
    "    ax[n].set_xlabel(\"x\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ce70765",
   "metadata": {},
   "source": [
    "## Hartree Fock code\n",
    "\n",
    "Here we use the two-body integrals calculated in the github link below because we just want to check that the SCF Hartree Fock works correctly. In order to do that, we calculate the energy of the He-H molecule with 2 electrons and 2 single basis wave functions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "95204582",
   "metadata": {},
   "outputs": [],
   "source": [
    "class integral_master():\n",
    "    \"\"\"\n",
    "    Calculates, stores and retrieves the values of the <pr|g|qs> integrals\n",
    "    \"\"\"\n",
    "    def __init__(self):\n",
    "        self.integral_list = None\n",
    "        self.TEI          = np.genfromtxt('https://raw.githubusercontent.com/adambaskerville/adambaskerville.github.io/master/_posts/HartreeFockCode/two_elec_int.dat') # Load two electron integrals\n",
    "        self.twoe         = {eint(row[0], row[1], row[2], row[3]) : row[4] for row in self.TEI} # Put in python dictionary\n",
    "        self.ENUC = np.genfromtxt('https://raw.githubusercontent.com/adambaskerville/adambaskerville.github.io/master/_posts/HartreeFockCode/enuc.dat',dtype=float, delimiter=',') # ENUC = nuclear repulsion, \n",
    "        Sraw = np.genfromtxt('https://raw.githubusercontent.com/adambaskerville/adambaskerville.github.io/master/_posts/HartreeFockCode/s.dat',dtype=None)                    # Sraw is overlap matrix, \n",
    "        Traw = np.genfromtxt('https://raw.githubusercontent.com/adambaskerville/adambaskerville.github.io/master/_posts/HartreeFockCode/t.dat',dtype=None)                    # Traw is kinetic energy matrix,\n",
    "        Vraw = np.genfromtxt('https://raw.githubusercontent.com/adambaskerville/adambaskerville.github.io/master/_posts/HartreeFockCode/v.dat',dtype=None)                    # Vraw is potential energy matrix\n",
    "\n",
    "        dim = 2 # dim is the number of basis functions \n",
    "        self.S = np.zeros((dim, dim)) # Initialize integrals, and put them in a Numpy array\n",
    "        self.T = np.zeros((dim, dim))\n",
    "        self.V = np.zeros((dim, dim))\n",
    "\n",
    "        for i in Sraw: self.S[i[0]-1, i[1]-1] = i[2] # Put the integrals into a matrix \n",
    "        for i in Traw: self.T[i[0]-1, i[1]-1] = i[2] # Put the integrals into a matrix\n",
    "        for i in Vraw: self.V[i[0]-1, i[1]-1] = i[2] # Put the integrals into a matrix\n",
    "\n",
    "        self.S            = symmetrise(self.S) # Flip the triangular matrix in the diagonal \n",
    "        self.V            = symmetrise(self.V) # Flip the triangular matrix in the diagonal\n",
    "        self.T            = symmetrise(self.T) # Flip the triangular matrix in the diagonal\n",
    "        self.Hcore        = self.T + self.V # Form core Hamiltonian matrix as sum of one electron kinetic energy, T and potential energy, V matrices\n",
    "        return\n",
    "\n",
    "    def get_1(self, p, q):\n",
    "        \"\"\"\n",
    "        Returns the value of the h_pq integrals\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        p, q: int\n",
    "            Indeces that specify the h_pq integral\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        float\n",
    "            Value of the h_pq integral\n",
    "        \"\"\"\n",
    "\n",
    "        return self.Hcore[p,q]\n",
    "\n",
    "    def get_2(self, p, q, r, s):\n",
    "        \"\"\"\n",
    "        Returns the value of the <pr|g|qs> integrals\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        p, q, r, s: int\n",
    "            Indeces that specify the <pr|g|qs> integral\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        float\n",
    "            Value of the <pr|g|qs> integral\n",
    "        \"\"\"\n",
    "\n",
    "        p, q, r, s = p+1, q+1, r+1, s+1\n",
    "        return self.twoe.get(eint(p, q, r, s), 0) \n",
    "\n",
    "\n",
    "def eint(a,b,c,d): # Return compound index given four indices using Yoshimine sort\n",
    "    if a > b: ab = a*(a+1)/2 + b\n",
    "    else: ab = b*(b+1)/2 + a\n",
    "    if c > d: cd = c*(c+1)/2 + d\n",
    "    else: cd = d*(d+1)/2 + c\n",
    "    if ab > cd: abcd = ab*(ab+1)/2 + cd\n",
    "    else: abcd = cd*(cd+1)/2 + ab\n",
    "\n",
    "    return abcd\n",
    "\n",
    "def symmetrise(Mat): # Symmetrize a matrix given a triangular one\n",
    "    return Mat + Mat.T - np.diag(Mat.diagonal())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "06403146",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E = -3.3453739 | N(SCF) = 1\n",
      "E = -2.4186182 | N(SCF) = 2\n",
      "E = -2.4396214 | N(SCF) = 3\n",
      "E = -2.4433698 | N(SCF) = 4\n",
      "E = -2.4440162 | N(SCF) = 5\n",
      "E = -2.4441272 | N(SCF) = 6\n",
      "E = -2.4441462 | N(SCF) = 7\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import lib\n",
    "\n",
    "############################################\n",
    "\n",
    "N_electrons = 2\n",
    "N_basis = 2\n",
    "\n",
    "normalized_wf = False\n",
    "\n",
    "max_iter_SCF = 100\n",
    "eps_SCF = 1E-5\n",
    "\n",
    "############################################\n",
    "\n",
    "# Previous calculations to HF\n",
    "integrals = integral_master()\n",
    "C = np.random.rand(N_basis, N_basis)\n",
    "C = np.zeros((N_basis, N_basis))\n",
    "rho = lib.density_matrix(C, N_electrons)\n",
    "\n",
    "    # Normalization of wave functions\n",
    "if not normalized_wf:\n",
    "    # calculate S matrix\n",
    "    S = integrals.S\n",
    "    SVAL, SVEC = np.linalg.eigh(S) \n",
    "    SVAL_minhalf = (np.diag(SVAL**(-0.5))) \n",
    "    X = np.dot(SVEC, np.dot(SVAL_minhalf, np.transpose(SVEC)))\n",
    "else:\n",
    "    S = np.eye(N_basis)\n",
    "\n",
    "# Self Consistent Field\n",
    "n_iterations = 0\n",
    "rho_old = np.zeros((N_basis, N_basis))\n",
    "\n",
    "while n_iterations < max_iter_SCF:\n",
    "    n_iterations += 1\n",
    "\n",
    "    F = lib.create_F_matrix(rho, integrals)\n",
    "    \n",
    "    if normalized_wf:\n",
    "        E, C = np.linalg.eigh(F)\n",
    "    else:\n",
    "        F_prime = np.conjugate(X.transpose()) @ F @ X\n",
    "        E, C_prime = np.linalg.eigh(F_prime)\n",
    "        C = X @ C_prime\n",
    "    \n",
    "    rho = lib.density_matrix(C, N_electrons)\n",
    "    total_E = lib.total_energy(rho, F, integrals)\n",
    "    total_E = total_E + integrals.ENUC\n",
    "\n",
    "    if lib.delta_rho(rho, rho_old) < eps_SCF:\n",
    "        break\n",
    "    \n",
    "    total_E_old = total_E\n",
    "    rho_old = rho\n",
    "\n",
    "    print(\"E = {:0.7f} | N(SCF) = {}\".format(total_E, n_iterations))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2720596e",
   "metadata": {},
   "source": [
    "The energy for He-H is −2.85994 (atomic units) which is \"close\" to the one obtained here. The main difference is that we have not optimized the distance between the He and H atoms and we are not in the minimum of energy. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}