HartreeFock_JosThijssen.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "accomplished-monitoring",
   "metadata": {},
   "source": [
    "### Reference \n",
    "\n",
    "[Computational Physics](https://www.cambridge.org/core/books/computational-physics/BEE73B0139D4A9993193B57CDC62096E#fndtn-information) by Jos Thijssen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ancient-recovery",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.linalg import eigh\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "distant-moral",
   "metadata": {},
   "source": [
    "# Variational calculation for the Hydrogen atom"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "premier-dividend",
   "metadata": {},
   "source": [
    "## Basis functions\n",
    "\n",
    "Use Gaussian type orbitals (GTOs). For example, an $s$-orbital centered at $\\bf{R}$ \n",
    "\n",
    "$$ \\chi_{\\alpha,\\bf{R}}^{(s)} = \\rm{e}^{-\\alpha (\\bf{r}-\\bf{R})^2 } $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "stock-tuition",
   "metadata": {},
   "outputs": [],
   "source": [
    "def χα(α, r):\n",
    "    return np.exp(-α*r*r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "prostate-slope",
   "metadata": {},
   "outputs": [],
   "source": [
    "αα = (13.00773, 1.962079, 0.444529, 0.1219492)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "vulnerable-spanish",
   "metadata": {},
   "outputs": [],
   "source": [
    "# r grid\n",
    "rr = np.linspace(0,8,100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "hybrid-soldier",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '$χ_α$')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEHCAYAAACumTGlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA+9UlEQVR4nO3deXhU5fXA8e+ZLZkkZA8hIYGEfREQRMFqFRcERaXaasWqtUUR64K1rtRdcNf6s9pat7pW3OuGogJacQUVEESUnSyQQAjZyDKZ9/fHTCBiErLcmTsh5/M888zMnTv3PQlhzry7GGNQSimlHHYHoJRSKjJoQlBKKQVoQlBKKRWkCUEppRSgCUEppVSQy+4A2is1NdXk5OTYHYZSSnUqX3311TZjTFpTr3XahJCTk8OSJUvsDkMppToVEdnY3GvaZKSUUgrQhKCUUipIE4JSSilAE4JSSqkgTQhKKaWAMCQEEXlCRIpEZEUzr4uIPCAia0RkuYiMCnVMSimlfi4cNYQngYktvH480D94mwb8MwwxKaWU2kvI5yEYY/4nIjktnDIZeNoE1uH+XEQSRSTDGFMYingWvfMi+S/cR8+Bh5GZ2Y/klGyievfGk5uLMy4uFEUqpVSnEAkT03oCmxs9zwse+1lCEJFpBGoR9OrVq12FFS1+gwM+34nj87nU7FWIq3t3vCNHEjPmEGLHjsWTm4uItKscpZTqbCIhITT1idvkrj3GmEeARwBGjx7drp19Dp3xKFdX38GfY//NRqeXpanD2LR+FQPKYjmqLpPkZcsonzcPAE/fviScOIn4E0/Ek53dnuKUUqrTiIRRRnlA40/bLKAgVIX1iI/m2+jDmN/77/wqJoabCj7korNuZtmEvpx36Ar+fctYer3zFuk3XI8zKZHi/3uAteOPY9PU86j45BN0hzml1P4qEhLCG8A5wdFGY4Gdoeo/ABARhmYmsLAkGaa+D043B/2wkOcmPcfFB17Mm+vf4qp19+E9/RRynn2WfgvmkzbjUqp/WM3mqeex/lenUPbuPE0MSqn9TjiGnT4PfAYMFJE8EZkqItNFZHrwlLnAOmAN8Cjwp1DHNDQznh+2llMblQzDT4dvX8axq5QLRlzA9WOv5+O8j5n23jTKastwZ2aSeuGF9Js/n4zZszF1deRfdhkbz5hC1ddfhzpUpZQKm5AnBGPMFGNMhjHGbYzJMsY8box52BjzcPB1Y4y5yBjT1xgzzBgT8iVMh2TGU1dv+LGoHA4+H3zVsPQ5AE4feDr3HHkPK7av4KqPrqLeXw+Aw+Mh8den0ueN18mYdSt1BQVsPPN35F9+Ob5t20IdslJKhVwkNBmF3dDMBABWFpRBjwOg16Gw+HHw+wE4Luc4rj3kWj4p+ISHlj70k/eKy0Xib35D33nvknrRRZS//wHrJp1I6X//q81ISqlOrUsmhNzUWLxuJ98VlAUOHHwe7FgPa+fvPue0Aadxav9TefTbR5m/cf7PruGIiSHtkovJfe1VPH36UHjNteRNvxBfSUm4fgyllLJUl0wITocwOKPbnoQw+GSI7Q5fPrr7HBFh5piZDEsdxsxFM1m3c12T14rq14/ezz1L+syZVH72Gesn/4rKz78Ix4+hlFKW6pIJAQLNRt8VluH3G3B54KBz4cf3oGT97nOinFHcN+4+PE4P139yPX7jb/Ja4nCQfM7Z5Lz4Ao64ODb94Q8UP/gQxt/0+UopFYm6cEKIp6LGx6aSqsCBUWcDBr5/+yfn9YjtwVUHX8Xy4uW8uPrFFq8ZPWgQua+8TMLJJ7PtwQfJnzEDf2VliH4CpZSyVhdOCI06lgESe0FKP1j/0c/OPbHPiYzNGMv9X9/P1sqtLV7XERNDxh23kz7zWsrnL2DDlDOpzcu3PH6llLJal00I/dPjcDmElQU79xzsMw42fAL1dT85V0S4YewN+Pw+7vjyjn1eW0RIPuccsh95hLotW9jw29+ya8VKi38CpZSyVpdNCNFuJ/26x+2pIQDkHgl1lZD386kQ2fHZXDjiQj7Y9AELNy1sVRlxhx9GzpznkSgPm845h8pPP7UqfKWUslyXTQgQmKD204TwS0CabDYCOGfoOfRN6Mt9X92Hz+9rVRlRffqQ8/wc3FlZbLpgOmVz51oQuVJKWa9LJ4ShmQlsq6ihqKw6cMCbBJkHwroPmzzf7XBz6ahL2VC2gdfXvN7qctzp3en97DPEjBhB/l+uoPS1/3Y4dqWUslqXTgiDenQDYE1xxZ6DfcZB3mKoqWjyPUdlH8WItBH8Y9k/qPZVt7osZ3w82Y8+QuyhYymcOZPSl1/uSOhKKWW5Lp0QUuI8AOyobNSJnHsk+H2wsen2fhHhslGXUVRVxPPfP9+m8hxeL1n/+Aexhx1G4XXXs2POC+2OXSmlrNalE0JSTCAhlFTV7jnYayw4o5rtRwAY3WM0h/c8nMe+fYyy2rJmz2uKIzqarIceJPbII9hy003afKSUihhdOiEkxrgBKK1slBDcXug1ptl+hAYzRs2grLaMJ1c82eZyHVFRZD3wADGHjqXwr3+l7L332nwNpZSyWpdOCFEuJ7EeJzuqfjrvgD7jYOsKqChu9r2Dkgcxvvd45nw/h4rapvsbWuKIiiL7wQfxDh9O/l+uoOLjRW2+hlJKWalLJwSAxBgPpY2bjAByxwXuN/yvxfdOPWAq5XXlvPTDS+0q2xEbS/Yj/yKqXz/yLrmEXcuWtes6SillhS6fEJJi3T/tQwDIGAGuaMhveUe0oalDGZMxhme+e4ba+toWz22OMz6eXo89iistjc3TL6R248Z2XUcppTpKE0KM5+dNRk4X9BgOBd/s8/1TD5hK8a5i3lz7ZrtjcKWk0OvRR8AYNp0/Dd/27e2+llJKtZcmhKaajAAyR0LBUghuodmcsRljGZw8mCdXPrl7u8328OTkkP3wP/EVFbH5wj/h37Wr3ddSSqn20IQQ42ZHZTMJoa4Stv3Y4vtFhD8O+yMbyjawcHPr1jhqjvfAA+l57z1Uf/stBddcq/spKKXCqssnhMQYD2XVPnz1e334Zo4M3Lei2Wh8r/Fkd8vmyZVPdjiebsccQ/crr6R83jy2PfjQvt+glFIW6fIJITk2MDmtdNde/Qip/cEd26qE4HQ4OXPQmSwrXsZ327/reEx/OJeEX5/Ktn/8g51vvb3vNyillAW6fELYPTlt734EhzMw2qgVCQHg5H4n43V5mfP9nA7HJCJk3Hgj3tEHUThzJru+XdHhayql1L50+YTQsHzFz0YaQaDZaMtyqN/3UtfxnnhO7HMic9fPpbS6tMNxicdD1gMP4EpNJe/SS3XkkVIq5DQhNKxn1FzHsq8air9v1bXOGHQGNfU1/HfNfy2JzZWcTM+/P0B9SQn5l/0ZU9dE0lJKKYtoQohtpskIoOeowH0rm40GJA3goPSDmLN6ToeGoDbmHTqUjFtupmrxYoruuceSayqlVFM0IbTUZJSUC1EJrU4IAFMGTSG/Ip9F+datTZQweTJJ55xNyVNPayezUipkunxCiPE48TgdTc9FcDggcwQUtLyERWNH9zqa7t7uPL+6bXsl7Ev6lVfiPeggCm+4gZq1ay29tlJKgSYERISkWDc7mmoygmDH8grw1bTqem6Hm1P6n8Kn+Z9SWFFoXZxuNz3vuxdHdDR5M2bgr6qy7NpKKQWaEIBm1jNqkDkS/HVQ1Pr5Baf0PwWA/679rwXR7eFOT6fnPXdTu3YdhTfdhDHG0usrpbo2TQgE5iI02akMbZqx3KBnXE/GZozltR9fs6xzuUHsL35B6iUXU/bGm7ovs1LKUmFJCCIyUURWi8gaEbmmidcTRORNEVkmIitF5A/hiKtBUoyn6WGnAIm9Ax3LW9o2OezU/qdSWFnIF4VfWBDhT6VOn07sLw5l6+zbqPmx5bWWlFKqtUKeEETECTwEHA8MAaaIyJC9TrsI+M4YMwIYB9wrIp5Qx9YgKdZDaXNNRiLQfXCbmowg0LmcEJXAq2tetSDCvUJyOMi8804csbHkX345/upqy8tQSnU94aghHAKsMcasM8bUAnOAyXudY4BuIiJAHFAC7Ht6sEWSYtyU7qprvk0+fQhs/Q7a0GbvcXo4qc9JzN80nx3VOyyKdA9XWhqZd95JzY9r2Hr7HZZfXynV9YQjIfQENjd6nhc81tiDwGCgAPgWmGGM+dnazyIyTUSWiMiS4uLm9ztuq6QYD/V+Q1l1Mzmo+xCo2Qll+W267qn9T8Xn93Vo85yWxB1+GCnnn0fpCy9Q9u68kJShlOo6wpEQpIlje3/VngAsBTKBA4EHRST+Z28y5hFjzGhjzOi0tDTLAtw9Oa25foT0oYH7rW1rNuqf1J/hqcN5bc1rIRsRlHbppUQPH07hDTdQV2jdMFelVNcTjoSQB2Q3ep5FoCbQ2B+AV03AGmA9MCgMsQF7lq9odi5C98GB+zb2IwBM7jeZNaVr+L6kdeshtZW43fS8+y7w+Si46mpMvbWjmpRSXUc4EsJioL+I5AY7is8A3tjrnE3AMQAikg4MBNaFITYgsEkO0HzHsjcJ4nu2KyFMyJmA2+HmjbV7/8jW8fTuTfp111G1eDHbH38iZOUopfZvIU8IxhgfcDEwD1gFvGiMWSki00VkevC0W4FfiMi3wHzgamPMtlDH1mDPekbN1BAgUEtoY5MRQEJUAuOyxzF3/Vzq/KFbrTThlF/RbeJEih94QPdPUEq1S1jmIRhj5hpjBhhj+hpjZgePPWyMeTj4uMAYc5wxZpgx5gBjzLPhiKtBcktLYDfoPgS2rYb6tn+on9TnJEqqS/is4LP2hrhPIkLGzTfhSk2l4Kqr8O/aFbKylFL7J52pDHSLduGQFpqMINCxXF8L29u+sNzhPQ8nKSoppM1GAM6EBDJvv43a9espuufekJallNr/aEIAHA4hMcazjyaj4Fy6opVtvr7b6eb43ONZuGkhZbVl7YyydWIPPZSkc85mx3PPUfHJJyEtSym1f9GEEJQU08KKpwBpA0Gc7epHADip70nU+mt5f8P77Yyw9bpffjmevn0pvHYm9aWlIS9PKbV/0IQQlBTjYUdlC01GrihI6deukUYAQ1OGkpuQG/JmIwBHdDSZd96Jr6SELbfOCnl5Sqn9gyaEoH02GUFwCYu2NxlBoNP3pD4n8XXR1xRU7D0Nw3reA4aS+qcLKXv7bZ3FrJRqFU0IQUkx7pY7lSHQj1C6EWoq2lXG8bnHA/DO+nfa9f62Sj3/fKIPOIAtN9+Mb/v2sJSplOq8NCEEJcd6KKmqbXmJiYaO5eL2zTrO6pbFiLQRzF0/t13vbytxu8m843b8lZUU3nijbqijlGqRJoSgxBgPtT4/u+paWPohPZgQ2tlsBHBC7gn8sOMHftwRnn0Movr1I23GDCo+mE/Zm6FZZE8ptX/QhBCUFNOwnlELzUaJOeCOaXfHMsBxOcfhFGfYmo0Aks/9Pd5Ro9gyazZ1RUVhK1cp1bloQghKit3HiqcADgekDoDi1e0uJ9WbytiMscxdPzdsTTjidJJ522xMTQ1bbtS9mJVSTdOEENSq9YwA0gZ1KCEAnNDnBPIr8llWvKxD12kLT04OaX++jIqFC7XpSCnVJE0IQYnBJqN9jjRKGwjlBVDd/hnHR2cfTZQzirfXvd3ua7RH8tlnB5qOZt+mTUdKqZ/RhBCU4A0khLLqViQEgG0/tLusOE8cR2YdyXsb38PnD9tOoYjTScbsWZjqarbcdLM2HSmlfkITQlB8dDAh7NrHB3RacN+edg49bXBC7gmUVJfw5ZYvO3SdtorKzQ2MOlqwgLK3wzP8VSnVOWhCCIp2O3A7Zd81hMTe4IzqcD/C4VmHE+uO5d3173boOu2R/PtziB4xnK2zZumENaXUbpoQgkSE+Gg3Zbv2kRCcLkjt3+GEEOWM4ujso/lg0wfUtWOPhY4Qp5PM2bPxV1ayZZaudaSUCtCE0Ei8101ZdSva9FMHdLjJCGBi7kTKa8v5tODTDl+rraL69SP1oosof+ddyt57L+zlK6UijyaERuK9raghQKAfoXQT1FZ1qLxDMw4lISqBdzaEb5JaYylT/0jUkMFsueVWXSZbKaUJobH4aBc7W5UQBgIGtnds+Qm3082xvY5l4aaFVPuqO3St9hC3m8xZs6jfsYOtd94V9vKVUpFFE0IjgSajVtYQoMP9CBBoNqryVfFx/scdvlZ7RA8ZQsp557HztdeoWKQ7rCnVlWlCaCTQqdyKPoTkPoHd0yzoRxidPprk6GRbRhs1SP3ThXhyc9lyww34Kytti0MpZS9NCI3Ee12tqyG4PJDS15Iagsvh4rjex/G/vP9RVdexPon2ckRFkTF7FnWFhRT97X5bYlBK2U8TQiPx0W5qfX6qW1oCu0HaQEsSAgSajarrq/ko7yNLrtceMaNGkfS737Hjueeo+uYb2+JQStlHE0IjrV6+AiB1IJSsA19Nh8sd2X0k3b3dbW02Auj+58twZfSg8Lrr8dfuY5E/pdR+RxNCI/ENCaG1Q09NPWxf2+FyHeLguJzjWJS/iIra9m3PaQVHbCwZN99M7dq1bH/4X7bFoZSyhyaERuKjXQDsbE3H8u5F7qxpNpqQM4Fafy0LNy+05HrtFffLXxJ/8klse/RRqn9o/wJ+SqnORxNCI/FtajLqD4hl/QjD04bTI7YH722wf9Zw+rXX4oyLo/C66zH1rehPUUrtFzQhNLJnxdNWJAS3F5J6WzL0FILNRr2PY1HBIspq27/XghVcSUmk//WvVC9fzo5nn7U1FqVU+GhCaCTeG2gyatV6RhDoWC62rlllYs5EfH4fCzfZ22wEED/pBOKOPJKi+/+P2rw8u8NRSoWBJoRG2lRDAEgbANvXgN+aZpUDUg8gMzaTdzfYO9oIAqu/9rjpRkSELTfcqJvpKNUFhCUhiMhEEVktImtE5JpmzhknIktFZKWI2DIgP9rtJMrlaENCGAT1NbBjgyXliwgTcibwecHn7KzZack1O8KdkUHaXy6n8tNP2fn663aHo5QKsZAnBBFxAg8BxwNDgCkiMmSvcxKBfwAnG2OGAqeFOq7mtHo9Iwg0GUGHttPc24TcCfiMjwWbFlh2zY5ImjIF78iRFN1+B75t2+wORykVQuGoIRwCrDHGrDPG1AJzgMl7nXMm8KoxZhOAMca2HeDjo12tW88IAk1GYFnHMsCQ5CFkxWUxb8M8y67ZEeJwkDHrVvxVVWy97Ta7w1FKhVA4EkJPYHOj53nBY40NAJJE5EMR+UpEzmnqQiIyTUSWiMiS4uLikATbphpCdALE9bC0Y3l3s1Hh5+yo3mHZdTsiqm9fUi6cTtncdyhfYH+Ht1IqNMKREKSJY3v3ULqAg4BJwATgehEZ8LM3GfOIMWa0MWZ0Wlqa9ZFC67bRbCxtgGWT0xpMyJlAvaln/qb5ll63I1LPO4+o/v3Zcsst1FfYN5taKRU64UgIeUB2o+dZQEET57xrjKk0xmwD/geMCENsP5PQ2m00GzQMPbVwFM6g5EH0ju8dMc1GAOLxkDHrVnxbt1J07712h6OUCoFwJITFQH8RyRURD3AG8MZe57wO/FJEXCISA4wBVoUhtp+J97Zy17QGaQOhthzKCy2LQUQ4rvdxfLnlS0qqSyy7bkd5R4wg+ZyzKX1+DlVffWV3OEopi4U8IRhjfMDFwDwCH/IvGmNWish0EZkePGcV8C6wHPgSeMwYsyLUsTWlocmo1ePuG9Y0smgJiwYTcibgN34+2PiBpdftqLRLL8WdmRlYEbWm4yu9KqUiR1jmIRhj5hpjBhhj+hpjZgePPWyMebjROXcbY4YYYw4wxtwfjriaEu914/MbdrVmTwQIydBTgAFJA8iJz4mItY0ac8TG0uOWW6hdv55t//yn3eEopSykM5X3sme2civ7EeK6B0YbWTj0FALNRhNzJ7J462K27Yqs8f9xhx9GwuTJbH/scapXW1szUkrZRxPCXvasZ9TKfgQRy9c0ajChd6DZ6P2N71t+7Y7qfs3VOOPjKfzrdRhfGzrhlVIRSxPCXhLasklOgxAMPQXol9SPfon9bN9JrSmupCR6XH8d1StWUPL0M3aHo5SygCaEvexuMmptDQECNYTKYqiyfkTQhJwJfFP0DVsrt1p+7Y7qNnEiccccQ/EDD1C7caPd4SilOkgTwl4aNslp89BTsLxjGQJLYhsM722MrM5lCK6IesMNiNtN4fU36IqoSnVymhD20rCNZqs7lSFkQ08BchJyGJQ8KCKWxG6KO7073a+6kqovv6T0pZfsDkcp1QGaEPbSra17IgAk9AKXNyQJAQLNRsuLl1NQsfcE78iQ+JvfEDNmDEV33U3dli12h6OUaidNCHvxuBx43c629SE4HJDaz/Khpw0m5EwAiKilLBoTETJuvQXj87Hlppu16UipTkoTQhPivW1YArtB2uCQ1RCyu2VzQMoBvLP+nZBc3wqeXr1Iu2wGFR9+SNlbb9sdjlKqHTQhNCGhLUtgN0gbCGV5UF0Wkpgm5k5kVckqNpZF7mie5LPPxjtiBFtnzdLNdJTqhNqUEERkXGjCiCzx0e62jTIC6D44cB+CkUawp9koEuckNBCnk4zbZuOvqmLLrNl2h6OUaqO21hBOFpH5InKBiDwbkogiQJs2yWmQNihwXxSaRVp7xPZgVPdRvLP+nYhuo4/q25fUiy6i/N13KXsv8obKKqWa19aEEANMA7YD+daHExnatI1mg6QccEWHrGMZ4ITcE1i7cy0/7AhNLcQqKVP/SNSQwWy5+RZ8OyJj1zel1L61NSGUEdjs5jWa3gltv9CuGoLDCan9Q5oQxueMxynOiO5cBhC3m8zbbqN+50623na73eEopVppnwkhuKlNgxuBUcCTQOg++WzW5j0RGqQNhqLQ/VqSo5MZmzmWdze8G9HNRgDRgwaROn06ZW++SfmCBXaHo5RqhdbUED4XkWEAxphdxpj7jDFnG2OeCHFstknwuvEbqKxt5Z4IDUI80ggCzUb5FfksK14WsjKskjrtfKIGDqTwxhupLy21Oxyl1D60JiFcADwnIlfs/YKIRO6Qlw5oWAI70kYaARydfTRRzqiIbzaCwD7MmbffRv2OUrberk1HSkW6fSYEY8xiAnscjxKRBSLyJxH5l4h825r3d0bx7Vm+AkI+0gggzhPHEVlHMG/DPHz+yN+HIHrIEFKnTWPn629QPn++3eEopVrQmj6EWcAKYDiwBbgeiAUmGmOOC2149ohvz54IEJaRRgDH5x7P9urtLN6yOKTlWCV1+gVEDRpE4Y036agjpSJYa77hnwscHNzr+EwCiaEbcI+IxIcyOLvs2ROhjd/AwzDSCOCIrCOIc8fx9rrOsUSEeDxk3nE79aWlbJ19m93hKKWa0ZqEMMAYs3vnF2NMsTFmMvAh8EWoArNTu/sQINBsFMKRRgBRziiO7X0sH2z6gGpfdUjLskr0oEGkXjidsrfeomyeTlhTKhK1pg+hqpnj/wJOsjyiCJDoDYy0La2qbfub0waFfKQRwIl9TqSyrpIP8z4MaTlWSp02jeghQ9hy00261pFSEahDncLGmDVWBRJJukW7cAiUVrWjhhCGkUYAo9NH093bnbfXdo5mIwhOWLvzDvyVlRTedFPEz6VQqqvZL0cJdZTDISTGeNjR3hoChHSkEYDT4eSEPiewKH8RO6o7T0dtVP/+pM2YQcUH89n5+ut2h6OUakQTQjMSY9ztqyGEaaQRwKQ+k/AZH+9t6Fxt8snn/h7v6IPYOms2dYWFdoejlArShNCMpPbWEBpGGoW4hgAwMGkg/RL78fb6ztNsBIFlsjNvvx3j91MwcybG77c7JKUUmhCalRTjZkd7aggA3YeEJSGICJP6TOKbom/YXL455OVZyZOdTfo1V1P12efseOYZu8NRSqEJoVmJMZ72jTICSD8AygugqmTf53bQpNxJAJ1mTkJjiaedRtxRR1F0733U/Pij3eEo1eVpQmhGoIbQ3oQwNHC/dYV1ATUjIy6Dg3sczBtr3+h0o3ZEhIxbb8ERF0f+VVdjatv5+1ZKWUITQjMSYzxU1/mprmvjiqcQqCEAbF1pbVDNmNx3MpvLN/NN0TdhKc9KrtRUMmbdSs2qVRT//e92h6NUl6YJoRlJMYHJae2qJcR1h5jUsNQQAMb3Ho/X5eWNtW+EpTyrdTv6aBJPO43tjz1O5ef75eR3pTqFsCQEEZkoIqtFZI2IXNPCeQeLSL2I/CYccbUkOTawntGOynZ0LIsEmo3CVEOIcccwvvd45m2Yxy7frrCUabX0a6/B07s3BVdfrXsnKGWTkCcEEXECDwHHA0OAKSIypJnz7gTmhTqm1kiM6cDyFQA9hgVGGvnb0eTUDpP7TqairoIFmzrn7mSOmBgy77kH3/btFN6os5iVskM4agiHAGuMMeuMMbXAHGByE+ddArwCFIUhpn3a02TUzqGn6UPBVw0l6yyMqnmje4wmMzaz0zYbAXgPGErajEspnzePna++Znc4SnU54UgIPYHGg+Tzgsd2E5GewCnAwy1dSESmicgSEVlSXFxseaCNJcUEmoxKOjrSaMu3FkXUMoc4OKnvSXxe+DlbK7eGpcxQSJk6lZgxY9gyaxY169bbHY5SXUo4EoI0cWzv9oD7gauNMS22rxhjHjHGjDbGjE5LS7MqvibtbjKqbGdCSB0I4gxbPwLAyX1Pxm/8vLnuzbCVaTVxOMi8604cUVHk/+Uv+HUoqlJhE46EkAdkN3qeBRTsdc5oYI6IbAB+A/xDRH4Vhtia5XE5iPU4299k5I4OLGERxoTQK74XB6UfxKs/vtqp2+Dd6elk3H4bNatWUXT3PXaHo1SXEY6EsBjoLyK5IuIBzgB+0tBtjMk1xuQYY3KAl4E/GWP+G4bYWtSh2coQmI8QxoQA8Ov+v2Zz+WaWbF0S1nKt1u2oo0g652x2PPMM5Qs6Z0e5Up1NyBOCMcYHXExg9NAq4EVjzEoRmS4i00NdfkckxXZgtjIE+hF2boLqndYFtQ/H9j6Wbu5uvPLjK2ErM1S6X3EFUUMGU3jtTOoK9q5UKqWsFpZ5CMaYucaYAcaYvsaY2cFjDxtjftaJbIw51xjzcjji2pfAiqftbDKCRjOWv7MmoFbwurxM6jOJ9ze8z86a8CWiUHB4PGTddx+mro78y/+CqevAv4VSap90pnILOt5kFL41jRr79YBfU+uv7ZQL3u3Nk5NDxqxb2bV0KUX33293OErt1zQhtKBDS2ADxGdCdGLYE8Kg5EEMSRnCKz++0qk7lxvEn3ACiVPOoOTxJyhfuNDucJTab2lCaEFijIey6jrq/e38UBUJzFjeEt6EAIHO5R92/MDK7eHt1A6V9GuuIWrIYAquuZbavHy7w1Fqv6QJoQVJMW6MgZ27OlBL6DE8MNKo3mddYK1wfO7xeF1eXv4hIrpjOswRFUXW3/4Gfj/5M2bgr6mxOySl9juaEFrQoRVPG2SOBN8u2Lbaoqhap5unG8fnHs/c9XMpqy0La9mh4undm8w7bqd65Uq2zr7N7nCU2u9oQmhBUmwHF7iDQEIAKAj/XgW/Hfhbdvl28ebazjtzeW/djjmGlPPPp/TFFynV9Y6UspQmhBY0rGfUriWwGyT3gah4WxLCkJQhDE8dzpzv5+wXncsN0mZcGljv6Oab2bVy/+gjUSoSaEJogSVNRg4HZIywJSEA/HbQb9lQtoEvt3xpS/mhIC4XPe+7F2dyMnmXXIJvxw67Q1Jqv6AJoQWJwRpCaUeGnkKg2WjLCvCFf6G2CTkTSIhK4IXVL4S97FBypaSQ9cAD1G/bTv6fL8f4wttpr9T+SBNCC+KiXLgc0rEaAgQSQn0NFK+yJrA2iHJGcWq/U1mwaUGnXha7Kd5hB9Dj5pup+vxziu651+5wlOr0NCG0QERIjPFYkxDAtmaj0wacht/4efnH/WMIamOJp/yKpLPOouTJJyl97b92h6NUp6YJYR+SYtwd61QGSMoJzFi2KSFkx2fzy6xf8uLqF6mp3//G76dffRUxY8ey5YYbqPrGnt+xUvsDTQj7kGRFDUEkUEvI/9qaoNrh7CFnU1Jdwtx1c22LIVTE7Sbr/r/hysgg75JLqSsstDskpTolTQj7kBjj7ninMgQSQtF3UFfd8Wu1w5geYxiQNIBnVj2zXw1BbeBMTCT7n//AVFez+U8X4a+stDskpTodTQj7YEkNAQIJwe8L+4Y5DUSEswafxY87fuTzws9tiSHUovr2ped991KzejX5V16FqW9xR1al1F40IexDYmyghtDhb9W7O5btazaa1GcSKdEpPPPdM7bFEGpxRxxB+syZVCxYQNFdd9sdjlKdiiaEfUiK8VBb76eqtoPfNhOyICYVCpZaEld7eJwefjvot3yc/zHrStfZFkeoJZ/1O5LOPpuSp56i5D//sTscpToNTQj7sHv5Cqs6lm2sIQCcPuB0PA4Pz6zaf2sJAOnXXE3cuHFsnTVb91BQqpU0IexDYkzDAncWdCxnjYaiVWHdY3lvKd4UTu53Mm+seYPiqmLb4gg1cTrpee89RA8eTP6fL2fXsmV2h6RUxNOEsA+WrGfUoNdYwMDmxR2/Vgf8cegf8RkfT3/3tK1xhJojNpbsfz2MKy2NzdMvpHbDBrtDUiqiaULYh+TYhiYjK2oIB4PDBZs+7fi1OiA7PpuJORN5YfUL7Kyxr7YSDq7UVHo9+ggAm86fhq94/60VKdVRmhD2YU+TkQU1BE9sYOXTjZ91/FodNHXYVHb5dvGfVft/p6snJ4fsfz2Mb/t2Np0/jfqy/WPDIKWspglhHxK9bhwCxeUWLfnQ61DI/wp89i4hMSBpAOOyx/HsqmeprNv/J3F5hw8n6+8PULN2LZunX4h/1y67Q1Iq4mhC2AeX00F6fDQFpRbNMO51aGDlU5vWNWrsvGHnUVZbtt/su7wvcYcdRs+772bX0qXkzZiBqQ3/cuRKRTJNCK2QkRBN4U6LvlH2OjRwv9HefgSAEWkjGNNjDP9e8W+q6qrsDics4idOoMdNN1L5v4/Jv+JK3UdBqUY0IbRCRqKXwp0W1RBiUyB1IGyyvx8B4KKRF7G9ejtzVs+xO5SwSTr9dLpfczXl771HwTXX6hIXSgVpQmiFzIRoCkp3WbcoXO9DYdMX4Lf/g2hk95Ec3vNwnljxBOW15XaHEzYp555L2p//TNlbb1F4ww0Yv9/ukJSynSaEVshI8FLj81sz9BQCzUY1OwOT1CLAxSMvZmfNTp797lm7Qwmr1AumkfqnP7HzlVfZcuNNmhRUl6cJoRUyE6MBKCi1uB8hQpqNhqYM5dhex/LUd09RWl1qdzhhlXrJxaRMv4DSl16i8PrrNSmoLk0TQitkJHgBrOtHSOwF8T0jomO5wUUHXkRVXRVPrHzC7lDCSkRImzGD1IsuYucrr1I486/ap6C6rLAkBBGZKCKrRWSNiFzTxOu/E5HlwdunIjIiHHG1VkZCoIZg2UgjkUAtYdNnECGb1fRL6scJfU7gP6v+Q2FF19pxTERIu+RiUi+5mJ3//S8FV16FqbOoeVCpTiTkCUFEnMBDwPHAEGCKiAzZ67T1wJHGmOHArcAjoY6rLVLjonA7xbq5CAC5v4TyQij+3rprdtClIy8F4P6v77c3EJukXXQR3a/4C2Vz55J3yaX4q+3Z3U4pu4SjhnAIsMYYs84YUwvMASY3PsEY86kxZkfw6edAVhjiajWHQ0iPt3AuAkC/8YH7H9+37podlBmXyTlDzmHu+rksL15udzi2SDnvPHrcdCMVH33E5gumU1+x/8/iVqpBOBJCT2Bzo+d5wWPNmQq809QLIjJNRJaIyJLiMC9SlpngpdDKGkJCT+g+FH58z7prWmDqsKmkelO5e/Hd++Xey62RdMYZZN51J1VLlrDp97/Ht22b3SEpFRbhSAjSxLEmP2lE5CgCCeHqpl43xjxijBltjBmdlpZmYYj7lpEYTYGVNQSA/sfCps+hOnIWW4t1x3LJyEtYWryUeRvm2R2ObRJOOomshx6kZu1aNpz5O2o3brQ7JKVCLhwJIQ/IbvQ8CyjY+yQRGQ48Bkw2xmwPQ1xtkpHgZWtZNX6/hd+a+x8H/jpY/5F117TA5L6TGZQ8iPu+uq/LLGnRlG7jxtH7qSfxl5WxYcqZ7Pr2W7tDUiqkwpEQFgP9RSRXRDzAGcAbjU8QkV7Aq8DZxpgfwhBTm2UmRlNXb9hWYeEqpdljICo+ovoRAJwOJ9cccg2FlYU8vPxhu8OxlXfECHo//x8cMTFsPPscyt6LrCY+pawU8oRgjPEBFwPzgFXAi8aYlSIyXUSmB0+7AUgB/iEiS0VkSajjaquGuQgFVs1FAHC6oc+4QEKIsPb6g9IP4tT+p/L0yqdZXbLa7nBsFZWbS84Lc4geOJD8S2ew7dFHu2z/itq/hWUegjFmrjFmgDGmrzFmdvDYw8aYh4OPzzPGJBljDgzeRocjrrbYPRfBqtnKDfofB+UFUPSdtde1wOUHXU5CVAK3fHYLftO1Z/C6UlLo9dSTxJ9wPMX33kfhzL/ir7F3TwulrKYzlVspMzEENQSAfscG7iNstBFAQlQCV4y+guXblvPS6pfsDsd2juhoMu+5J7D+0WuvsfHsc6jbutXusJSyjCaEVkqKcRPlclhfQ4jPgB7D4McPrL2uRU7scyJjMsZw/9f3s6Vyi93h2E4cDtIuvYSeD/wfNWvWsP43v6Hqa/s3O1LKCpoQWklEyLRyX4TG+o2HzZ9DVYn11+4gEeHGQ2/Eb/xct+i6Lt901CD+uOPImfM8Dm8MG885h5Knn9Z+BdXpaUJog4yEEMxFABj6K/D74LvXrb+2BbK7ZXPVwVfxxZYv+M+q/9gdTsSIHjCA3JdfIu7II9l62+3kz7iM+ooKu8NSqt00IbRBhtWzlRv0GA5pg2D5i9Zf2yKn9j+VcVnj+NtXf2Nt6Vq7w4kYzvh4sh78O92vvJLy+fNZf+qvdb6C6rQ0IbRBZmI0ReXV+OotbjYRgeGnw6ZPYUdkzogVEW78xY3EeeK49uNrqa3XDeobiAgpU/9I72eexvjq2DDlzMDQVN1bQXUymhDaICPBi9/A1vIQDDccdlrg/tvIHc2T6k3lxkNvZFXJKu5afJfd4UScmFGj6PPaa3Q79liK772PTef+gbr8fLvDUqrVNCG0QUZiiOYiQGDTnF6/gOUvRNwktcaO7nU05w49lxdWv8Cba9+0O5yI40xIoOff7iNj9iyqV6xg3cmTKX35Ze1wVp2CJoQ2yAzFbOXGhp8O236AwmWhub5FZoyawej00dzy2S38sCMiVxqxlYiQ+Otfk/vGG0QfcACF113P5gsuoK7gZ0t4KRVRNCG0QUhrCBAYbeT0RHTnMoDL4eLuI++mm6cbf174Z3bW7LQ7pIjkyepJr38/Qfpf/0rV4iWsO/EkSp55VrfoVBFLE0IbxEe7SYpxs644RJumeJMCS1l8+xLU+0JThkVSvancO+5eCisLuWzhZdrJ3AxxOEg++yz6vPkm3oMOYuvs2Ww480x2rVxpd2hK/YwmhDYanpXIsrzS0BVw4JlQWQTfR377/MjuI7n1sFtZsnUJ132ik9Za4snqSfYj/yLz7rupy8tnw29Oo/Dmm6kvLbU7NKV204TQRgdmJ/LD1nIqa0L0DX7AREjuC4vuj+jO5QaT+kxixqgZvLP+Hf7+zd/tDieiiQgJJ51I33fmknTWWZS+8CJrJx5PyXPPYerq7A5PKU0IbXVgdiJ+A9/mh6jd3OGEX1wChUth/f9CU4bFph4wldMGnMZj3z7GUyufsjuciOeMj6fHX2eS++orRA0YwNZbZ7Fu8q8oX7hQRyMpW2lCaKMR2YkALN1cGsJCpkBsd/jk/tCVYSERYeaYmRzX+zjuWXIPT6982u6QOoXoQYPo9dSTZD30INTXk3fhn9h49tlULYm47UBUF6EJoY2SYz30Tolh6abS0BXijoax02HtgogfgtrA5XBxxxF3ML73eO5ecjfPfPeM3SF1CiJCt2OOoc9bb5J+w/XUbdzExrPOZtO0aexavtzu8FQXowmhHQ7MTgxtDQFg9FTwdINPHghtORZyO9zcecSdjO89nrsW38Ujyx/RJpBWEreb5DPPpO978+h+5RVUL1vOhtN/y6bzzqfq66/tDk91EZoQ2uHA7ES2lFWzJVQT1AC8iTD6XFj5KmxbE7pyLNaQFE7scyJ//+bv3PL5Lfj8kT2ENpI4vF5Spk6l7/z5dL/iL1R/9x0bz/wdG353FuULFuj6SCqkNCG0w4G7+xF2hLagQy8BTxzMvaJTjDhq4Ha4ue3w2zhv2Hm8/MPLzFg4g6q6KrvD6lSccbGknHce/T54n/SZ1+IrLCTvTxexbtKJlPznP9RXhGgujOrSNCG0w+CMeNxOYenmEM/Q7ZYOx9wA6xbCildCW5bFRIQZo2Zw/djrWZS/iClvT2Fd6Tq7w+p0HDExJJ9zDn3fm0fmvffgiIlh6y23smbcOLbMvo2atboUubKOJoR2iHY7GZIRH/oaAsDoP0LmSJg3E3aVhr48i50+8HT+Nf5flNaUcsbbZ/D2urftDqlTEpeLhEmTyHn5JXLmPE/cUUexY84c1k06kQ1n/o7SV1/DX6W1MNUxmhDa6cDsRL7N20m9P8RNOQ4nnPg3qCyGBbNCW1aIjM0Yy4snvsjg5MFc8/E1XLfoOl3/qJ1EBO+BB9Lz7rvo/+FCul95BfUlJRTOnMkPh/+SgquvpuKTT3S9JNUumhDa6cBeiVTW1vNjUXnoC8scCYdMg8WPwYZFoS8vBNJj03lswmOcP+x83lr3Fqe8fgofbv7Q7rA6NVdKCilTp9Lnnbn0fvYZEiadQPmChWyeeh4/jhvHlltupWrxYu2IVq0mnXVY4OjRo80SGyfwrN9WyVH3fMgdpw7jjEN6hb7A6jJ47JhATeG8+ZDSN/RlhsjK7Su5btF1rCldw/je47n8oMvJ6pZld1j7BX9NDRULF1L2zrtUfPQRproaZ2oq3Y4aR9zRRxN76KE4oqPtDlPZSES+MsaMbvI1TQjtY4xh5K3vc1i/VB46c1R4Ci1ZB48eAzHJcN4HgdVRO6m6+joeX/E4T6x4Ap/fx1lDzuL8YefTzdPN7tD2G/7KSio++ojyDz6g4qP/4a+sRKKjiRlzCHFHHEHc4Yfj7tULEbE7VBVGmhBCZPbb3/H4ovV8cPmR9EmLC0+hGz+Dp0+G7DFw1qvg8oSn3BDZWrmVB755gDfWvkE3dzemDJ7CWYPPIim68ya7SOSvraXqiy+p+OgjKj7+H3UbNwHgyswgduyhxI4dQ8zBB+POyLA5UhVqmhBCpLi8hl/etYBJwzK59/QR4St42Qvw2jTIPRJOezJQY+jkVm1fxaPfPsr7G9/H6/JySr9TOH3g6fRN7LxNY5GsdsMGKj79lKrPPqfyiy/wl5UB4M7MxDv6ILwHHoh3xAiiBw5EXC6bo1VW0oQQQre+9R1PfrqBBX85kt4pseEreOnz8OalkJAFZ74Iqf3DV3YIrS1dy+PfPs47G97B5/dxUPpB/Lr/rzkq+yjiPGGqhXUxpr6emtWrqVryFVVLllD19dfUb9sGgERHEz1oENFDhwZugwcR1bcv4uncNdOuTBNCCBWVVfPLuxZy8ohM7j4tjLUEgE1fwJwzob4OJt0Dw06D/aQ9ePuu7by+9nVeWv0SeRV5eBwejsg6gvG9x/OLzF+QGJ1od4j7LWMMdfkF7Fq2lOrly9m1ciXV363CNMxzcLmI6tOHqH79iOrfD0+/fkTl5uLp1UsTRSegCSHEbnpjJc98vpEPrxhHdnJMeAsv3QQv/h4Kvoasg2HiHZDV5L91p+Q3fpYXL+ed9e/w3sb32LZrGw5xcEDqARyWeRij0kcxPHU4Me4w/967GFNfT+3GjdR8/z3V36+mZvVqatasoS4/f89JDgfurCw8vXrh6ZWNu1cvPFlZuHv2xJ2VhbObDhiIBJoQQmxrsJZwSE4yf58ykqTYMH9L8vth2fMw/2ao2Ar9joWRZ8HAE8AVFd5YQqjeX8+K7StYlL+IRXmLWLl9JQaDS1wMTB7IkJQhDE4ZzKCkQfRJ7EOsO4xNeF2Uv7KSmnXrqN2wgdr166lZv566jZuo3bwZf/lP5+g44uJw9UjH3SMDV3p3XN27405Px5Waiis1FWdqKq6UFBxer00/Tddge0IQkYnA/wFO4DFjzB17vS7B108AqoBzjTEtrvkbSQkB4PkvN3HD6ytIiY3ib789kEP7poQ/iJpy+Owf8PXTUJYXGJbafwLkHA45h0FS7n7TpARQXlvO0qKlfFP0DcuKl7GqZBXltXs+hNJj0slJyCErLousbln0jOtJ95juu29Rzv0nWUYaYwz1paXU5RdQl5dHXX4+dVu24NtSSF3hFnxFRfi2bQt8mdmLeL24kpJwJiXhTEwM3Cck4EyIxxEfjzM+AUe3OJxxcTjiuuGIiw08jo1FvF4dRrsPtiYEEXECPwDjgTxgMTDFGPNdo3NOAC4hkBDGAP9njBnT0nUjLSEArMjfyaXPf8P67ZX86sCeHJKbzIisRPqnx+F2hnFSuL8e1n0YqDWsXQhVgQ5CohMhdUDglpQDcd2hWw+ISYWobhAdD55YcHnB2flGlhhjyK/IZ3XJataXrWdd6To2lG0gvyKfkuqSn50f544jOTqZpOgkEqISiPfEE++JJ84TR5w7jlh3LF6Xd/ct2hVNlDOKKGcUHqcHt8ONx+nBJS5cjj03pzj1Q6kVTH09vm3b8W0rpn7bNnzbtuHbXkJ9SQn1O0rw7dhBfelO6nfsoH7nzkCNY1+fVyI4YmKQGC8ObwwOrxeH14t4o3FEe5HoKBxR0UhUFI7oKMTjQTwN9x7E4w7cu917bi434nYhrsANpyvw3OmE4DFxOAKPHQ5wOn/63OHYfRyHI/C34XSCiC1/J3YnhEOBm4wxE4LPrwUwxtze6Jx/AR8aY54PPl8NjDPGFDZ33UhMCACVNT5uf2cVby0vpLRqz8bpUS4HcVEuot1OnA7B6RBEoPGfQ+M/Dsv+TIyhl38zw+tX0Kd+A9n+fHr5N5NsSlt8Wz0OanFTjwufOKnHhR/BjwM/gsGBCQSNIXAD2PPXtOdYq0Nt44/WFrsEilyw3QklzsD9zoabA8odUOGESgdUCZgO/kd1GIMDcJhAtdgBiAncO3afE7wPPhf2/LuLafS40XWF5v82mj3eOVuFf0aMIboGvDUQHbx5ayCq1hBVC9G1EFUHnob7OvDUGTx14PYFbp46cPnAXR947qoP3OziBxAwjW/suUf2ehx8bcPBGZz2yIJ2ldlSQgjH18CewOZGz/MI1AL2dU5P4CcJQUSmAdMAevUKw3IR7RAb5WLWr4Zx6+QD2Li9iqWbS9lcUkVFrY+Kah+76urx+w31BvyNk/FPHlr9P3gIqxnC6kZHXP5autXvIKG+hNj6nXj9lXjrK4kyVbj9dbhNLS5Ti5N6nMaHw/hwYBATSAk0pAFjdge/5wNsz7HWkpCmg4BuwVuOAXzBWxP8GGox7HL4qRFDDYH7uoYbgXsfBp8YfAL1mMBNAu+vB/xi8EPwZnb/x/YH/4Ub/5b8AjQ63mDvxz/9LZkmz2v6DIvZVQGKCtz8BNqWLVnf1Ric9QRvjR77A48dfnDUG5z+4GM/OPym0ePAczGBx2ICyV78gS8GEnzckOQDxwPlimmU/E3D8eB5Df94u88xu19z9OxhxU/+M+FICE396ez9d9qaczDGPAI8AoEaQsdDCx0RISc1lpxU7dhUSnUO4WjYzgOyGz3PAgracY5SSqkQCkdCWAz0F5FcEfEAZwBv7HXOG8A5EjAW2NlS/4FSSinrhbzJyBjjE5GLgXkE+teeMMasFJHpwdcfBuYSGGG0hkCz4B9CHZdSSqmfCsvYQmPMXAIf+o2PPdzosQEuCkcsSimlmqY7pimllAI0ISillArShKCUUgrQhKCUUiqo0652KiLFwMZ2vj0V2GZhOFaJ1LggcmPTuNpG42qb/TGu3saYtKZe6LQJoSNEZElza3nYKVLjgsiNTeNqG42rbbpaXNpkpJRSCtCEoJRSKqirJoRH7A6gGZEaF0RubBpX22hcbdOl4uqSfQhKKaV+rqvWEJRSSu1FE4JSSimgCyYEEZkoIqtFZI2IXGN3PAAi8oSIFInICrtjaUxEskVkoYisEpGVIjLD7pgARCRaRL4UkWXBuG62O6bGRMQpIt+IyFt2x9JARDaIyLcislREImbvWRFJFJGXReT74N/ZoREQ08Dg76nhViYil9kdF4CI/Dn4N79CRJ4XkWhLr9+V+hBExAn8AIwnsCnPYmCKMeY7m+M6AqgAnjbGHGBnLI2JSAaQYYz5WkS6AV8Bv4qA35cAscaYChFxA4uAGcaYz+2Mq4GIXA6MBuKNMSfaHQ8EEgIw2hgTUZOsROQp4GNjzGPB/VJijNnHht9hFPzMyAfGGGPaOxHWqlh6EvhbH2KM2SUiLwJzjTFPWlVGV6shHAKsMcasM8bUAnOAyTbHhDHmf0CJ3XHszRhTaIz5Ovi4HFhFYK9rW5mAiuBTd/AWEd9sRCQLmAQ8ZncskU5E4oEjgMcBjDG1kZQMgo4B1tqdDBpxAV4RcQExWLyzZFdLCD2BzY2e5xEBH3CdgYjkACOBL2wOBdjdLLMUKALeN8ZERFzA/cBVBPaBjyQGeE9EvhKRaXYHE9QHKAb+HWxie0xEIm0T8jOA5+0OAsAYkw/cA2wCCgnsLPmelWV0tYQgTRyLiG+WkUxE4oBXgMuMMWV2xwNgjKk3xhxIYP/tQ0TE9qY2ETkRKDLGfGV3LE04zBgzCjgeuCjYTGk3FzAK+KcxZiRQCUREvx5AsAnrZOAlu2MBEJEkAi0auUAmECsiZ1lZRldLCHlAdqPnWVhc5drfBNvoXwGeM8a8anc8ews2MXwITLQ3EgAOA04OttfPAY4WkWftDSnAGFMQvC8CXiPQfGq3PCCvUe3uZQIJIlIcD3xtjNlqdyBBxwLrjTHFxpg64FXgF1YW0NUSwmKgv4jkBrP/GcAbNscUsYKdt48Dq4wx99kdTwMRSRORxOBjL4H/KN/bGhRgjLnWGJNljMkh8Le1wBhj6Te49hCR2OCgAIJNMscBto9oM8ZsATaLyMDgoWMAWwcs7GUKEdJcFLQJGCsiMcH/m8cQ6NezTFj2VI4UxhifiFwMzAOcwBPGmJU2h4WIPA+MA1JFJA+40RjzuL1RAYFvvGcD3wbb6wFmBvfItlMG8FRwBIgDeNEYEzFDPCNQOvBa4DMEF/AfY8y79oa02yXAc8EvaOuAP9gcDwAiEkNgNOIFdsfSwBjzhYi8DHwN+IBvsHgJiy417FQppVTzulqTkVJKqWZoQlBKKQVoQlBKKRWkCUEppRSgCUEppVSQJgSllFKAJgSllFJBmhCUCpHgxDmlOo0uNVNZqVATkZcIrKg7EpgPzLI3IqVaTxOCUtYaRmDtp6PsDkSpttKlK5SySHA7w01ApjHGZ3c8SrWV9iEoZZ2hwBeaDFRnpQlBKesMA5bbHYRS7aUJQSnraEJQnZr2ISillAK0hqCUUipIE4JSSilAE4JSSqkgTQhKKaUATQhKKaWCNCEopZQCNCEopZQK+n+0BQwNqW5y5QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# store wavefunctions in a matrix\n",
    "χrα = np.array( [ [ χα(α, r) for α in αα] for r in rr ] )\n",
    "\n",
    "plt.plot(rr, χrα)\n",
    "plt.xlabel(r\"$r$\")\n",
    "plt.ylabel(r\"$χ_α$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "attractive-trading",
   "metadata": {},
   "source": [
    "## Variational wavefunction\n",
    "\n",
    "$$\n",
    "\\Psi(\\mathbf{r}) = \\sum_\\alpha C_\\alpha \\chi_\\alpha(\\mathbf{r})\n",
    "$$\n",
    "\n",
    "with variational parameters $C_\\alpha$. \n",
    "\n",
    "Writing the Hamiltonian in this basis one obtains the generalized eigenvalue equation\n",
    "$$\n",
    "\\mathbf{H} \\mathbf{C} = E \\, \\mathbf{S} \\mathbf{C}\n",
    "$$\n",
    "\n",
    "where the Hamiltonian and matrix elements are given by\n",
    "\n",
    "$$\n",
    "H_{pq} = \\langle{\\chi_p}|H|\\chi_q\\rangle, \\qquad \n",
    "S_{pq} = \\langle{\\chi_p}|\\chi_q\\rangle\n",
    "$$\n",
    "\n",
    "Note that we don't need to normalize the wavefunctions."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "widespread-lingerie",
   "metadata": {},
   "source": [
    "## Overlap integrals and matrix elements\n",
    "\n",
    "![3.29](pngs/3_29.png)\n",
    "\n",
    "Atomic units: Energy unit is Hartree."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "id": "incomplete-invention",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_Spq(αp,αq): # overlap\n",
    "    return ( np.pi / (αp+αq) )**1.5\n",
    "\n",
    "def calc_Tpq(αp,αq): # kinetic\n",
    "    return  3 * np.pi**1.5 * αp * αq / (αp+αq)**2.5\n",
    "\n",
    "def calc_Apq(αp,αq): # Coulomb\n",
    "    return -2 * np.pi / (αp+αq)    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "id": "confused-reset",
   "metadata": {},
   "outputs": [],
   "source": [
    "# form Hamiltonian and overlap matrices\n",
    "HM = np.array( [ [calc_Tpq(αp,αq) + calc_Apq(αp,αq) for αp in αα] for αq in αα ] )\n",
    "SM = np.array( [ [calc_Spq(αp,αq)                   for αp in αα] for αq in αα ] )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "id": "purple-proposal",
   "metadata": {},
   "outputs": [],
   "source": [
    "EE, UU = eigh(HM, SM)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "id": "derived-library",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "$E_0=$ -0.499278405667485\n"
     ]
    }
   ],
   "source": [
    "# print Hydrogen atom ground state energy\n",
    "print(r\"$E_0=$\"+f\" {EE[0]}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "id": "unable-plumbing",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1533d7fd7d00>"
      ]
     },
     "execution_count": 105,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEHCAYAAAC9TnFRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAArwElEQVR4nO3dd3iV9f3/8ec7yclOCCEJKwlLlD1DRLSA4kAcOLCC1FqrIqLtr8Pd1jpaa6fWakv5Or5aByJS9Gu1VgVEQZSgqOw9QhghQPbO5/dHIkUIcAjJuU+S1+O6uMw59+fc53U+l+TFPY855xARETmeEK8DiIhI86DCEBERv6gwRETELyoMERHxiwpDRET8EuZ1gIZKSkpyXbt29TqGiEizsmzZsr3OueSGvLbZFkbXrl3JysryOoaISLNiZlsb+lrtkhIREb+oMERExC8qDBER8UuzPYYhIs1bZWUl2dnZlJWVeR2lRYqMjCQ1NRWfz9do61RhiIgnsrOziYuLo2vXrpiZ13FaFOcceXl5ZGdn061bt0Zbr3ZJiYgnysrKaNeuncqiCZgZ7dq1a/StNxWGiHhGZdF0mmJum29hFO7yOoGISKvSvAsjb6PXKUSklRo3bhwHDhw45piHH374G49HjBjR6Dm+973vMXv27EZfb32ab2EY8MnfvU4hIq2Mc46amhreeustEhISjjn28MJYvHhxEyZres23MKLawucvQOkBr5OISDN011138de//vXg4/vvv58HHniAMWPGMGTIEPr378/rr78OwJYtW+jduzfTpk1jyJAhbN++na5du7J3714ALrvsMoYOHUrfvn2ZMWMGAHfffTelpaUMGjSIyZMnAxAbGwvUls4dd9xBv3796N+/P6+88goACxYsYPTo0UyYMIFevXoxefJkvv5W1AcffJBhw4bRr18/pkyZghffltp8T6uNSYHKbPjseTjzh16nEZGT8MD/rWRVTkGjrrNPp3h+eUnfoy6fOHEiP/rRj5g2bRoAs2bN4t///jc//vGPiY+PZ+/evQwfPpxLL70UgLVr1/Lss89+o2S+9swzz5CYmEhpaSnDhg3jyiuv5JFHHuGJJ55g+fLlR4yfM2cOy5cv54svvmDv3r0MGzaMkSNHAvD555+zcuVKOnXqxJlnnsmiRYs466yzuO2227jvvvsAuPbaa3nzzTe55JJLTnaaTkjz3cLwRUHXb8GnM6C6yus0ItLMDB48mD179pCTk8MXX3xB27Zt6dixI/feey8DBgzg3HPPZceOHezevRuALl26MHz48HrX9fjjjzNw4ECGDx/O9u3bWb9+/THf+6OPPmLSpEmEhobSvn17Ro0axdKlSwHIzMwkNTWVkJAQBg0axJYtWwCYP38+p59+Ov3792fevHmsXLmy8SbDT813CwNg+DSYOQlWvwH9rvA6jYg00LG2BJrShAkTmD17Nrt27WLixIm8+OKL5ObmsmzZMnw+H127dj14LUNMTEy961iwYAHvvfceH3/8MdHR0YwePfq41z8ca3dSRETEwZ9DQ0OpqqqirKyMadOmkZWVRVpaGvfff78nV8g33y0MgFMvgLbdYMmRm4giIsczceJEZs6cyezZs5kwYQL5+fmkpKTg8/mYP38+W7ce/07g+fn5tG3blujoaNasWcOSJUsOLvP5fFRWVh7xmpEjR/LKK69QXV1Nbm4uCxcuJDMz86jv8XU5JCUlUVRUFLCzog7XvAsjJBSG3wLZS2H7Uq/TiEgz07dvXwoLC+ncuTMdO3Zk8uTJZGVlkZGRwYsvvkivXr2Ou46xY8dSVVXFgAED+MUvfvGN3VZTpkxhwIABBw96f+3yyy9nwIABDBw4kHPOOYff/e53dOjQ4ajvkZCQwE033UT//v257LLLGDZsWMM/9EkwL460N4aMjAyXlZUF5UXwpz5wyhi46lmvY4mIn1avXk3v3r29jtGi1TfHZrbMOZfRkPU17y0MgIhYGPpdWPU6HNjudRoRkRar+RcGQObNtf/9dIa3OUREWrCWURgJadD7Elj2XO0uKhERaXQtozAAzrgVyvNh+UteJxERaZFaTmGkZULnDPjkb1BT43UaEZEWp+UUBsAZ02DfJlj/jtdJRERanGZbGPuKK6isPmxLovd4iE+FxX/xJpSICDB37lxWrVrldYxG12wLY8eBUi54dCH/XrHrv5fZh4bVHsvYugi2LTn2CkREmogKI8h0aRdNSIgx9YVlTJj+MVlb9tUuGHodRLeDD//obUARCXovvPACmZmZDBo0iJtvvplPPvmEAQMGUFZWRnFxMX379mXFihUUFRXVe9tzgOeff/7gVdvXXnstixcv5o033uCOO+5g0KBBbNzYcr7orVlf6b3kk0+ZvSybP727jj2F5ZzXpz13jT2NU9ZMh3m/gps/hI4DvI4qIvX4xlXIb98Nu75q3Dfo0B8ufOSY73/nnXcyZ84cfD4f06ZNY/jw4axbt46ysjJKS0tJTU3lnnvuoaqqipKSkm/c9nz9+vWsWrWKK664gkWLFpGUlMS+fftITEzke9/7HhdffDETJkxo3M90ghr7Su9mfbfasNAQJmamc+mgTjzz0Wamf7CJ8x9dyD3nnMtNEY/DR3+Cq/7X65giEoTef/99li1bdvC+TKWlpaSkpHDfffcxbNgwIiMjefzxx4Hau8vee++9LFy4kJCQkIO3PZ83bx4TJkwgKSkJgMTERM8+TyAEpDDMbCzwZyAUeMo598hhy0cDrwOb656a45x70N/1R4eHcds5PZmUmc59r6/k1+/ncOHQyaSunA5nr4ekno3zQUSkaRxjS6CpOOe47rrr+M1vfvON53ft2kVRURGVlZWUlZURExNz1NueO+cws4Bn90qTH8Mws1DgSeBCoA8wycz61DP0Q+fcoLo/fpfFodrFRvCHqwbSPSmGKeuH4cIi4KPHGh5eRFqsMWPGMHv2bPbs2QPAvn372Lp1K1OmTOGhhx5i8uTJ3HXXXQBHve35mDFjmDVrFnl5eQfXARAXF0dhYaEHn6ppBeKgdyawwTm3yTlXAcwExjfVm0WFh/Lo1YNYVxTFwrhx8OVM3ZRQRI7Qp08ffvWrX3H++eczYMAAzjvvPJ577jnCwsK45ppruPvuu1m6dCnz5s076m3P+/bty89+9jNGjRrFwIED+clPfgLUfs/G73//ewYPHqyD3if0BmYTgLHOuRvrHl8LnO6cu+2QMaOB14BsIAe43Tl3xPcPmtkUYApAenr60GN9ucnj76/n5Xc/ZlHUjwkZ9n0Y9/vG+1AictJ0e/Om1xxvb17fDr7DW+ozoItzbiDwF2BufStyzs1wzmU45zKSk5OP+abTRvegfVoP5rpv4ZY9D0V7GhBdRES+FojCyAbSDnmcSu1WxEHOuQLnXFHdz28BPjNLOpk3DQsN4dGrB/E/1ZfiqitwH+trXEVETkYgCmMp0NPMuplZODAReOPQAWbWwepONTCzzLpceSf7xt2SYrj24nP5V3UmlUtmQOn+k12liDSi5nodWHPQFHPb5IXhnKsCbgPeAVYDs5xzK81sqplNrRs2AVhhZl8AjwMTXSN92kmZaWSlXU94dTF75z/ZGKsUkUYQGRlJXl6eSqMJOOfIy8sjMjKyUdfbrK/0zsrK8mtsbmE5a/44lgG2gag7VhEeHdfE6UTkeCorK8nOzqasrMzrKC1SZGQkqamp+Hy+bzzfaq/09ldyXAQbz7mTNvMm8tHsP3DWdx/wOpJIq+fz+ejWrZvXMeQENNubD56o4SMvZFXUEPpsfJp9+/d5HUdEpNlpNYUBEHvh/SRaIV+8GvjbEIiINHetqjDSB4xiZdwIhu74Bzm7dnodR0SkWWlVhQGQdOlDxFsJq2b/2usoIiLNSqsrjPY9M1jZ9lzOyJ3Fps2bj/8CEREBWmFhAHS+/EEiqWDj3Ie8jiIi0my0ysJISO/L2g4XM/LAG3y1erXXcUREmoVWWRgAXa94kBBz7Py/h3SlqYiIH1ptYUS3787G1Cs4u/jffPLZMq/jiIgEvVZbGADdr7yfGgul8J2HqanRVoaIyLG06sIIb9uZbT2u4Zzyecxf9KHXcUREglqrLgyAHpf9nHKLxBY8QkVVjddxRESCVqsvjJC4ZPb0uZ5zqhfx7vvveB1HRCRotfrCAOhyyV0UWDztl/yK0vIqr+OIiAQlFQZgUQnsz/wpGW4FC/7vea/jiIgEJRVGnS7n30pOWBq9V/yO/MJir+OIiAQdFcbXQn1UnfsQXdnJ0tl/8DqNiEjQUWEcIv30y1gTPZSMLTPIzd3ldRwRkaCiwjiUGXGXPkI8xayb9Uuv04iIBBUVxmE698pkWbuLGLbnVXI2rvQ6johI0FBh1KPLlQ9TSRi5c+/2OoqISNBQYdQjpXMXlqVex8DChWz9/F2v44iIBAUVxlH0v+pn7HaJ1Lz9M6jRLUNERFQYR9E2IYEve/0/ulWsZfOC//U6joiI51QYxzDi8mmsojtxH/0aV17kdRwREU8FpDDMbKyZrTWzDWZ21CPJZjbMzKrNbEIgch1PTGQ4WzPvI6lmL5vn3O91HBERTzV5YZhZKPAkcCHQB5hkZn2OMu63QFDdMva8C8bzju9c0tc+Q9kOnWYrIq1XILYwMoENzrlNzrkKYCYwvp5xPwBeA/YEIJPfwkJDSLzsNxS5SHJfuQ30/d8i0koFojA6A9sPeZxd99xBZtYZuByYfqwVmdkUM8sys6zc3NxGD3o0w/qeytvtp5BW8Bl5H/8jYO8rIhJMAlEYVs9zh/8z/THgLudc9bFW5Jyb4ZzLcM5lJCcnN1Y+v4yaeDtfuFPwvX8flB4I6HuLiASDQBRGNpB2yONUIOewMRnATDPbAkwA/mpmlwUgm986JcayZsgviak6QM4/f+51HBGRgAtEYSwFeppZNzMLByYCbxw6wDnXzTnX1TnXFZgNTHPOzQ1AthMyftxFzPVdSId1L1CZ/ZnXcUREAqrJC8M5VwXcRu3ZT6uBWc65lWY21cymNvX7N6ZIXyjtLn6QPBfP/lm3Qc0x96CJiLQoAbkOwzn3lnPuVOdcD+fcr+uem+6cO+Igt3Pue8652YHI1RCjBp7CnKRbSClYScHip72OIyISMLrS+wSZGeddfRtLavoQNu9BKArc2VoiIl5SYTRA95Q4vhp0H77qEvLm3ul1HBGRgFBhNNCki87jhdDxtNswh6q1QXVxuohIk1BhNFBsRBidx9/P+prOlL52G5Tlex1JRKRJqTBOwvkDu/Bmt58RXZ7Lvrl3eR1HRKRJqTBO0vev/jYvhV5C4pqXqVj3vtdxRESajArjJLWJ9pF+xa/ZWNOR0tm3Qnmh15FERJqECqMRjOqXzr97/Jy48l3k/vMer+OIiDQJFUYjufbqq3k1bBzJa/5BxYYPvI4jItLoVBiNJD7SR9qER9hak0Lxq7dARbHXkUREGpUKoxGN6J3Oez1/TtvyHezSHW1FpIVRYTSyid+ezJywcaSsfpayjR95HUdEpNGoMBpZTEQYaVf9lmyXRMkrN0JZgdeRREQahQqjCQw7LZ2F/R4mvnw321+Y5nUcEZFGocJoIhOvmMBrsdeQlv1/7P7oea/jiIicNBVGEwkLDWHkjY/wOacR+96dlOdu8jqSiMhJUWE0oY5t4yi7ZDrVDnY9ey1UV3kdSUSkwVQYTeyMoUNY0PMeupSsYM2sX3gdR0SkwVQYATBu0g9YEDmGnmv+xo4v5nkdR0SkQVQYARAWGkKvG/5OjqUQNvdmygr3eR1JROSEqTACpENyMnvOe4J2NXtZ+9SN4JzXkURETogKI4CGnnk+H6fdxMD891k293Gv44iInBAVRoCdcd2v+SpiMP2WP8T65R96HUdExG8qjAAL8/lIvfElDlg8MXO/z949u7yOJCLiFxWGB9omd6L4smdIcnlse/o7VFTq+gwRCX4qDI90HzSaNYPuZUj5UhY+fbfXcUREjisghWFmY81srZltMLMjfjua2Xgz+9LMlptZlpmdFYhcXhtw2U9Z0W4s5+x8igX/mul1HBGRY2rywjCzUOBJ4EKgDzDJzPocNux9YKBzbhDwfeCpps4VFMzofdPT7AjvwoBPf8qXK1d4nUhE5Kj8Lgwz69jA98gENjjnNjnnKoCZwPhDBzjnipw7eGFCDNBqLlIIjYylzXUzCbcaQmdfx+59+V5HEhGp14lsYfwawMwmm9kiM7vIz9d1BrYf8ji77rlvMLPLzWwN8C9qtzKOYGZT6nZZZeXm5p5A9OAWn9qb/Asep6/bwGczplJWWe11JBGRI5xIYRyo++/5wFnA5X6+zup57ogtCOfcP51zvYDLgIfqW5FzboZzLsM5l5GcnOzn2zcPnc+4ik2n3siFZW/x+lMP4XQluIgEmRMpjDAz+zmwrW73UbGfr8sG0g55nArkHG2wc24h0MPMkk4gW4vQfeLv2JJ4Flfu+jNzXnvB6zgiIt9wIoXxU2ARdbumgDA/X7cU6Glm3cwsHJgIvHHoADM7xcys7uchQDiQdwLZWoaQULpMeYncyC6c+9WdvP/hR14nEhE56EQKo8o5N985VwbgnLvVnxc556qA24B3gNXALOfcSjObamZT64ZdCawws+XUnlF1tWul+2Qssg3tbpqLCw2nx3s3sHztRq8jiYgAYP7+Xjazp4ApwHRqi+Zl59z7TZjtmDIyMlxWVpZXb9/kCtYtIvKl8XxJT9rf+jZpyQleRxKRFsDMljnnMhry2hPZwih2ztUABcCPgQkNeUPxT/ypZ3Lg3D+RwSpWzLiB/JIKryOJSCt3IoWRAOCcu905VwiUNEkiOSjlrO+S3f8HXFj5Hm/9/V4qq2u8jiQirdiJFMa/zWy6mfU0s25Ap6YKJf+VevmDbO94AVcfeIqXn/+7TrcVEc8ctzDM7GUzu4va6zBeAm4H7gUeaNpoAkBICGnX/y+7Y3tz1Zb7eHXuHK8TiUgr5c8WxnRqdz9dCTwOpAA/cc6tacpgcojwaDpMnUtRRDLnL/8Bb70/3+tEItIK+VMY6cAHwNS6mwPORlsXAWdx7UmY8i8IjWDQwu/z4dLPvY4kIq2MP4VxKvArYI2ZfU7trUGuMLMxZtay7s8R5HxJ3Qi/fi7xIeV0evMalq/d4HUkEWlF/CmMx51zlzrnTgHOAV4AooHJwFtNGU6OFJ02kOpvv0Sa5RLy8kQ27tjtdSQRaSX8KYz/mNkWM3sb+CW1twh52Tn3fefcsKaNJ/Vp03s0+RfNoC8byX3q2+zULdFFJACOWxjOucFAD2rPjloIPEFtaYiHkoddwa5vPcJwt5zVf7uW/OJyryOJSAvn13UYzrlq59xK59wc59ybdfeHEo91HnMzWwfdwTmVH7Doie9TXFbpdSQRacEC8p3e0nS6jP8Zm069gXGlb7LwLzdQVqEuF5GmocJo7szoPumPbOhxHRcWv86Cv0yhvFKlISKNT4XREphxynf+zLqu32Fs4WsseOIWKqv0Na8i0rhUGC2FGade9wRr0idyQf4s5j95K9W6WaGINCIVRktiRq/rp7Oy81Wcv/9l5v/1B9SoNESkkagwWhoz+t4wg686XM65eS/wwd9/hKtRaYjIyVNhtEQhIfSb8jRfpIzn7D3PsfBvt2lLQ0ROmgqjhbKQUAZMfZbPUq5gVO6LLHniOqoqdZ2GiDScCqMFs5BQBk99mqVp1zNi/xss//NVVJSXeR1LRJopFUYLZyEhDLvhMT455UdkFM1nzWOXUFZS6HUsEWmGVBitxOnfeYCP+95Hv5KlbH1sLCUFeV5HEpFmRoXRipxx1U/5JOOPdCtfze7Hz6UoL8frSCLSjKgwWpkzLrmBz8/8Gx0qsznw13PZn73W60gi0kyoMFqh08+/mhVjniOmKh+eOo+crz7wOpKINAMqjFZq2Mhx7LjyDQqJot1rV7JxwQteRxKRIBeQwjCzsWa21sw2mNnd9SyfbGZf1v1ZbGYDA5Grtes3YCjc8C7rQk6hx4JbWTv7QXDO61giEqSavDDMLBR4ErgQ6ANMMrM+hw3bDIxyzg0AHgJmNHUuqZWelk6nH77DhxGjOG3FH1n7P9dDtS7wE5EjBWILIxPY4Jzb5JyrAGYC4w8d4Jxb7JzbX/dwCZAagFxSp11CG4b99DXeSpjMaTn/ZONj46gpOeB1LBEJMoEojM7A9kMeZ9c9dzQ3AG/Xt8DMpphZlpll5ebmNmJEiQz3ccEPn+T1LveSXrCMXY+OpGSnzqASkf8KRGFYPc/Vu6PczM6mtjDuqm+5c26Gcy7DOZeRnJzciBEFIDTEGH/9Xbw7dDqRFXm4v49m19K5XscSkSARiMLIBtIOeZwKHHHFmJkNAJ4CxjvndBmyh8Zd+m02Xf4WW2lPh39dx+bZvwDdIl2k1QtEYSwFeppZNzMLByYCbxw6wMzSgTnAtc65dQHIJMeRMWgg8dPe573wc+i24nE2P3GpjmuItHJNXhjOuSrgNuAdYDUwyzm30symmtnUumH3Ae2Av5rZcjPLaupccnypKe0486evMrv9D0nNW0zuoyMoyV7hdSwR8Yi5ZnrefUZGhsvKUq8EgnOOt96cQ2bWj4m1cvLPf5QOI67xOpaINICZLXPOZTTktbrSW47LzLjokivZNuFt1pFOh//cwqZnbsBVFHsdTUQCSIUhfhvavy/tfziPubHfpvu22ez6wwiKt3/pdSwRCRAVhpyQDolxXPKTGbze/wnCyvcT9vQYtv3nL7qliEgroMKQExYaYoy/8lpyJr3L8pA+pC/+ORufvJKa4v3Hf7GINFsqDGmwgb1Oo9ft/+GfSTeTnruAfX/KJG/lfK9jiUgTUWHISWkTHcFlt/6W+Wc+T0kVtH31ctY9d5sOiIu0QCoMOWlmxvnnX4y7eRHvRl/EqZv/wa7fZZK76kOvo4lII1JhSKPp0imF825/gXeGTMdVlpH4yiWsfv5HuMpSr6OJSCNQYUijCgkxLrh0ElVTF7Eg+gJ6b3qWHb89nT1rP/Y6moicJBWGNIn0jh04+46ZvD/kSXyVBSS+NI6vnv0BVaWFXkcTkQZSYUiTCQkxxlz6HSpvXsxHsRfQf+vz5P1+MBs+etXraCLSACoMaXKpnTox6vaXWTLqRYpdJKe8dyNf/elS9u/a6nU0ETkBKgwJCDNj+NkX0/6OT1mQegs98xfjm346y2b9hpqqKq/jiYgfVBgSUDHR0Yy+8RF2fmcBGyL6MHTVI2x+5HTWfPqu19FE5DhUGOKJbj37MfCu9/hk6B+IrdpPr7cmkPXHy8nZqu/PEglWKgzxjIWEcPolNxF3++csSbuRfgUfkvjMCD6a8WPy8w94HU9EDqPCEM9Fx7Zh+A1/pPCmJaxuM5Kzcp6h7NHBfDj7CSp1fEMkaKgwJGgkp57C4J/MYdOlr1Hsa8e3VvyMDQ8P56P/zKa6RrdPF/GaCkOCTvch59Ltnk9YmfkI7dx+zlp8A188PIrFC/9Dc/1KYZGWQIUhQclCQuk77haS7lnBygF30716MyPmXcXi31zEp0s/VnGIeECFIUEtJDyKvlfcQ+wdX7Gy5y0MrljG0DcvZP7vvk3W8uUqDpEAsub6Fy4jI8NlZWV5HUMCrCJ/DxvnPECPrTMx5/ggagxRZ9/BiMxhmJnX8USCnpktc85lNOi1KgxpjsrytrLl9d/QbdtswlwVC8JHETLqDkaOOJPQEBWHyNGoMKTVqjyQw+Y3HiFt00wiXAUf+EZQMeKnnD3ybMLDtMdV5HAqDGn1qgtz2fLmb+m49gWiKeUjG8Le/jcx+oIrSYiJ8DqeSNBQYYjUcSX72Pr2YySu/F/ia/JZ5bqyqut3GXbRDXRJSfA6nojnTqYwArLNbmZjzWytmW0ws7vrWd7LzD42s3Izuz0QmaRlsuhEul75IPH3rGXnyN/SLqKGCVsfxPfkYGb/5U6Wrd2iM6tEGqjJtzDMLBRYB5wHZANLgUnOuVWHjEkBugCXAfudc3843nq1hSF+qanhwJf/omDeo6QXLKPIRTI/4mxCMm9i9LdGERMR5nVCkYA6mS2MQPxtyQQ2OOc2AZjZTGA8cLAwnHN7gD1mdlEA8khrEhJCwqBLSBh0CaXblrHnnce4YMfbhH/0Nlkf9mZz92sYcsF36NEh0eukIkEvELukOgPbD3mcXffcCTOzKWaWZWZZubm5jRJOWo+o9KF0v+kf+O5YS3bGPXQNz+eqzb8g7m+DmfPHacz/5HMqqmq8jikStAJRGPWdFN+g/WDOuRnOuQznXEZycvJJxpLWymLakXrx3STds5L8K16mpF0/Lit8iZFvnc3SX53Na/94gg0793odUyToBGKXVDaQdsjjVCAnAO8rcmwhIbQZMI42A8ZRnbeZ7fP+h95rZnHmxp+xb8NveDN6DL6M73LWWaN1rEOEwBTGUqCnmXUDdgATgWsC8L4ifgtt142uVz0MNQ9xYOV/2P/h05y/51+Ef/g6KxZ2Z0PHi+l05mSG9j1NV5JLqxWQ6zDMbBzwGBAKPOOc+7WZTQVwzk03sw5AFhAP1ABFQB/nXMHR1qmzpKSpueI8ti98jtAvXqJz2XqqXAhLQwawt9ulnDpqEqd16eR1RJETpgv3RJpYec4Ktn3wPAkbXye5ahdlzscn4adT0utK+o+8nNTktl5HFPGLCkMkUJyjYP0icj58no7Zb9PGFVDoovgs8nQqel5M31FX0Cm5ndcpRY5KhSHihepK9nzxDnuXzqbzrvdp4woocRF8HjGM8lMvptfICXRK0dl8ElxUGCJeq65i91fvkfvpq3Te+R5t3QHKnY+vwgdQmH4unU+/gp49T9N3dojnVBgiwaSmmp0rPmDPp7NJ2TmPjtU7AVhvXclpP5q2gy6ld8YofGE6VVcCT4UhEqycY9/WFWxbMofoze/So2wFoebIow0b4obhuo+h+xmXkNIh7fjrEmkEKgyRZqLkwB7WL/onNevepWv+J7Sl9szxDaE9yG1/Jm36jaVnxhh84ZEeJ5WWSoUh0gy5mmo2f7WY3OVvEbdjIT3LV+OzakpdOBuj+lOaeibJ/c+lS78RWKjP67jSQqgwRFqAgvx9rP/4X5RvWEDHfZ/SrWYbAEVEsSVmEJXpZ9Kh/xg6nDZMBSINpsIQaYF27tjGlmX/pmbTQlIPZNGF2oPnJUSyLbov5R0zSewzms59zyIkMtbjtNJcqDBEWjjnHJs3rWfHl/Op2foxHfM/55SarYSYo5JQssNPoTB5EFHdhpPafyRRKT1Ap/BKPVQYIq2Mc44du3azefl8KjYtol3e55xavZ5oKwfggLVhV1w/qjsNJfG0EXQ47XQsWl8SJSoMEQH2F5awfsWn5K//mPBdy0gtXkUP23Fw+e7Qjuxr0wc6DaZdz0ySe2Zi0boHVmujwhCRI1TXODZs287O1R9Ttu0z4vK+Ir18PWm25+CYPaEd2Bd3GjXt+xHfdRDtT83El9hFu7NaMBWGiPilvKqajVu2s3PtEiq3LSNm/2o6l2+gK7sIsdrfBUUWw56oHpS17YWvU18Suw4isdsA7dJqIVQYItJgVdU1bNmZS866ZRRvW44vdyXJxevo5rYTb6UHx+0LSSQvugcViacS3qE3bdP70q5rfywmSVskzYgKQ0QaXV5hGVs2r2f/li+o3LmSqP3rSC7bRHeXTZRVHBxXYHHsjehCSZseWFJPYjr1IqlLb2I79ISwCA8/gdRHhSEiAeGcI7eglO1bNnBg21dU7V5L+IENJJRsJq06myT775dk1mDkhqRwICqN8vhuWLvuRHc4hbadT6Vt555YeIyHn6T1OpnC0O0yRcRvZkZKm2hSBg6AgQO+saysspp1OTvI27qakp3rcHkbiSzcQmLJNroWrSJ+Zwms+O/4PGvLvvBOlMSkUR2fhi+xC9Htu9G2Uw/aduiG+XQ/rWCjwhCRRhHpC+XULunQJR244BvLSsur2Lgrh7ztaynZvZGafZvwFWwjrnQHKXlZtM97h9At/93bUeOMvJC27PN1oDiyA1WxnbA2qUS0Syc2pQuJnboTn9gBCwkJ8Kds3VQYItLkoiLC6NElnR5d0oHzvrHMOUdBUSm7czaRn7OR8r2bqdm/HV9hNrFlOSQVrCIlfyEROVXfeF2585EXkkh+WBLFkSlURnXAxXXE17YzkYmdiUvqTEJKGm3aJOiLqxqJCkNEPGVmtImLps1p/eC0fvWOqaqqYveeHPbv2kzR7q1U7t8G+TvwlewmqnwPHYrW0K5gEVF7Ko54bZGLYl9IAgWh7SgJT6IiKpnq6CQstj3hbdoTmdCBmMSOtEnuRNv4OHyh2mo5GhWGiAS9sLAw2ndKp32n9KMPco6i/Dz27dxC0d5syg7kUJ2/CyvaRVjJHqLK99K5bC1tSj4hNq+03lUUuCh2WjyFIQmUhCVQFt6WyshEaqKSCIlpR2hsEuFxSUS1SSY2MYX4hHYkREcSHtY6SkaFISItgxmxCUnEJiQBxz4JqKa8hKJ9OeTvzaFk3y7K83dRXbAbSvIILd2Lr3w/KZW5xBSvp01RAT6q6l1PlQvhALEUEEtxSByloXGU+eKpDG9DdUQCLjIBi0wgLCaBsJhEImMTiYxPJDq+HXFx8cRF+YjyhTabXWYqDBFpdUIioonveArxHU85/mDncGX5lBXspXD/HkoO7KG8YC+Vhbm44jxc6X5Cyw4QXnGAuMoCoiqyiSktJJbiY6620oVSSBS5xFBs0ZSFxFIWGktFWCxVvlhqfHG4iFiIiMci4wmNiscXHY8vKp7w6HgiY9oQFdeGmOhYYiJ9RPtCCQlp2uJRYYiIHIsZFpVAVFQCUe39KJivVVdBeQHlRXmU5udRUriP8oI8Kov3U1W8H1eWD2UFWEUBoRUFJFQWEl61k8jyYiLLSoh2JYRw/OvkqlwIxUSykyjKiKQsJIpyi6IyNIrK0Giqw6KoDoumJiwafNEnMREqDBGRphEaBtGJREQnEpHSk4QTfb1zUFGMKy+govgAJYUHKCs6QGVJARWlBVSVFFBTXkhNWSGUF0FlMSGVxYRVlRBZVYKvOg9fZTYR5aVEuDIiKSeUmpP6SAEpDDMbC/wZCAWecs49cthyq1s+DigBvuec+ywQ2UREgpIZRMRiEbFExHciouNJrs85qCqHB6IavIomP7RvZqHAk8CFQB9gkpn1OWzYhUDPuj9TgL81dS4RkVbFDE7y6vlAnAuWCWxwzm1yzlUAM4Hxh40ZDzzvai0BEszsZPtUREQaUSAKozOw/ZDH2XXPnegYzGyKmWWZWVZubm6jBxURkaMLRGHUd57X4Yf+/RmDc26Gcy7DOZeRnJzcKOFERMQ/gSiMbCDtkMepQE4DxoiIiIcCURhLgZ5m1s3MwoGJwBuHjXkD+K7VGg7kO+d2BiCbiIj4qclPq3XOVZnZbcA71J5W+4xzbqWZTa1bPh14i9pTajdQe1rt9U2dS0RETkxArsNwzr1FbSkc+tz0Q352wK2ByCIiIg3TOm6xKCIiJ02FISIiflFhiIiIX1QYIiLiF6s93tz8mFkhsNbrHH5IAvZ6HcIPytm4mkPO5pARlLOxneaci2vIC5vz7c3XOueO/bVaQcDMspSz8Shn42kOGUE5G5uZZTX0tdolJSIiflFhiIiIX5pzYczwOoCflLNxKWfjaQ4ZQTkbW4NzNtuD3iIiEljNeQtDREQCSIUhIiJ+CfrCMLOxZrbWzDaY2d31LDcze7xu+ZdmNiRIc442s3wzW1735z4PMj5jZnvMbMVRlgfLXB4vZzDMZZqZzTez1Wa20sz+Xz1jPJ9PP3MGw3xGmtmnZvZFXc4H6hkTDPPpT07P57MuR6iZfW5mb9azrGFz6ZwL2j/U3g59I9AdCAe+APocNmYc8Da139o3HPgkSHOOBt70eD5HAkOAFUdZ7vlc+pkzGOayIzCk7uc4YF2Q/r/pT85gmE8DYut+9gGfAMODcD79yen5fNbl+AnwUn1ZGjqXwb6FkQlscM5tcs5VADOB8YeNGQ8872otARLMrGMQ5vScc24hsO8YQ4JhLv3J6Tnn3E7n3Gd1PxcCqznye+g9n08/c3qubo6K6h766v4cfkZOMMynPzk9Z2apwEXAU0cZ0qC5DPbC6AxsP+RxNkf+z+7PmKbmb4Yz6jZl3zazvoGJdkKCYS79FTRzaWZdgcHU/mvzUEE1n8fICUEwn3W7UJYDe4B3nXNBOZ9+5ATv5/Mx4E6g5ijLGzSXwV4YVs9zh7e5P2Oamj8ZPgO6OOcGAn8B5jZ1qAYIhrn0R9DMpZnFAq8BP3LOFRy+uJ6XeDKfx8kZFPPpnKt2zg0CUoFMM+t32JCgmE8/cno6n2Z2MbDHObfsWMPqee64cxnshZENpB3yOBXIacCYpnbcDM65gq83ZV3tNxD6zCwpcBH9EgxzeVzBMpdm5qP2l/CLzrk59QwJivk8Xs5gmc9D8hwAFgBjD1sUFPP5taPlDIL5PBO41My2ULt7/Bwze+GwMQ2ay2AvjKVATzPrZmbhwETgjcPGvAF8t+6o/3Ag3zm3M9hymlkHM7O6nzOpnfu8AOc8nmCYy+MKhrmse/+ngdXOuT8dZZjn8+lPziCZz2QzS6j7OQo4F1hz2LBgmM/j5vR6Pp1z9zjnUp1zXan9XTTPOfedw4Y1aC6D+m61zrkqM7sNeIfaM5Gecc6tNLOpdcunU/td4eOADUAJcH2Q5pwA3GJmVUApMNHVna4QKGb2MrVncCSZWTbwS2oP2gXNXPqZ0/O5pPZfcdcCX9Xtzwa4F0g/JGcwzKc/OYNhPjsCz5lZKLW/YGc5594Mtr/rfuYMhvk8QmPMpW4NIiIifgn2XVIiIhIkVBgiIuIXFYaIiPhFhSEiIn5RYYiIiF9UGCIi4hcVhoiI+EWFIdJE6i7uEmkxgvpKb5HmxsxepfYuoIOB94FfeZtIpPGoMEQaV39q79t0ttdBRBqbbg0i0kjMLBLYBnRyzlV5nUeksekYhkjj6UvtV12qLKRFUmGINJ7+wJdehxBpKioMkcajwpAWTccwRETEL9rCEBERv6gwRETELyoMERHxiwpDRET8osIQERG/qDBERMQvKgwREfHL/wfrlkZpsTSKnQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the ground state wavefunction\n",
    "\n",
    "plt.plot(rr, np.abs( χrα@UU[:,0] )      , label=\"variational\")\n",
    "plt.plot(rr, np.exp(-rr)/np.sqrt(np.pi) , label=\"exact\")\n",
    "    \n",
    "plt.xlim([0,4])\n",
    "plt.xlabel(r\"$r$\")\n",
    "plt.ylabel(r\"$\\psi_{GS}$\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "optimum-match",
   "metadata": {},
   "source": [
    "# The Helium atom\n",
    "\n",
    "Interaction matrix elements\n",
    "\n",
    "![4.17](pngs/4_17.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "id": "terminal-bailey",
   "metadata": {},
   "outputs": [],
   "source": [
    "αα = (0.298073, 1.242567, 5.782948, 38.47497)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "id": "differential-attendance",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_Qprqs(αp,αr,αq,αs): # overlap\n",
    "    return 2*np.pi**2.5/( (αp+αq)*(αr+αs)*np.sqrt(αp+αq+αr+αs) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "id": "accessory-condition",
   "metadata": {},
   "outputs": [],
   "source": [
    "HM = np.array( [[   calc_Tpq(αp,αq) + 2 * calc_Apq(αp,αq) for αp in αα] for αq in αα ] )\n",
    "SM = np.array( [[   calc_Spq(αp,αq)                       for αp in αα] for αq in αα ] )\n",
    "QT = np.array( [[[[ calc_Qprqs(αp,αr,αq,αs)               for αs in αα] for αq in αα ]\n",
    "                                                          for αr in αα] for αp in αα ] )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "greenhouse-horse",
   "metadata": {},
   "source": [
    "![4.18](pngs/4_18.png)\n",
    "![4.20](pngs/4_20.png)\n",
    "![4.21](pngs/4_21.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "id": "floating-liverpool",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "$E_0=$ -2.85516023441995\n"
     ]
    }
   ],
   "source": [
    "tol = 1e-6\n",
    "Niter = 100 \n",
    "EGS = 1e3\n",
    "\n",
    "CV = np.array( [1/np.sqrt(len(αα))]*len(αα) )\n",
    "\n",
    "for it in range(Niter):\n",
    "    \n",
    "    FM = HM + np.einsum('prqs,r,s',QT,CV,CV)\n",
    "    \n",
    "    EE, UU = eigh(FM, SM)\n",
    "    \n",
    "    CV = UU[:,0]\n",
    "    \n",
    "    E0 = 2*np.einsum('p,pq,q', CV, HM, CV) + np.einsum('prqs,p,q,r,s', QT, CV, CV, CV, CV)\n",
    "    \n",
    "    #print(f\"{it:2}, {EGS:10.8}, {E0:10.8}\")\n",
    "    if np.abs(EGS-E0)<tol:\n",
    "        break\n",
    "        \n",
    "    EGS = E0\n",
    "\n",
    "# print Helium atom ground state energy\n",
    "print(r\"$E_0=$\"+f\" {EGS}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "initial-stress",
   "metadata": {},
   "source": [
    "# The H2+ ion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "fundamental-dodge",
   "metadata": {},
   "outputs": [],
   "source": [
    "def F0(t):\n",
    "    from scipy.special import erf\n",
    "    if t==0:\n",
    "        return 1.\n",
    "    sqt = np.sqrt(t)\n",
    "    return  np.sqrt(np.pi)*erf(sqt)/2/sqt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "vocal-enlargement",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.0, 0.9999966666766666)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F0(0), F0(1e-5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "id": "lucky-trademark",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_Rp(α,β,RA,RB):\n",
    "    return (α*RA+β*RB)/(α+β)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "id": "conditional-relations",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_SαβRARB(α,β,RA,RB):\n",
    "    αpβ=α+β\n",
    "    return (np.pi/αpβ)**1.5*np.exp(-α*β/αpβ*(RA-RB)**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "id": "stock-temple",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_TαβRARB(α,β,RA,RB):\n",
    "    αβ=α*β\n",
    "    αpβ=α+β\n",
    "    αβoαpβ = αβ / αpβ\n",
    "    RAB2 = (RA-RB)**2\n",
    "    return αβoαpβ * (3-2*αβoαpβ*RAB2) * (np.pi/αpβ)**1.5 *np.exp(-αβoαpβ*RAB2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "id": "automated-acquisition",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_AαβRARB(α,β,RA,RB,RC,Z=1):\n",
    "    αpβ=α+β\n",
    "    RP = calc_Rp(α,β,RA,RB)\n",
    "    return -2*np.pi*Z/αpβ*np.exp(-α*β/αpβ*(RA-RB)**2)*F0(αpβ*(RC-RP)**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "id": "negative-basic",
   "metadata": {},
   "outputs": [],
   "source": [
    "αα = (13.00773, 1.962079, 0.444529, 0.1219492)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "id": "greenhouse-limitation",
   "metadata": {},
   "outputs": [],
   "source": [
    "RR = [0, 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "id": "ranging-explanation",
   "metadata": {},
   "outputs": [],
   "source": [
    "ααRR = [ [α,R] for R in RR for α in αα ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "id": "alien-thirty",
   "metadata": {},
   "outputs": [],
   "source": [
    "SM = np.array( [ [calc_SαβRARB(α,β,RA,RB) for α,RA in ααRR] for β,RB in ααRR ] )\n",
    "\n",
    "HM = np.array( [ [calc_TαβRARB(α,β,RA,RB)+\n",
    "                  calc_AαβRARB(α,β,RA,RB,RR[0],Z=1)+\n",
    "                  calc_AαβRARB(α,β,RA,RB,RR[1],Z=1) for α,RA in ααRR] for β,RB in ααRR ] )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "id": "million-million",
   "metadata": {},
   "outputs": [],
   "source": [
    "EE, UU = eigh(HM, SM)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "id": "handmade-contrast",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "$E_0=$ -1.4424552955482246\n"
     ]
    }
   ],
   "source": [
    "# print H2+ ground state energy\n",
    "print(r\"$E_0=$\"+f\" {EE[0]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "united-preference",
   "metadata": {},
   "source": [
    "# The H2 molecule"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "id": "sunset-marking",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_QαβRARB(α,β,γ,δ,RA,RB,RC,RD,Z=1):\n",
    "    αpγ=α+γ\n",
    "    βpδ=β+δ\n",
    "    RP = calc_Rp(α,β,RA,RC)\n",
    "    RQ = calc_Rp(α,β,RB,RD)\n",
    "    return 2*np.pi**2.5/( αpγ*βpδ*np.sqrt(αpγ+βpδ) ) *\\\n",
    "             np.exp(-α*γ/αpγ*(RA-RC)**2 - β*δ/βpδ*(RB-RD)**2)*\\\n",
    "             F0(αpγ*βpδ/(αpγ+βpδ)*(RP-RQ)**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "id": "outside-nebraska",
   "metadata": {},
   "outputs": [],
   "source": [
    "αα = (13.00773, 1.962079, 0.444529, 0.1219492)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "id": "collected-classic",
   "metadata": {},
   "outputs": [],
   "source": [
    "RR = [0, 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "id": "emotional-marine",
   "metadata": {},
   "outputs": [],
   "source": [
    "ααRR = [ [α,R] for R in RR for α in αα ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "id": "gothic-lunch",
   "metadata": {},
   "outputs": [],
   "source": [
    "SM = np.array( [ [calc_SαβRARB(α,β,RA,RB) for α,RA in ααRR] for β,RB in ααRR ] )\n",
    "\n",
    "HM = np.array( [ [calc_TαβRARB(α,β,RA,RB)+\n",
    "                  calc_AαβRARB(α,β,RA,RB,RR[0],Z=1)+\n",
    "                  calc_AαβRARB(α,β,RA,RB,RR[1],Z=1) for α,RA in ααRR] for β,RB in ααRR ] )\n",
    "\n",
    "QT = np.array( [ [ [ [ calc_QαβRARB(α,β,γ,δ,RA,RB,RC,RD,Z=1)\n",
    "               for α,RA in ααRR] for β,RB in ααRR ] for γ,RC in ααRR] for δ,RD in ααRR ] )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "id": "choice-greek",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "$E_0=$ -1.0780816790168988\n"
     ]
    }
   ],
   "source": [
    "tol = 1e-6\n",
    "Niter = 100 \n",
    "EGS = 1e3\n",
    "\n",
    "X = RR[1]-RR[0]\n",
    "\n",
    "CV = np.array( [1/np.sqrt(len(ααRR))]*len(ααRR) )\n",
    "\n",
    "for it in range(Niter):\n",
    "    \n",
    "    FM = HM + np.einsum('prqs,r,s',QT,CV,CV)\n",
    "    \n",
    "    EE, UU = eigh(FM, SM)\n",
    "    \n",
    "    CV = UU[:,0]\n",
    "    \n",
    "    #E0 = CV.T@(FM+HM)@CV+1/X\n",
    "    E0 = 2*np.einsum('p,pq,q', CV, HM, CV) + np.einsum('prqs,p,q,r,s', QT, CV, CV, CV, CV)\n",
    "    \n",
    "    #print(f\"{it:2}, {EGS:10.8}, {E0:10.8}\")\n",
    "    if np.abs(EGS-E0)<tol:\n",
    "        break\n",
    "        \n",
    "    EGS = E0\n",
    "\n",
    "# print Helium atom ground state energy\n",
    "print(r\"$E_0=$\"+f\" {EGS+1/X}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "id": "hybrid-change",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.09275446, 0.16665171, 0.11588497, 0.02540915, 0.09329479,\n",
       "       0.16302778, 0.12491446, 0.01671247])"
      ]
     },
     "execution_count": 128,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "CV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "suffering-liquid",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "coated-motel",
   "metadata": {},
   "source": [
    "# Car-Parrinello for HF orbitals\n",
    "\n",
    "Start with the Langrangian\n",
    "\n",
    "![9.20](pngs/9_20.png)\n",
    "\n",
    "which leads to \n",
    "\n",
    "![9.21-22](pngs/9_21_22.png)\n",
    "\n",
    "In terms of the time-dependent coefficients of the wavefuctions\n",
    "\n",
    "![9.18](pngs/9_18.png)\n",
    "\n",
    "the first equation above becomes\n",
    "\n",
    "![9.23](pngs/9_23.png)\n",
    "\n",
    "Update coefficients according to\n",
    "\n",
    "![9.30](pngs/9_30.png)\n",
    "\n",
    "Do it in two steps: First calculate\n",
    "\n",
    "![9.31](pngs/9_31.png)\n",
    "\n",
    "Including friction this equation becomes\n",
    "\n",
    "$$\n",
    "(\\mu + \\gamma h/2)\\tilde{\\mathbf{C}} = 2\\mu \\mathbf{C} - (\\mu-\\gamma h/2) \\mathbf{C'} - 4 \\mathbf{F}\\mathbf{C} h^2\n",
    "$$\n",
    "\n",
    "Then find $\\lambda$ from \n",
    "\n",
    "![9.32](pngs/9_32.png)\n",
    "\n",
    "and add the missing term in Eq. 9.30."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "id": "sudden-arlington",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9999999999999998"
      ]
     },
     "execution_count": 130,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "CV = np.array( [1/np.sqrt(len(ααRR))]*len(ααRR) )\n",
    "np.sum(CV**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "id": "solar-apparel",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0000000000000002  should be equal to 1\n"
     ]
    }
   ],
   "source": [
    "tol = 1e-6\n",
    "Niter = 100 \n",
    "EGS = 1e3\n",
    "\n",
    "CV = np.array( [1/np.sqrt(SM.sum())]*len(ααRR) )\n",
    "\n",
    "# check condition \n",
    "print(CV.T@SM@CV, \" should be equal to 1\")\n",
    "\n",
    "CVold = CV.copy()\n",
    "\n",
    "h = 0.1\n",
    "h2 = h*h\n",
    "h4 = h2*h2\n",
    "\n",
    "μ = 1.\n",
    "γ = 1.\n",
    "\n",
    "μmγ = μ - 0.5*γ*h\n",
    "μpγ = μ + 0.5*γ*h\n",
    "\n",
    "Evals = []\n",
    "\n",
    "for it in range(Niter):\n",
    "    \n",
    "    FM = HM + np.einsum('prqs,r,s',QT,CV,CV)\n",
    "    \n",
    "    EE, UU = eigh(FM, SM)\n",
    "\n",
    "    E0 = CV.T@(FM+HM)@CV+1/X\n",
    "    \n",
    "    # CV = UU[:,0]\n",
    "    # instead update CV with the Verlet algorithm\n",
    "    \n",
    "    # calculate Ctilde\n",
    "    CVt = ( 2*μ*CV - μmγ*CVold - 4*FM@CV*h2 ) / μpγ\n",
    "    \n",
    "    # solve for lambda\n",
    "    #a = np.dot(CV,SM@CV)*h4\n",
    "    a = (SM@CV).T@SM@SM@CV   \n",
    "    b = -CVt.T@SM@(SM@CV)-(SM@CV).T@SM@CVt  \n",
    "    c = CVt.T@SM@CVt-1.\n",
    "    \n",
    "    λ = np.roots([a, b, c] )[1]\n",
    "    #print(a,b,c,λ)\n",
    "    \n",
    "    # add the missing term\n",
    "    CVt -= λ*SM@CV\n",
    "\n",
    "    #print(f\"{it:2}, {EGS:10.8}, {E0:10.8}\")\n",
    "    if np.abs(EGS-E0)<tol:\n",
    "        break\n",
    "\n",
    "    # update variables\n",
    "    CVold = CV\n",
    "    CV = CVt\n",
    "\n",
    "    EGS = E0\n",
    "    Evals.append(EGS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "id": "underlying-tuner",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1533d7b59220>]"
      ]
     },
     "execution_count": 132,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAD7CAYAAACFfIhNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAqUklEQVR4nO3de5zcdX3v8ddnbnuZ2fvmshsuCSRCEjCBrsi9AonFVI1QbfGCWNvy6OMhVbHag7Wn2vb0eNf29GGpodrSYw9UEQoiCiEKKMgl3HLbhCQkkGQ3e81m77eZ7/ljfjM72czszuws2Vnm/Xw8xrn9fr/5kJh5z/f3/X2/X3POISIiAuCb6wJERKRwKBRERCRJoSAiIkkKBRERSVIoiIhIkkJBRESS8goFM6s1s81mtte7r8mw3a1mttPMdpjZXWZWmsv+IiJyauTbUrgN2OKcWwFs8Z6fwMyWAJ8Empxz5wF+4IZs9xcRkVMnkOf+G4F3eI/vBB4D/keGzykzszGgHGjJcf8T1NfXu6VLl86sYhGRIvX88893OucWTLVNvqGwyDnXCuCcazWzhZM3cM4dMbNvAK8DQ8AjzrlHst0/naVLl7J169Y8SxcRKS5m9tp020x7+sjMHvX6AibfNmZZRA3xFsEyoBEIm9lHstl30nFuNrOtZra1o6Mj191FRCQL07YUnHPrMr1nZm1m1uD9ym8A2tNstg444Jzr8Pa5F7gU+AGQzf6JOjYBmwCampo0YZOIyBsg347mB4CbvMc3Afen2eZ14GIzKzczA64BmnPYX0RETpF8Q+ErwHoz2wus955jZo1m9hCAc+4Z4B7gBWC795mbptpfRETmhs3HqbObmpqcOppFRHJjZs8755qm2kYjmkVEJEmhICIiSUUVCr/Y3cY/P7ZvrssQESlYRRUKT7zSye2/3D/XZYiIFKyiCoWa8hB9I+OMRWNzXYqISEEqqlCoLg8CcHxobI4rEREpTEUZCj2Do3NciYhIYSqqUKgpDwHQM6iWgohIOkUVComWwjGFgohIWkUVComWwjGdPhIRSauoQiHZ0ayWgohIWkUVCpGSAAGfqaUgIpJBUYWCmVFdHlSfgohIBkUVCgBVZUGOD6mlICKSTtGFQk15iGMDaimIiKRTdKFQXR5Sn4KISAZFGApBTXMhIpJB0YVCTXlQLQURkQyKLhSqy0MMj8UYHovOdSkiIgWnCEMhMSmeTiGJiExWdKGgqS5ERDIrulCYmBRPoSAiMlnxhUJZvKWg+Y9ERE5WdKFQE9b02SIimRRfKKhPQUQko6ILhdKgn5KATwPYRETSKLpQgMT8R2opiIhMVpShoOmzRUTSyysUzKzWzDab2V7vvibDdrea2U4z22Fmd5lZqff6l8zsiJm95N025FNPtqrLg/SoT0FE5CT5thRuA7Y451YAW7znJzCzJcAngSbn3HmAH7ghZZNvO+fWereH8qwnKzXlIXrUpyAicpJ8Q2EjcKf3+E7gfRm2CwBlZhYAyoGWPD83L9XlIbUURETSyDcUFjnnWgG8+4WTN3DOHQG+AbwOtALHnXOPpGxyi5ltM7PvZzr9NNvip4/GcM6dio8TEZk3pg0FM3vU6wuYfNuYzQd4X/QbgWVAIxA2s494b98OnA2sJR4Y35ziODeb2VYz29rR0ZHNR2dUUx5kPOboHxnP6zgiIm82gek2cM6ty/SembWZWYNzrtXMGoD2NJutAw445zq8fe4FLgV+4JxrSznWHcCDU9SxCdgE0NTUlNdP/GpvAFvP4BgVpcF8DiUi8qaS7+mjB4CbvMc3Afen2eZ14GIzKzczA64BmgG8IEm4DtiRZz1ZqS7TpHgiIunkGwpfAdab2V5gvfccM2s0s4cAnHPPAPcALwDbvc/c5O3/NTPbbmbbgKuAW/OsJys14YmWgoiITJj29NFUnHNdxH/5T369BdiQ8vyLwBfTbHdjPp8/UzWaPltEJK2iHNFcVaaWgohIOkUZClqSU0QkvaIMhaDfR0VJQKePREQmKcpQAKjS/EciIicp2lDQ/EciIicr2lDQ9NkiIicr4lDQpHgiIpMVbSjUeJPiiYjIhKINheryEL3DY0RjmilVRCSheEOhLIhzcFydzSIiSUUbCjXhxAA29SuIiCQUbSgkps/WFUgiIhOKNxTK1FIQEZmsaEOhplyT4omITFb0oaD5j0REJhRtKFSUBvCZWgoiIqmKNhR8PqOqLEjPkFoKIiIJRRsKED+FpKuPREQmFHUoaPpsEZETFXUo1JSH1KcgIpKiqEOhukyT4omIpCrqUIiUBugfGZ/rMkRECkZxh0JJgIGRcZzTTKkiIlDkoRAuCTAec4yMx+a6FBGRglDUoVBRGgDQKSQREU9Rh0I45IXCsEJBRASKPBQiaimIiJyguEOhRKEgIpIqr1Aws1oz22xme737mgzbfcrMdpjZTjP7dK77v1ESoTCgUBARAfJvKdwGbHHOrQC2eM9PYGbnAX8CXASsAd5tZiuy3f+NFFZLQUTkBPmGwkbgTu/xncD70myzEnjaOTfonBsHHgeuy2H/N4yuPhIROVG+obDIOdcK4N0vTLPNDuBKM6szs3JgA3B6Dvu/YZItBV19JCICQGC6DczsUWBxmre+kM0HOOeazeyrwGagH3gZyPlb2MxuBm4GOOOMM3LdPa3yoB8z9SmIiCRMGwrOuXWZ3jOzNjNrcM61mlkD0J7hGN8Dvuft87+Bw95bWe3vHWMTsAmgqalpVual8PmMcChAn0JBRATI//TRA8BN3uObgPvTbWRmC737M4Drgbty2f+NlJj/SERE8g+FrwDrzWwvsN57jpk1mtlDKdv92Mx2AT8BPuGcOzbV/qdSuMSvjmYREc+0p4+m4pzrAq5J83oL8Q7lxPMrctn/VIqUBukfic5lCSIiBaOoRzQDREr89A9roR0REVAoeH0KaimIiIBCgXCJVl8TEUko+lCoUCiIiCQVfSgkWgpaklNERKFApDRAVEtyiogACoXk9Nl9mv9IREShoDUVREQmFH0oaE0FEZEJRR8KFQoFEZGkog8FrakgIjKh6EMh4q2+NjCqUBARUSjo6iMRkSSFgq4+EhFJKvpQKA/Fl+RUR7OIiEIBMyMS0vxHIiKgUADinc26+khERKEAxC9L1dVHIiIKBSDe2ayrj0REFApAYvU1hYKIiEKBeCioo1lERKEAeH0KWqdZREShAFBRGqBveGyuyxARmXMKBSBc4mdgNKolOUWk6CkUgEhJkGjMMTymJTlFpLgpFIBIiR/If6qL7z6+n0d2Hp2NkkRE5oRCgYnps/MNhX9+bD/f+eW+2ShJRGROKBSAcCj/mVL7R8Y5PjTG9iPHOT6kTmsRmZ/yCgUzqzWzzWa217uvybDdp8xsh5ntNLNPp7z+JTM7YmYvebcN+dQzU4mWQj6jmlt7hgCIOXj61a5ZqUtE5FTLt6VwG7DFObcC2OI9P4GZnQf8CXARsAZ4t5mtSNnk2865td7toTzrmZHZWFPhsBcKAE/t68y7JhGRuZBvKGwE7vQe3wm8L802K4GnnXODzrlx4HHgujw/d1YlQiGfPoUWLxTOWVTBk/vVUhCR+SnfUFjknGsF8O4XptlmB3ClmdWZWTmwATg95f1bzGybmX0/0+mnN9pshULAZ7zvgiXsa++nrXd4tsoTETllpg0FM3vU6w+YfNuYzQc455qBrwKbgZ8DLwOJb9/bgbOBtUAr8M0p6rjZzLaa2daOjo5sPjprs3H1UUvPMIurSrliRT0AT+3XKSQRmX+mDQXn3Drn3HlpbvcDbWbWAODdt2c4xveccxc6564EuoG93uttzrmocy4G3EG83yFTHZucc03OuaYFCxbk/l86hbKgH5/l16dwpGeIxuoyVjVUUl0e5Ml9OoUkIvNPvqePHgBu8h7fBNyfbiMzW+jdnwFcD9zlPW9I2ew64qeaTjkzI5znmgpHjg2xpLoMn8+49Ow6ntrXqWkzRGTeyTcUvgKsN7O9wHrvOWbWaGapVxL92Mx2AT8BPuGcO+a9/jUz225m24CrgFvzrGfG8llTIRpzHO0dprG6FIBLz66n5fgwBzoHZrNEEZE3XCCfnZ1zXcA1aV5vId6hnHh+RYb9b8zn82dTPmsqtPcNE405llSXA3DZ8ni/wpP7uzhrQWTWahQReaNpRLMnnEcoJC5HTbQUltaV01hVqvEKIjLvKBQ8FaUzD4XDx+KhsKS6DIj3UVy6vJ7fvNpFLKZ+BRGZPxQKnnBo5n0KLT3xMQkNXigAXL68np7BMXa19s5KfSIip4JCwRMpDdA/w6uPWnqGqCoLJgfBAVx6dh0AT+oUkojMIwoFTz4dzS09Q8lTRwkLK0tpqCplT1vfbJQnInJKKBQ8iVCYydiCxMC1yeojJXQPjM5GeSIip4RCwRMuCRBzzGhJziM9QyzxrjxKVRcJKRREZF5RKHiSayqM5LZATu/wGH3D42lbCnXhErr6FQoiMn8oFDyJdZoHRqI57dfqXXmUNhQiITr7RzTdhYjMGwoFT6QkCJDzFUiJgWtLatK1FEKMjMcYGM0taERE5opCwRP2Wgq5XoF0pOfEgWup6iIlAHTrFJKIzBMKBU9FoqUwg1AI+o0FXgCkqouEAOgcGMm/QBGRU0Ch4Akn+xRyP320uKoUn89Oeq8uHA8FdTaLyHyhUPBMXH2Ueyg0Vp186ggmTh919aulICLzg0LBkzh9lHtLYThtJzOktBQ0VkFE5gmFgqc06MNnuV19NB6NcbR3OG0nc/yYfiIlAZ0+EpF5Q6HgMbOc5z9q6xshGnNpxygk1IZDdKmjWUTmCYVCilxDYWJxncyhUBcJqaUgIvOGQiFFpDS3NRWSA9fSzHuUUBcuUZ+CiMwbCoUUuS7JeSSLlkJ9JKSrj0Rk3lAopMj19NGRY0PUlAcpDwUyblMbjs+UqmU5RWQ+UCikqCgN0DuU/SypLRnWUUhVFylhPOboHc5t9lURkbmgUEhRUx7i2GD2X96tx4dpyDBwLaE+orEKIjJ/KBRS1IVD9AyOEs3yVE9H3wgLKk6e8+jEYyZGNSsURKTwKRRS1IRDxBxZnUIaj8boHhxlgdcSyKQ2Of+ROptFpPApFFLU5jAtxbHBMZyD+mlaCvXJmVLVUhCRwqdQSJEIhWOD03+Bd3q//OvTTJmdqsY7ptZUEJH5QKGQoqbc+wLP4ld9tqEQ9PuoLg9qqgsRmRcUCikSLYXcQmHqPoXEcdXRLCLzQV6hYGYfMLOdZhYzs6YptrvWzPaY2T4zuy3l9Voz22xme737mnzqyVdOodAX36ZumpYCQH24JBkiIiKFLN+Wwg7geuCJTBuYmR/4DvAuYBXwQTNb5b19G7DFObcC2OI9nzOlQT/lIT/HsmwphPw+Kkszj2ZOqIuEsgoaEZG5llcoOOeanXN7ptnsImCfc+5V59wocDew0XtvI3Cn9/hO4H351DMbasqz+wLv6B+hPhLC7ORlOCeri4Q0eE1E5oVT0aewBDiU8vyw9xrAIudcK4B3vzDTQczsZjPbamZbOzo63rBi6yIhurO6+mh02stRE2rDJRwbHGU8Gsu3PBGRN9S0oWBmj5rZjjS3jdPtmzhEmtdynh3OObfJOdfknGtasGBBrrtnraY8lNXpo67+kWmvPEqoj4Rwjpym0BARmQvTnhB3zq3L8zMOA6enPD8NaPEet5lZg3Ou1cwagPY8PytvteEQ+zv6p92us3+E1Y2VWR0zMdVF98DotNNiiIjMpVNx+ug5YIWZLTOzEHAD8ID33gPATd7jm4D7T0E9U6oNT99SiMUcXf2jWbcU6iKa6kJE5od8L0m9zswOA5cAPzWzh73XG83sIQDn3DhwC/Aw0Az80Dm30zvEV4D1ZrYXWO89n1O14RADo1GGx6IZtzk+NMZ4zGUfCmFNdSEi88P011NOwTl3H3BfmtdbgA0pzx8CHkqzXRdwTT41zLbEqOZjg6MZp8VOjDmoy2LgWny7xEypaimISGHTiOZJasNBYOoBbB3el/uCLFsK1WVBfJbdoLhsDY1G+dHWQ4yMZ27RiIjkKq+WwptRrdcpfGwg85VCnd6UFdlekurzGbXhkuR++eofGefj//4czx7oZizq+NDbz5iV44qIqKUwSaKlMNUEdl1ZToaXqi4cmpXTR8cHx/jwvz7D868dI1IS4Nf73rgxGyJSfBQKkyT7FKY41dPZP4LfZ1SXBbM+7mxMddHZP8INdzxNc0svt3/4Qt513mKe3NeV9UpxIiLTUShMUl0ewgy6pxho1tk3Sl04hM83/RQXCXWRkrymuhgei/IH3/0NBzr7+d7Hmnjn6sVcvqKe40Nj7Gw5PuPjioikUihMkmgBdE9x+qgzh9HMCXXhUF4zpb7w+jH2dwzw9fev4YoV8RHdly2vB+BXeztnfFwRkVQKhTRqwqFpOppHsr4cNaEuHKJveHzGVwvtaukF4OKz6pKv1UdKWNlQya8VCiIySxQKadSFpz7/39k/mvXlqMljRqa/qmkqza19LKgoOWmajCtW1PP8a8cYGtWlqSKSP4VCGlNNn+2ci0+bneMcRomWxUxPITW39rKq4eS5li5fXs9oNMYzB7pmdFwRkVQKhTRqw5mnz+4bGWd0PJbVMpypElNdzKSzeXQ8xt72PlamCYWLltUSCvh0CklEZoVCIY3EpHjOnXypZ2Kt5Zw7mvOY6mJ/Rz9jUcfKhoqT3isN+nnb0hp+vU+hICL5UyikURsOMR5z9A6Pn/Re5wwGrsHE6aOZjFVIdDJnmqr78uUL2H20j/a+4ZyPLSKSSqGQxlQD2Dr7ZhYKFSUBQn4fHX25txSaW3spCfhYWhdO+/4VK+KXpj61T/0KIpIfhUIatYlf9Wn6FSZaCrn1KZgZi6tKaT2e+6/5Xa29nLu4goA//V/XqoZKasqDGq8gInlTKKRR67UUutNMYNfRP4pZ/BRTrhqqSmk9PpTTPs45mlt703YyJ/h8xqXL6/n1vo60/SAiItlSKKSR+MLP1FKoKQ9l/NU+lcbqMlp6cmspHO0d5tjgGKumWfrziuX1tPWOsK99+qVERUQyUSikkQiFTH0KuZ46SmioKuVo73BOE9g1t8Y7madqKQC8bVktAC8d6plRbSIioFBIqzzkJxTwpb1SaCbzHiU0VpcRjbmcOpsTVx6du/jky1FTLasLUx7ys9PbXkRkJhQKaZgZtRlGNXcNjOYRCqUAtOTQr9Dc2scZteVUlE49TbfPZ6xsqEyGiIjITCgUMqgJhziWrk+hb+YthcSaz6059CvsyjC9RTqrGirZ1dpLTOsriMgMKRQySDcp3tBolIHRaM4zpCY0eqHQ0pNdS2FgZJyDXQPT9ickrG6spH9knEPHBmdUn4iIQiGDmjShkBijkOsMqQmVZQHCIX/Wp492H+3DOdJOb5HO6sYqgFnpVxiPxnitayDv44jI/KJQyKC2PHhSKHQkBq5VzKylYGY0VJdlffooceXRdJejJqxYFMHvs1lZie1LP9nJVd94jGde1ShpkWKiUMigNlxC7/A4Y9FY8rWZTnGRqqGqNOuWwq7WXipLAyypLstq+9KgnxULI3l3Nr/eNcjdzx7CAZ+6+6Up16sWkTcXhUIGteH41T6pnc2dM5whNdWSHAawJUYym2W/FvSqxsq8Tx/945a9+H3G925qontglM/+6GWNlBYpEgqFDGqSA9gmVkpLTHs9045miF+B1Nk/Mu2ynNGYY3dr+jUUprKqoZL2vpEZTbwHsK+9n/tePMxHLzmTq89dxOc3nMuW3e3825MHZ3Q8EZlfFAoZJKe6GEhtKYxQWRqgJOCf8XEbvLEKR6eZGO/17kGGxqJZX46akOhs3tU6s9bCPzz6CqVBP3/622cD8LFLl7Ju5UK+/LNmth/Ov69CRAqbQiGD9KEwmvMynJNNXJY6dSjsznJ6i8kSITKTzubm1l4e3NbKxy9bllwUyMz4+vvXUB8p4VP/9aLGQIi8yeUVCmb2ATPbaWYxM2uaYrtrzWyPme0zs9tSXv+SmR0xs5e824Z86plNyZlSvT6F4bEozUd78+pPgIlRzdPNltp8tA+fxa8oykVVeZDTaspm1K/wrc2vUFka4E+uPOuE12vCIT77znN4tWOA7UfUWhB5M8u3pbADuB54ItMGZuYHvgO8C1gFfNDMVqVs8m3n3Frv9lCe9cyampRJ8cajMW75fy9yoHOAj1+2NK/jNmQ5gK25tZdl9WFKg7mfqlrdmPt0F9sO97B5Vxs3X3kWVWUnT6lxzcqF+H3Go81tOdcjIvNHXqHgnGt2zu2ZZrOLgH3OuVedc6PA3cDGfD73VAj6fVSUBujqH+EL9+3g0eY2vvSe1Vx7XkNexy0L+akpD9IyTZ/C7qO9nJvjqaOEVQ1VHOwaoH/k5OVEM7nn+cOUBf187LJlad+vLg/RdGYNm3cpFETezE5Fn8IS4FDK88Peawm3mNk2M/u+mdVkOoiZ3WxmW81sa0dHxxtV6wlqwyHuef4w/7X1EH929XJuunTprBy3sbqM1ilaCn3DYxzqHmLlNDOjZrK6sRLnJvolpuOcY0tzO5evqCdSEsi43fpVi9h9tI9D3ZpGQ+TNatpQMLNHzWxHmlu2v/bTXWSf6K28HTgbWAu0At/MdBDn3CbnXJNzrmnBggVZfnR+asMhBkajfPCi0/nM+rfM2nEbqqYeq/BKWx8A5y6eWUth9ZJEZ3N2ofBKWz9Heoa45tyFU263buUiAJ1CEnkTy/yz0OOcW5fnZxwGTk95fhrQ4h07+e1iZncAD+b5WbNq3cpFvGVhBX+38bycBpBNp7G6lGcOZJ4+ornVC4Us5zyabHFlKTXlwaz7FRJf8ldNEwpL68OsWBjh0eY2/jDDaSYRmd+mDYVZ8BywwsyWAUeAG4APAZhZg3Ou1dvuOuId1wXjE1ctf0OO21BVRt/wOH3DY2nXSdh9tJeKkuynt5jMzFjdWMXO1uyuFPrF7nbOX1LFosrSabddt2oRdzzxKseHxtJ2SIvI/JbvJanXmdlh4BLgp2b2sPd6o5k9BOCcGwduAR4GmoEfOud2eof4mpltN7NtwFXArfnUM19MXJaa/hTS7tY+zm2oyKt1srqxkleO9p8wd1M63QOjvPD6Ma6eppWQsG7lIsZjjsf2tM+4NhEpXHm1FJxz9wH3pXm9BdiQ8vwh4KTLTZ1zN+bz+fNVY/XEZalvWXTiKSLnHLuP9nH9hUvS7Zq1VY2VjEZjNLf28tbTqjNu99iedpyLX3KajbWnV1MfCfFoczsb1+ZXo4gUHo1ongMNVZlbCoePDdE/Mj7jTuaEy5bX4/cZP9txdMrttuxuZ0FFCed502NMx+8zrjl3EY/taWd0fOpWiIjMPwqFObCoshSfkfay1N1H8+tkTqiPlHDp2XX85OWWjDOcjkVjPLGng6vPWYjPl/2pqnWrFtE3PM5zB7tnXN/IeJSXDvXww+cO8b8e3MVN33+WO554dcbHE5HZcSo6mmWSoN/HwopSjqS5LDUxtuCcRfmFAsB71zTyuXu28eKhHi484+QhIM8d7KZvZJyrszx1lHD58npKAj4272rjsuX1Odc1Mh7lA//yG7Z5E+yVBHwsqizl8Vc6GI3G3rAOfhGZnloKc6ShujTt/Ee7j/ZxZl054SkGkWXrd85bTCjg44GXWtK+/4vmdkJ+H5fn+MVeFvJzxYp6Nu9qm9EEed94eA/bDh/nbzeu5rHPvoNdf3stj332HbxvbSNff3gP//GbgzkfU0Rmh0JhjjRWlaXtU2g+2su5MxzJPFllaZCrz1nIT7e3Ek3z5b1ldzsXn103owDauHYJR3qGeGJvbqPLf723kzt+dYAbLz6Tj16ylKX1Yfw+w+czvv6BNaxftYi/vn8nP37+cM41iUj+FApzpLG6lJaeoRPO9w+NRjnYOZB3J3Oq965tpKNvhKcnrbX8akc/BzoHph3FnMnvrF5MfaSEHzz9Wtb7HBsY5TM/fInlCyP85YaVJ70f9Pv4pw9ewGXL6/jcPS/zyM6pO8mz9WpHP1/+WTP3PH+YgRzmgxIpRgqFOdJQVcbIeOyE9Rr2tvcRc7Ayz07mVFefu5BISeCEU0ixmOPLP9uNz7K/FHWyUMDHDW87nV/sbufINDO+QvxS29vu3caxwVH+8Ya1lIXSz/5aGvSz6cYmzl9SxWd++DIHOwdmVB/Agc4BPvNfL7HuW4+z6YlX+eyPXuZtf/8on/nhSzy1v1NLjIqkoVCYI+kGsO1uzW/Oo3RKg37euWoRP9vRmlwC9Ks/383mXW38z3ev4rSa8hkf+4NvPwOAu555fdptf7j1EA/vbONzv3NOcnW4TMIlAb7z4Qvx+4xb7nph2qVLJxsei/L5e7ex7luP89COVv74irN49i/Xcc+fXsLGtY1s3tnGh+54hj//0cvTDu4TKTYKhTmSGMCW+iu7+WgvZUE/Z9TO/Is6nfesbaR3eJwnXunk7mdf57tPvMpHLj6Dj+U56+uS6jKuPnchdz93aMoxC4e6B/nbn+zikrPq+OPLz8q4XarTasr5+vvfyo4jvXz5od1Z1zQwMs4f/ttz3P3cIT56yZk88RdX8ZcbVrKgooSmpbV8+fq38txfreOT16zg3heO8Ed3bs1pinGRNzuFwhxJLLbz6K625FTUu1v7OGdxRU5jBrJx+fJ6asqDfHvzK/zVf+/gihX1fOk9q2dlkr+PXHwmnf0jPJzh/H8s5vjzH72Mz4xv/P6anP7b3rl6MR+/bBn//tRBfr6jddrtjw+NceP3nuHZg9186/fX8MX3rGZhxcnzOZUG/Xxm/Vv46u+dz5P7Orlh029o75t6fQuRYqFQmCN14RAXnFHNj54/zBVf+yXrv/U4Lx3qmdX+hISg38eG8xvY5a3m9p0PX0jAPzt/9VeuWMAZteUZO5y//+QBnj3QzRffu3pGE/zd9q5zWXNaFZ+7Zxuvd2Vex6Grf4QP3fE0248c5zsfupDrLjht2mP/wdvO4F8/2sT+9gF+7/ancl6tLtVYNMbAyDjHB8fo7B+hd3hsxscSmUs2Hzvbmpqa3NatW+e6jFlxoHOAX+xu5xe723ju4DG++YE1vGdN46x/zr72fr78UDNfeu9qTp/l01P/8vh+vvKz3Wy+9UpWpAy629vWx+/+06/57bcsYNONvzXjlsmh7kE2/J9fAXDrurdw4yVnEvRCLRpzPLithW9tfoWjx4f57o2/xTvOya3z/OVDPfzxf2zl+OAYf3HtOXz8smVZtWg6+kZ4aHsrP3m5ha2vHTvhPTO44PRq1q9azPpVi1i+MLe1tjMZj8Y40jPEyHiMhRUlVJUFZ3Vad3lzM7PnnXNNU26jUCgczrl5+Q+8e2CUi7+8hWtXL+bPrl7OGXXl+My4/p+f4kjPEI/ceiX1kZK8PmN/Rz9femAnv9rbyVsWRfir311F6/Ehbn9sPwe7BlmxMMLfX3c+Fy2rndHxu/pHuO3e7Wze1cbly+v55u+vSTuVeFf/CI/sauOn21p5an8nMQfnLq7gmpULqSwNEvD7CPqN7oFRtjS3s/1IfNT22QvCXH/habzvgiVZt5hiMcf2I8d5/JUOth3u4dXOAQ51DzIWnfg3Gwr4WBAp4eyFEd6+rJaLltXy1tOqKAnkvra3vPkpFOSU+cJ92/lP7yokn0FtuITO/hH+5SMX5r2udYJzjs272vi7n+7iUHe8g/68JZXcctVy3rlqcd59Mc457nr2EH/34C4CfmPt6dUsqS7jtJoySoN+frmnnadf7SYacyytK+c9axp5z5rGk2a6TdV6fIhHd7Xxk22tPHugGzO4eFkdG85fzNkLIyytC7O4shSzeMtjX3s/+zr6eeG1Yzyxt5PugVHMYPmCCGcviLC0PsxZ9WFKQ37ae4fp6BuhvW+EXS297PFW7CsJ+Dh/SRXnJW+VNFaXURrwE/QbZkYs5jg+NEbXwAid/aN09o/Q2Rd/3DUwwui4wzHx3VBZGqQ2HKI2HKIuHKK+ooT6SAn1kRCRkkDaHzPOOUbGY/QOjdE1MEqXd+zugVEGR6MMjo4zOBplZDxG0GeEAj5CAR+lAT+VZUGqyoJUlcfvK0sDREqCVJQGKA/5M/54cs4xHnOMRx1jsRjjUUfM+44z4muN+M0IBoyg30fAZzP+IeacYyzqiMYcURf/nFjMYRjmA78Zfp8R8MXvC+EHn0JBTplozLGz5TgHOgfY3zHAqx39LF8Y4dPrZm8Z04ThsSj//eIRGqrLuHJF/az/Y9vf0c8/bdnLga5BjhwborN/BICz6sNsOL+BDec3sHIG610c6h7kvhePcO8LhzmY0j8SCvgoCfjoG564Cqo+EuLKFQv47XMWcPnyeuqyaGkdGxjluYPdPHOgm5cP9bCrtZfB0RMv5/X7jNKAj+HxWNpR7vFADyVbGmbgHPQOj51QX6qSgI/SoJ+g10ry+4zB0Sh9w2MntGom8/uM8qCfkqCPsahjdDzGaDR9XZMlvmgDPsNnlgyA8RlMuxLy6g4GfBNBAcm/32TQxBzj0RjjMcdYNDblf9tUn5NoTQZ8PgLen5fPDDPwmSXHzzjvf5xXQ8yBw/H196/h4rPqcv7vBIWCyKwYGo1yfGiMRZUlsxJAzjmO9AzxWtegdxtgaCzKWfVhViyqYPnCCAsr8v+saMxxoHOAnS3H6egbYXgsytBYlKHRGGUhH3XhEuoiIerCJdRXhKiPlFBTHsKfocU1Mh7l2EBK66JvhM7+EboGRhkZizIWc4x5YVNe4qeiNP7LvqI0SJ3XwqiLlFAbDhEu8RPy+9L+N46Ox+gdHuP40Bg9g2P0Do3RNxJfqbBveJzBkXGiziVDIOZc8ss8kLw3gilfuvE/9/jxE1/uY9EYo14YTTyPh4vztnc4/BY/XsDnw++1aBJf6okveZ/3he4zS36JJ1oQ0ahLftZY1PusRMB4LZmYw7uPtzQSK9snwslnE49vvvIsVjbMbCyTQkFERJKyCQVdkioiIkkKBRERSVIoiIhIkkJBRESSFAoiIpKkUBARkSSFgoiIJCkUREQkaV4OXjOzDiD7xYFPVA90zmI5s62Q6yvk2qCw6yvk2qCw6yvk2qCw65tc25nOuQVT7TAvQyEfZrZ1uhF9c6mQ6yvk2qCw6yvk2qCw6yvk2qCw65tJbTp9JCIiSQoFERFJKsZQ2DTXBUyjkOsr5NqgsOsr5NqgsOsr5NqgsOvLubai61MQEZHMirGlICIiGRRVKJjZtWa2x8z2mdltc1zL982s3cx2pLxWa2abzWyvd18zh/Wdbma/NLNmM9tpZp8qlBrNrNTMnjWzl73a/qZQakup0W9mL5rZgwVY20Ez225mL5nZ1gKsr9rM7jGz3d7//y4phPrM7Bzvzyxx6zWzTxdCbSk13ur9m9hhZnd5/1Zyqq9oQsHM/MB3gHcBq4APmtmqOSzp34FrJ712G7DFObcC2OI9nyvjwJ8751YCFwOf8P68CqHGEeBq59waYC1wrZldXCC1JXwKaE55Xki1AVzlnFubcrliIdX3j8DPnXPnAmuI/znOeX3OuT3en9la4LeAQeC+QqgNwMyWAJ8Empxz5wF+4Iac63POFcUNuAR4OOX554HPz3FNS4EdKc/3AA3e4wZgz1z/uaXUdj+wvtBqBMqBF4C3F0ptwGneP76rgQcL7e8WOAjUT3qtIOoDKoEDeP2dhVZfSj3vBJ4spNqAJcAhoBYIAA96deZUX9G0FJj4A0s47L1WSBY551oBvPuFc1wPAGa2FLgAeIYCqdE7PfMS0A5sds4VTG3APwB/AcRSXiuU2iC+FvwjZva8md3svVYo9Z0FdAD/5p1++1czCxdQfQk3AHd5jwuiNufcEeAbwOtAK3DcOfdIrvUVUyikW41cl15Nw8wiwI+BTzvneue6ngTnXNTFm/GnAReZ2XlzXBIAZvZuoN059/xc1zKFy5xzFxI/lfoJM7tyrgtKEQAuBG53zl0ADDD3p9pOYGYh4L3Aj+a6llReX8FGYBnQCITN7CO5HqeYQuEwcHrK89OAljmqJZM2M2sA8O7b57IYMwsSD4T/dM7d671cUDU653qAx4j3zxRCbZcB7zWzg8DdwNVm9oMCqQ0A51yLd99O/Jz4RQVU32HgsNfyA7iHeEgUSn0QD9MXnHNt3vNCqW0dcMA51+GcGwPuBS7Ntb5iCoXngBVmtsxL+huAB+a4pskeAG7yHt9E/Dz+nDAzA74HNDvnvpXy1pzXaGYLzKzae1xG/B/D7kKozTn3eefcac65pcT/P/YL59xHCqE2ADMLm1lF4jHxc847CqU+59xR4JCZneO9dA2wiwKpz/NBJk4dQeHU9jpwsZmVe/9+ryHeSZ9bfXPZWTMHHTEbgFeA/cAX5riWu4if9xsj/uvoj4A64h2Ue7372jms73Lip9e2AS95tw2FUCPwVuBFr7YdwF97r895bZPqfAcTHc0FURvxc/Yve7ediX8HhVKfV8taYKv39/vfQE2h1Ef8woYuoCrltYKozavlb4j/QNoB/F+gJNf6NKJZRESSiun0kYiITEOhICIiSQoFERFJUiiIiEiSQkFERJIUCiIikqRQEBGRJIWCiIgk/X+Zle8GNucE/gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(Evals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "id": "accompanied-memorial",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.0781674703685713"
      ]
     },
     "execution_count": 133,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Evals[-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "north-packet",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}