HF-Roothaan-equations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "79b12ed9",
   "metadata": {},
   "source": [
    "# Solving the many-electron Schrodinger Equation numerically\n",
    "\n",
    "$\\renewcommand{\\vec}[1]{\\mathbf{#1}}$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c482dc7",
   "metadata": {},
   "source": [
    "## Systems\n",
    "\n",
    "- H atom\n",
    "- He ion\n",
    "- He\n",
    "- H2O\n",
    "\n",
    "## Methods\n",
    "\n",
    "HF + DFT\n",
    "\n",
    "## Motion\n",
    "\n",
    "Geometry + MD\n",
    "\n",
    "## Visualisation\n",
    "\n",
    "wavefunctions\n",
    "densities...\n",
    "\n",
    "## something more? \n",
    "\n",
    "periodic?\n",
    "NaCl + water\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d59b8cf6",
   "metadata": {},
   "source": [
    "You have seen the Hartree-Fock equations introduced in Nanophysics. \n",
    "\n",
    "They give a way to solve the Schrodinger equation accounting for the electron-electron interaction within a mean-field approximation. That is, we do not account for where the electrons actually are relative to each other at any instant instead we just allow the electron to interact with the average distribution of the other electrons. \n",
    "\n",
    "Going beyond the mean-field introduces *correlation* into the motion of the different electrons - and get increasingly complicated."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7246a5ac",
   "metadata": {},
   "source": [
    "## Hartree-Fock equations\n",
    "\n",
    "the many-body wavefunction in HF is given by\n",
    "\n",
    "$$\n",
    "\\Psi(\\vec{x}_i, \\vec{x}_2, \\ldots \\vec{x}_N) = \\frac{1}{\\sqrt{N_e!}} \\begin{vmatrix}\n",
    "\\psi_1(\\vec{x}_1) & \\psi_2(\\vec{x}_1) &\\cdots & \\psi_{N/2}(\\vec{x}_1) \\\\\n",
    "\\psi_1(\\vec{x}_2) & \\psi_2(\\vec{x}_2) &\\cdots & \\psi_{N/2}(\\vec{x}_2) \\\\\n",
    "\\vdots & \\vdots & \\vdots & \\vdots \\\\\n",
    "\\psi_1(\\vec{x}_N) & \\psi_2(\\vec{x}_N) &\\cdots & \\psi_{N/2}(\\vec{x}_N) \\\\\n",
    "\\end{vmatrix}\n",
    "$$\n",
    "\n",
    "We plug this *ansatz* into the many-body Hamiltonian and chug through a bunch of tedious calculations (see notes), introduce a basis set and and eventually it emerges that\n",
    "\n",
    "$\\vec{FC} = \\vec{SCE}$\n",
    "\n",
    "i.e. a set of matrix equations called the Roothaan-Hall equations."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86ed54c4",
   "metadata": {},
   "source": [
    "## Roothaan-Hall equations\n",
    "\n",
    "the Hartree-Fock equations are a set of non-linear integro-differential equations to be solved subject to a set of conditions (correct number of electrons and orthogonal wave-functions). This is essentially impossible to treat analytically.\n",
    "\n",
    "We can make progress by introducing *basis functions* into the system of equations and, using a trick we used to develop a tight-binding model in Condensed Matter Physics, move from differential equations to matrix equations.\n",
    "\n",
    "To do this we will need to introduce *matrix elements* of the basis functions with respect to each other and of the Hartree-Fock functional."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf09e85c",
   "metadata": {},
   "source": [
    "### Approximations\n",
    "\n",
    "Fundamental approximations that we are going to introduce:\n",
    "\n",
    "- Born-Oppenheimer approximation separating motion of nuclei and electrons;\n",
    "- Independent electron approximation: write our many-body wavefunction as *single determinant* built up of single particle wavefunctions;\n",
    "- Write our single particle functions as Linear Combination of Atomic Orbitals (basis-functions)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b15c3fe",
   "metadata": {},
   "source": [
    "## Solving the Roothaan-Hall equations using PYSCF\n",
    "\n",
    "PYSCF is a python library for solving the many-body Schrodinger equation. It has many methods but includes Hartree-Fock method as a starting point for many of them.\n",
    "\n",
    "The website fro PySCF is https://pyscf.org/index.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "7883c783",
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can access submodules from pyscf like this\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pyscf"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "882addb7",
   "metadata": {},
   "source": [
    "As a super quick example let's run a HF model of a water molecule"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "7fa88dc3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# here we create a H2O molecule and run a HF calculation on it using a minimal basis set (STO3G)\n",
    "mol_h2o = pyscf.gto.M(atom = 'O 0 0 0; H 0 1 0; H 0 0 1', basis = 'cc-pVDZ')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bec552d7",
   "metadata": {},
   "source": [
    "`mol_h2o` now is our atomistic model of the water molecule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e6177a24",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "converged SCF energy = -76.0167894720692\n"
     ]
    }
   ],
   "source": [
    "rhf_h2o = pyscf.scf.RHF(mol_h2o)\n",
    "e_h2o_DZ = rhf_h2o.kernel()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31ea14da",
   "metadata": {},
   "source": [
    "### Converging basis sets\n",
    "\n",
    "In principle, the basis sets used are approximations and we should be able to estimate how well we are doing - we can increase the quality of the basis sets by using more gaussian functions allowing more flexibility in the variational solution of the HF equations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "50dda972",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "converged SCF energy = -76.0167894720692\n",
      "converged SCF energy = -76.0456257970043\n",
      "converged SCF energy = -76.0530065175422\n",
      "converged SCF energy = -76.0551821620884\n"
     ]
    }
   ],
   "source": [
    "my_basis_sets = ['cc-pVDZ', 'cc-pVTZ', 'cc-pVQZ', 'cc-pV5Z']\n",
    "energy = []\n",
    "for basis in my_basis_sets:\n",
    "    mol_h2o = pyscf.gto.M(atom = 'O 0 0 0; H 0 1 0; H 0 0 1', basis = basis)\n",
    "    rhf_h2o = pyscf.scf.RHF(mol_h2o)\n",
    "    energy.append(rhf_h2o.kernel())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a1a9f8c6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaIAAAEKCAYAAABQRFHsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA1E0lEQVR4nO3deXxV1bn/8c83EwkzSJiNqMggEBBStdYBFdS2KohItWrtvf3V6m3v7a11rLXXtrbFsb+29rZS7622/rQO4FSwVKlDq1ZlCoOAoKImhHkewhCe3x97BQ7xJDnASfZJ8rxfr/06Zw9r72fneHxYa6+zlswM55xzLi5ZcQfgnHOuZfNE5JxzLlaeiJxzzsXKE5FzzrlYeSJyzjkXK09EzjnnYhVLIpL0uKS5YVkuaW7CvmJJb0paKGm+pPwk5TtLelHS0vDaKWw/QtLLkrZKur9GmVckLUm4btcGv1HnnHP1yonjomb2per3ku4FNoX3OcAjwJVmVirpCGB3klPcDMwws4mSbg7rNwGVwG3A4LDUdLmZzUzrzTjnnDsssTbNSRIwAXgsbDoHmGdmpQBmts7MqpIUHQM8HN4/DIwNx28zs38QJSTnnHNNQCw1ogSnAavMbGlY7weYpOlAIfAnM7srSbluZlYBYGYVB9HM9ntJVcBk4A5LYViJLl26WJ8+fVI8vXPOOYBZs2atNbPCVI5tsEQk6SWge5Jdt5rZs+H9ZeyvDVXHcyrwGWA7MEPSLDObkYaQLjezckntiBLRlcAfaon9auBqgKKiImbO9NY855w7GJI+SvXYBktEZjaqrv3hedA4YETC5jLgVTNbG46ZBgwHaiaiVZJ6hNpQD2B1CvGUh9ctkh4FTqSWRGRmk4BJACUlJT4Yn3PONaA4nxGNAhabWVnCtulAsaTWIVGdAbybpOxzwFXh/VXAs0mO2UdSjqQu4X0ucD6w4DDjd845lwZxPiO6lAOb5TCzDZLuA94BDJhmZlMBJD0I/Db0epsIPCHpa8DHwCXV55C0HGgP5EkaS9QB4iNgekhC2cBLwO8a9O6cc86lRD4NRN1KSkrMnxE559zBCc/3S1I51kdWcM45F6u4u283S8/MKefu6UtYsXEHPTsWcMO5/Rl7Qq+4w3LOuYzkiSjNnplTzi1T5rNjd/Q73PKNO7hlynwAT0bOOZeEN82l2d3Tl+xLQtV27K7i7ulLYorIOecymyeiNFuxccdBbXfOuZbOE1Ga9exYcFDbnXOupfNElGY3nNufgtzsA7blZosbzu0fU0TOOZfZvLNCmlV3SLh7+hLKN+4gLzuLrCw48ejOMUfmnHOZyX/QWo/D/UHrx+u2c94vXmN4USf+8K8nkpWlNEbnnHOZyX/QmkGKjmjNrV8cyD+WreWRt1IejNY551oMT0SN4MsnFnFGv0J+Nm0xH67dFnc4zjmXUTwRNQJJ3HlxMbnZ4vonS6na682hzjlXzRNRI+neIZ8fjx3MrI82MOm1D+IOxznnMoYnokZ04dCefGFId37+4nssXrk57nCccy4jeCJqRJK4Y+wQ2hfk8p3HS9m1Z2/cITnnXOw8ETWyzm3ymDhuCIsqNvPLGUvjDsc552LniSgGo47vxiUjevPfryxjzscb4g7HOediFUsikvS4pLlhWS5pbsK+YklvSlooab6k/CTlO0t6UdLS8NopbB8taVYoN0vSWQllRoTtyyT9UlKsvyy97YLj6dGhgO8+UcqOXVX1F3DOuWYqlkRkZl8ys2FmNgyYDEwBkJQDPAJcY2aDgJHA7iSnuBmYYWbHATPCOsBa4AIzGwJcBfwxocxvgKuB48JyXppv66C0z8/l7vHFfLB2G3f+ZXGcoTjnXKxibZoLtZIJwGNh0znAPDMrBTCzdWaWrLowBng4vH8YGBuOn2NmK8L2hUC+pFaSegDtzexNi8Y0+kN1mTid0rcLXz2lDw+9sZw3lq2NOxznnItF3M+ITgNWmVn1U/t+gEmaLmm2pBtrKdfNzCoAwmvXJMdcDMwxs51AL6AsYV9Z2Ba7m84bwDFd2nDDU/PYXJms8uecc81bgyUiSS9JWpBkGZNw2GXsrw1BNBr4qcDl4fUiSWcfwrUHAXcC36jelOSwWoc3kHS1pJmSZq5Zs+ZgL39QCvKyuWfCUCo27eDHz7/boNdyzrlM1GDTQJjZqLr2h+dB44ARCZvLgFfNbG04ZhownOg5UKJVknqYWUVodludcN7ewNPAV8zs/YTz9k4o3xtYQS3MbBIwCaLRt+u6j3QYXtSJa0cey69ffp9zB3Vn1PHdGvqSzjmXMeJsmhsFLDazxCaz6UCxpNYhUZ0BJKsmPEfUGYHw+iyApI7AVOAWM3u9+uDQfLdF0snhudRXqstkim+f3Y+BPdpz85T5rN+2K+5wnHOu0cSZiC7lwGY5zGwDcB/wDjAXmG1mUwEkPSipem6LicBoSUuB0WEd4FtAX+C2hO7h1c+PrgUeBJYB7wMvNNSNHYq8nCzumzCUTTt28f1n5uPzRDnnWgqfGK8ehzsx3sH671eWcddflvCLS4cxZlhG9KdwzrmD5hPjNWHfOP1Yhhd15LZnFrByU2Xc4TjnXIPzRJRhsrPEvROGsbvKuGnyPG+ic841e56IMtDRXdpwyxcG8Op7a3js7U/iDsc55xqUJ6IMdcVJR3Fq3y7cMfVdPl63Pe5wnHOuwXgiylBZWeKu8cVkZ/n04s655s0TUQbr2bGA2y8YxNvL1/O///gw7nCcc65BeCLKcOOG92L08d24+69LeG/VlrjDcc65tPNElOEk8bNxQ2jbKofrnpjL7iqfXtw517x4ImoCurRtxU8vGsyC8s3c/7dlcYfjnHNp5YmoiThvcA8uOqEX97+8jHllG+MOxznn0sYTURNy+4WDKGzbiuueKKVyt08v7pxrHjwRNSEdCnK5a3wxy1Zv5Z7pS+IOxznn0sITURNzer9Crji5iP95/UP++cG6uMNxzrnD5omoCfreFwZS1Lk11z9Zytade+IOxznnDosnoiaodV4O914ylPKNO/jJVJ9e3DnXtHkiaqJK+nTm6tOP4bG3P+HlxavrL+CccxnKE1ETdt3ofvTv1o6bJs9j43afXtw51zTFkogkPZ4wlfdySXMT9hVLelPSQknzJeUnKd9Z0ouSlobXTmH7aEmzQrlZks5KKPOKpCVJphBvslrlZHPvhKGs37aL255dGHc4zjl3SGJJRGb2JTMbZmbDgMnAFABJOcAjwDVmNggYCexOcoqbgRlmdhwwI6wDrAUuMLMhwFXAH2uUu7z6umbWLNqzBvfqwLfPPo7nS1fw53kr4g7HOecOWqxNc5IETAAeC5vOAeaZWSmAma0zs2S/3BwDPBzePwyMDcfPMbPq/xsvBPIltWqg8DPGtSOPZeiR0fTiq7f49OLOuaYl7mdEpwGrzGxpWO8HmKTpkmZLurGWct3MrAIgvCZrZrsYmGNmOxO2/T40y90WkmCzkJOdxb2XDGX7ripumTzfpxd3zjUpDZaIJL0kaUGSZUzCYZexvzYEkAOcClweXi+SdPYhXHsQcCfwjYTNl4cmu9PCcmUd5a+WNFPSzDVr1hzs5WPRt2tbbjxvADMWr+bJmWVxh+OccylrsERkZqPMbHCS5VnY9zxoHPB4QrEy4FUzW2tm24FpwPAkp18lqUc4Tw9g3/MeSb2Bp4GvmNn7CfGUh9ctwKPAiXXEPsnMSsyspLCw8ND+ADH4l1P6cPIxnfnRn9/lk/U+vbhzrmmIs2luFLDYzBL/+T4dKJbUOiSqM4Bkv9h8jqgzAuG1Orl1BKYCt5jZ69UHS8qR1CW8zwXOBxak93bil5Ul7h4/FIAbniplr08v7pxrAuJMRJdyYLMcZrYBuA94B5gLzDazqQCSHpRUEg6dCIyWtBQYHdYBvgX0BW6r0U27FTBd0rxw3nLgdw14b7E5snNrbjt/IP/8YD0PvbE87nCcc65e8gfbdSspKbGZM2fGHcZBMTO+9vBMXl+2lqn/cRp9u7aNOyTnXAsjaZaZldR/ZPy95lwDkMTEcUMoyMvmu0+WssenF3fOZTBPRM1U1/b53DF2MKWfbOQ3r7xffwHnnIuJJ6Jm7PzinlwwtCe/mLGUBeWb4g7HOeeS8kTUzP14zCA6t8nju0+UsnOPTy/unMs8noiauY6t87jz4mKWrNrCfS++F3c4zjn3KZ6IWoAzB3TlshOPZNJrHzBz+fq4w3HOuQN4Imohbv3i8fTuVMB3nyxlm08v7pzLIJ6IWoi2rXK4Z/xQPl6/nZ+9sCjucJxzbh9PRC3IScccwdc+dzSP/PNjXnuvaQzm6pxr/lJORJLaSMpuyGBcw7v+3P7RSN1PzWPTjmRzDjrnXOOqNRFJypL0ZUlTJa0GFgMVYQrvuyUd13hhunTJz83mvglDWbN1Jz98zqcXd87Fr64a0cvAscAtQHczO9LMuhLN5fNPYKKkKxohRpdmxb078s0z+zJlTjl/WbAy7nCccy1cTh37RpnZp9puzGw9MBmYHKZUcE3Qv5/Vl78tXsWtT8+npE8nurRt9jOqO+cyVF01onaSOicsnWpOr50sUbmmITc7i/smDGPLzj18b4pPL+6ci09diWgWMDO8zgJmA6vDFOB9GiE218D6dWvH9ef046/vrmLK7PK4w3HOtVC1JiIzO9rMjgmv1Ush8N/AbxsvRNeQvnbqMZzYpzO3P7eQFRt3xB2Oc64FOujfEZnZFKDr4VxU0uMJM6gulzQ3YV+xpDdD77z5kvKTlO8s6UVJS8Nrp7D9xITzlkq6KKHMiHC+ZZJ+WbOZsaXKzhL3XDKUKjNufGqeTy/unGt0B52IJLU9lHKJzOxLZjbMzIYRdXyYEs6dAzwCXGNmg4CRQLLnUDcDM8zsOGBGWAdYAJSE854HPBDOCfAb4GrguLCcdzj30JwUHdGaW784kH8sW8sjb30UdzjOuRam1l5zkq5LsrkTcCFwfzouHmolE4CzwqZzgHlmVgpgZutqKTqGKEkBPAy8AtxkZtsTjskHLFynB9DezN4M638AxgIvpOM+moMvn1jEXxeu4qfTFnHacYUc3aVN3CE551qIOnvN1VjaAiuBK8zsd2m6/mnAKjNbGtb7ASZpuqTZkm6spVw3M6sACK/7mgolnSRpITCfqGa1B+gFlCWULwvbXCCJOy8uJi87i+8+MZcqb6JzzjWSWmtEZvbDZNsl5Uu6xMyerOvEkl4CuifZdauZPRveXwY8ViOeU4HPANuBGZJmmdmMuq5VI+63gEGSBgIPS3oBSPY8qNb/00q6mqgZj6KiolQv3eR175DPj8cO5tt/mssDr73Pv43sG3dIzrkWIKVnPZKyJX0+NGl9BHypvjJmNsrMBidZng3nzAHGAY8nFCsDXjWztaGZbRowPMnpV4Xmtupmt9VJrr8I2AYMDuftnbC7N7CijtgnmVmJmZUUFhbWd6vNyoVDe/KFId35+Yvvsahic9zhOOdagDoTkaTTJf0WWA78H6JnOEeb2fg0XHsUsNjMEpvMpgPFklqHRHUG8G6Sss8BV4X3VwHVye3o6s4Jko4C+gPLQ/PdFkknh+dSX6ku4w4kiTvGDqFDQR7XPVHKrj174w7JOdfM1TXoaRkwEXgdON7MLgZ21OgQcDgu5cBmOcxsA3Af8A4wF5htZlNDPA9KKgmHTgRGS1oKjA7rEDXrlYbu4E8D/2Zma8O+a4EHgWXA+3hHhVp1bpPHxHFDWFSxmV/M8OnFnXMNS7UN7SLpF0Q9y+YDjxLVIOab2TGNFl0GKCkpsZkzZ8YdRixueLKUybPLeOraUxhe1CnucJxzTUh4vl9S/5F1j6zwbaAPUQ3lTOA9oFDShPBbItfM/eCC4+nRoYDrnyhlx66quMNxzjVTdT4jssjfzOzrREnpy0S1pOUNHpmLXbv8XO4eX8wHa7dx518Wxx2Oc66ZSnmEBDPbbWbPm9mXgSMbMCaXQU7p24WvntKHh95YzhvL1tZfwDnnDlJdnRWel3RBLXMO9ZD0I0n/2oCxuQxx03kDOKZLG254ah6bK33mD+dcetVVI/o60cgHiyW9I2mapL9J+hB4AJhlZv/bKFG6WBXkZXPPhKFUbNrBj59P1pveOecOXV0jK6wEbgRuDPMP9QB2AO+lsQu3ayKGF3Xi2pHH8uuX3+fcQd0ZdXy3uENyzjUTKT0jMrPlZvammc31JNRyffvsfgzs0Z6bp8xn/bZdcYfjnGsmDms6B9ey5OVkcd+EoWzasYvvP+PTizvn0sMTkTsoA3u05zuj+zFt/kqeK611uD7nnEtZvYlI0vmSPGG5fb5x+rEML+rIbc8sYOWmyrjDcc41cakkmEuBpZLuClMruBYuO0vcO2EYu6uMmybP8yY659xhqTcRmdkVwAlEA4X+XtKbkq6W1K7Bo3MZ6+gubbjlCwN49b01PPr2x3GH45xrwlLtNbcZmAz8iagb90XAbEn/3oCxuQx3xUlHcWrfLvxk6iI+Wrct7nCcc01UKs+ILpD0NPA3IBc40cw+DwwFrm/g+FwGy8oSd40vJjtLXP9kqU8v7pw7JKnUiC4Bfm5mxWZ2t5mtBgi/J/Ihflq4nh0LuP2CQbyzfAP/848P4g7HOdcEpfKM6Ctm9lot+2akPyTX1Iwb3otzju/GPdPf471VW+IOxznXxKTSNLdF0uYayyeSnpbUoibJc8lJ4qfjhtA2P4frnpjL7iqfXtw5l7pUmubuA24AegG9iZ4L/Y6o48IhDXoq6XFJc8OyPEztXb2vOPTMWyhpvqT8JOU7S3pR0tLw2ilsPzHhvKWSLkoo84qkJQn7ux5K7C65Lm1b8dOLBrOgfDP3/21Z3OE455qQVBLReWb2gJltMbPNZjYJ+IKZPQ4c0vzRZvYlMxtmZsOIeuNNAZCUAzwCXGNmg4CRQLJ5B24GZpjZccCMsA6wACgJ5z0PeCCcs9rl1detftbl0ue8wT246IRe3P/yMuaVbYw7HOdcE5FKItobpgfPCsuEhH2H1U1KkoAJwGNh0znAPDMrBTCzdWaWbI7qMcDD4f3DRLPGYmbbzWxP2J5/uPG5g3f7hYMobNuK654opXK3Ty/unKtfKonocuBKYDWwKry/QlIB8K3DvP5pwCozWxrW+wEmabqk2ZJurKVcNzOrAAiv+5rZJJ0kaSEwn6hmtSeh3O9Ds9xtIQm6NOtQkMtd44tZtnor90xfEnc4zrkmoNb5iAAkZQPXmtkFtRzyjzrKvgR0T7LrVjN7Nry/jP21oep4TgU+A2wHZkiadTC988zsLWBQGI7oYUkvmFklUbNceRgRYjJRQv1DLbFfDVwNUFRUlOqlXXB6v0KuOLmI/3n9Q0Yd342Tjzki7pCccxmszhpRaBYbcSgnNrNRZjY4yfIs7HseNA54PKFYGfCqma0Nv1OaBgxPcvpVknqE8/Qgqq3VvP4iYBswOKyXh9ctwKPAiXXEPsnMSsyspLCw8OBv3vG9LwykqHNrrn+ylK0799RfwDnXYqXSNDdH0nOSrpQ0rnpJw7VHAYvNrCxh23SgWFLrkKjOAJLNTf0ccFV4fxVQndyOru6cIOkooD+wXFKOpC5hey5wPlHHBtdAWuflcO8lQynfuIOfTPXpxZ1ztUslEXUG1gFnAReE5fw0XPtSDmyWw8w2EHUXfweYC8w2s6kAkh6UVBIOnQiMlrQUGB3WIWrWKw3dwZ8G/s3M1gKtgOmS5oXzlhN1QXcNqKRPZ64+/Rgee/sTXl7snRSdc8nJh/CvW0lJic2cOTPuMJqsnXuquPBXr7Nh+y7++p3T6dg6L+6QnHONIDzfL6n/yNRGVugnaYakBWG9WNL3DzdI1zK0ysnm3glDWb9tF7c9uzDucJxzGSiVprnfAbcQflhqZvOImtWcS8ngXh349tnH8XzpCv48z6cXd84dKJVE1NrM3q6xzbtBuYNy7chjGXpkR77/zAJWb/bpxZ1z+6WSiNZKOpYwSoGk8UBFg0blmp2c7CzuvWQoO3ZVcfOU+T69uHNun1QS0TeBB4ABksqB/wSubcigXPPUt2tbbjpvAH9bvJonZn4SdzjOuQyRynxEH5jZKKAQGGBmp5rZ8gaPzDVLXz2lDycf05kfPf8un6zfHnc4zrkMkEqvuVaSvgx8G/iOpB9I+kHDh+aao6wscff4oUjR9OJ7fXpx51q8VJrmniUa7XoP0ZA51Ytzh+TIzq257fyBvPXhen7/xvK4w3HOxazOQU+D3mZ2XoNH4lqUCSVHMn3hKu76y2LO6FdI365t4w7JOReTVGpEb0ga0uCRuBZFEhPHDaEgL5vvPlnKHp9e3LkWK5VEdCowK0yzPS9M3z2voQNzzV/X9vncMXYwpZ9s5DevvB93OM65mKTSNPf5Bo/CtVjnF/dk+sJV/GLGUs4c0JXBvTrEHZJzrpGl0n37I+BI4Kzwfnsq5ZxL1Y/HDKJzmzy++0QpO/f49OLOtTSpdN/+L+AmovHmAHKBRxoyKNeydGydx50XF7Nk1Rbue/G9uMNxzjWyVGo2FwEXErpsm9kKoF1DBuVanjMHdOWyE49k0msfMHP5+rjDcc41olQS0S6LBgarHmuuTcOG5FqqW794PL07FfDdJ0vZ5tOLO9dipJKInpD0ANBR0teBl/DZTV0DaNsqh3vGD+Xj9dv52QuL4g7HOddIUumscA/wFDAZ6A/8wMx+dTgXlfS4pLlhWR6m9q7eVyzpTUkLQ1fx/CTlO0t6UdLS8Nqpxv4iSVslXZ+wbUQ43zJJv5Skw7kH1zBOOuYIvva5o3nknx/z2ntr4g7HOdcIUur9ZmYvmtkNZna9mb14uBc1sy+Z2TAzG0aU4KYASMoh6ghxjZkNAkYSJuSr4WZghpkdB8wI64l+DrxQY9tvgKuB48Lio0VkqOvP7U/frm258al5bNqe7ON3zjUnsXbDDrWSCcBjYdM5wDwzKwUws3Vmlqw/7xjg4fD+YWBswjnHAh8ACxO29QDam9mb4XnXHxLLuMySn5vNfROGsmbrTm5/3qcXd665i/v3QKcBq8xsaVjvB5ik6ZJmS7qxlnLdzKwCILx2hX0dKW4Cfljj+F5AWcJ6WdjmMlRx745868y+PD2nnL8s8HkYnWvOGiwRSXpJ0oIky5iEwy5jf20IopEeTgUuD68XSTr7IC77Q+DnZra1ZjhJjq11/gFJV0uaKWnmmjX+nCIu3zqrL4N7ted7Ty9gzZadcYfjnGsgtQ7xI2k+yf9nLcDMrLiuE4fJ9GoVngeNA0YkbC4DXjWzteGYacBwoudAiVZJ6mFmFaHZbXXYfhIwXtJdQEdgr6RKoudQvRPK9wZW1BH7JGASQElJiU+YE5Pc7CzumzCM83/1D7739HwmXTkC72PiXPNT11hz5zfwtUcBi80ssclsOnCjpNbALuAMoo4HNT0HXAVMDK/PApjZadUHSLod2Gpm94f1LZJOBt4CvgIcVs8/1zj6dWvH9ef046fTFjN5djnjR/Suv5BzrkmpNRGFceUa0qUc2CyHmW2QdB/wDlFtbJqZTQWQ9CDwWzObSZSAnpD0NeBj4JIUrnct8BBQQNSjrmavOpehvnbqMbz07mp++NxCPnvsEfTqWBB3SM65NFLUiayOA6JaxK+AgUAekA1sM7P2DR9e/EpKSmzmzJlxh9HifbxuO+f94jWGF3XiD/96IllZ3kTnXCaTNMvMSlI5NpXOCvcTdSpYSlSb+D94s5ZrZEVHtObWLw7kH8vW8shbDV1Zd841plR/0LoMyDazKjP7PXBmw4bl3Kd9+cQizuhXyE+nLeLDtdviDsc5lyapJKLtkvKAuZLukvQdwAc+dY1OEndeXExedhbffWIuVXu9Q6NzzUEqiejKcNy3iKaCOJKo27Vzja57h3x+PHYwsz/eyAOv+fTizjUHqSSisWZWaWabzeyHZnYdDd+127laXTi0J18Y0p2fv/geiyo2xx2Oc+4wpZKIrkqy7atpjsO5lEnijrFD6FCQx3VPlLJrz964Q3LOHYZaE5GkyyQ9Dxwt6bmE5RVgXaNF6FwSndvkMXHcEBZVbOYXM3x6ceeasrpGVngDqAC6APcmbN8CzGvIoJxLxajju3HJiN785pX3OXtgN4YXdaq/kHMu49RaIzKzj8zsFTP7LLAYaBeWMjPzeZxdRvjBBcfTo0MB1z9Ryo5dyWYMcc5lunqfEUm6BHibaBidCcBbksY3dGDOpaJdfi53jy/mg7XbuPMvi+MOxzl3COpqmqv2feAzZrYaQFIh8BLR9OHOxe6Uvl346il9eOiN5Yw+vhuf69sl7pCccwchlV5zWdVJKFiXYjnnGs1N5w3gmC5tuOHJUjZX+vTizjUlqSSUv4QZU78q6avAVHzkapdhCvKyuWfCUFZuruRHz78bdzjOuYNQbyIysxuAB4BiYCgwycxqm8LbudgML+rEtSOP5alZZbz47qq4w3HOpSiVzgp3mtkUM7vOzL5jZk9LurMxgnPuYH377H4M7NGeW6bMY91Wn17cuaYglaa50Um2fT7dgTiXDnk5Wdw3YSibduzm+88soL75tpxz8atrZIVrJc0H+kual7B8iP+g1WWwgT3a853R/XhhwUqeK10RdzjOuXrUVSN6FLgAeC68Vi8jzOyKw7mopMclzQ3LcklzE/YVS3pT0kJJ8yXlJynfWdKLkpaG10419hdJ2irp+oRtr0haknDdrodzDy6zfeP0Yxle1JHbnlnAyk2VcYfjnKtDXSMrbDKz5WZ2WRhloXpZf7gXNbMvmdkwMxsGTAamAEjKAR4BrjGzQcBIIFlf3JuBGWZ2HDAjrCf6Ocl79l1efd0aXdJdM5OdJe6dMIzdVcZNk+d5E51zGSzW3wNJEtFoDY+FTecA88ysFMDM1plZsnFbxgAPh/cPA2MTzjkW+ABY2DBRu6bi6C5tuOULA3j1vTU8+vbHcYfjnKtF3D9MPQ1YZWZLw3o/wMLvlmZLqq2beDczqwAIr10BJLUBbgJ+WEu534dmudtCEkxK0tWSZkqauWbNmkO5L5chrjjpKE7t24WfTF3ER+t8enHnMlGDJSJJL0lakGQZk3DYZeyvDUE05NCpwOXh9SJJZx/EZX8I/NzMtibZd7mZDSFKfqcRzTyblJlNMrMSMyspLCw8iMu7TJOVJe4aX8xeM0bf9xp9bp7K5yb+jWfmlMcdmnMuSGWsuUNiZqPq2h+eB40DRiRsLgNeNbO14ZhpwHCi50CJVknqYWYVknoA1c97TgLGS7oL6AjslVRpZvebWXmIa4ukR4ETgT8c1k26JuHtD9dTtdfYXRU9JyrfuINbpswHYOwJveIMzTlHvE1zo4DFZlaWsG06UCypdUhUZwDJxmt5jv0zx14FPAtgZqeZWR8z6wP8X+CnZna/pBxJXQAk5RJNdb6gAe7JZaC7py/Zl4Sq7dhdxV0+WrdzGSHORHQpBzbLYWYbgPuAd4C5wGwzmwog6UFJJeHQicBoSUuJfnA7sZ5rtQKmS5oXzlsO/C49t+Ey3YqNO5Jv31TJT6ctYsnKLY0ckXMukbxba91KSkps5syZcYfhDsPnJv6N8iTJKD83iz1Vxp69RnHvDlw8vDcXDu1JpzZ5MUTpXPMiaZaZldR/ZPy95pxrcDec25+C3OwDthXkZjNxXDFvfe9sfnD+8eypMv7ruYWc+NOXuPaRWcxYtIo9VXtjiti5lsVrRPXwGlHz8Myccu6evoQVG3fQs2MBN5zb/1MdFd5dsZnJs8t4Zk4567btokvbVlx0Qk8uHtGbAd3bxxS5c03TwdSIPBHVwxNRy7O7ai+vLFnDU7M+Ycai1ezZawzp1YHxI7zpzrlUeSJKI09ELdu6rTt5rnQFT80qY+GKzeRmi7MHdGP8iN6c0b+Q3Gxv3XYuGU9EaeSJyFX7dNNdHmOH9WJ8iTfdOVeTJ6I08kTkatpdtZdXl6zhqVllzFi8it1VxuBe7Rk/vDcXDutFZ2+6c84TUTp5InJ1Wb9tF8/NLeep2WUsKI+a7s4a0JXxI45kpDfduRbME1EaeSJyqVpUsZnJs8p4Zm45a7fu4og2eYw9oRfjR/RmYA9vunMtiyeiNPJE5A7W7qq9vPZe1HT30qKo6W5Qz/aMH9GbMd5051oIT0Rp5InIHY4N23bt63U3v3yTN925FsMTURp5InLpsnhl1HT39JwVrN26kyPa5DFmWNR0d3xPb7pzzYsnojTyROTSLVnT3fE9qpvuenJE21Zxh+jcYfNElEaeiFxD2rBtF8/Pi5ru5pVtIieruumuN2cO6OpNd67J8kSURp6IXGNZsnILk2eXMWV2+b6muwuH9WT8iN4M6tkh7vCcOyieiNLIE5FrbHuq9vLa0tB09+5qdlXtZWBC010Xb7pzTYAnojTyROTiVN10N3lWGaWh6e7M6qa7/l3Jy/GmO5eZMj4RSXoc6B9WOwIbzWxY2FcMPAC0B/YCnzGzyhrlOwOPA32A5cAEM9sgqQ+wCFgSDv2nmV0TyowAHgIKgGnAty2Fm/dE5DLFe6u2MHlWGVPmlLNmy046t8ljjDfduQyV8YnogACke4FNZvYjSTnAbOBKMyuVdARRkqqqUeYuYL2ZTZR0M9DJzG4KiejPZjY4yXXeBr4N/JMoEf3SzF6oLz5PRC7T7Knay9+XruWpWWW8+O4qdlXtZUD3dowf0ZuxJ/TypjuXEZpMIpIk4GPgLDNbKukLwJfN7Ip6yi0BRppZhaQewCtm1r+2RBSOednMBoT1y0L5b9QXoycil8k2bt/F8+EHs9VNdyP7R013Zw3wpjsXn4NJRDkNHUw9TgNWmdnSsN4PMEnTgULgT2Z2V5Jy3cysAiAko64J+46WNAfYDHzfzP4O9ALKEo4pC9uca9I6ts7jys/24crP9mHpqi08FXrdvbRoFZ1a5+77weygnu2J/t3nXOZpsEQk6SWge5Jdt5rZs+H9ZcBjNeI5FfgMsB2YEbLqjBQvWwEUmdm68EzoGUmDgGTfwFqrgpKuBq4GKCoqSvHSzsXruG7tuOXzA7nhnP78fVnUdPfoWx/z0BvL9zXdjRnWi8J23nTnMkuDJSIzG1XX/vA8aBwwImFzGfCqma0Nx0wDhgM1E9EqST0SmuZWh2vuBHaG97MkvU9UyyoDeieU7w2sqCP2ScAkiJrm6rlV5zJKTnYWZ/bvypn9u0ZNd/MqeGpWGXdMXcTPXljMmf0LQ9NdN2+6cxkhzqa5UcBiM0tsMpsO3CipNbALOAP4eZKyzwFXARPD67MAkgqJOjFUSToGOA74wMzWS9oi6WTgLeArwK8a6L6cyxgdW+dx5clHceXJR+1runt6djkvLVrtTXcuY8TWWUHSQ0Tdq39bY/sVwC1ETWfTzOzGsP1B4LdmNjP0pnsCKCLq7HBJSDYXAz8C9gBVwH+Z2fOhfAn7u2+/APy7d992LdGeqr38IzTd/fXdVezas9eb7lzaNZlec02BJyLXnG3avnvfWHdzP9lIdpYY2S803Q3sSquc7LhDdE2UJ6I08kTkWoplq7fw1Kxynp5TxqrNO+nYOpcxQ3syfsSRDO7lTXfu4HgiSiNPRK6lqdpr+5rupi9cya49e+nfLTTdndCTru3y4w7RNQGeiNLIE5FryTZt382f50dNd3M+3t90d/GI3pztTXeuDp6I0sgTkXORZau3hmkqoqa7DgW5+8a6G9KrgzfduQN4IkojT0TOHShZ012/bm2jse6G9aJre2+6c56I0soTkXO127RjN38O01TMDk13Z4Red95017J5IkojT0TOpeb9NVujaSpml7NycyUdCnK5cGjUdFfc25vuWhpPRGnkici5g1O113g9oelu5569HNc1arq76ARvumspPBGlkSci5w7d5srdTA1j3c36aANZIjTdHcnZA7vylwUruXv6ElZs3EHPjgXccG5/xp7gA+M3B56I0sgTkXPp8f6arUwJ01RUbKqkdV42O/fspWrv/v8HFeRm87NxQzwZNQOeiNLIE5Fz6VW113jj/bVc/YdZ7Nhd9an9BbnZ/Mvn+tCjQz7dOxTQo0M+PTrk07lNnj9nakKa0sR4zrkWJjtLnHZcIZVJkhDAjt1VTHrtA/bsPfAfyXk5WXRvn0/3Dvn0TEhS0XoB3Tvkc0SbPLKyPFk1NZ6InHOx6NmxgPKNOz61vVfHAv5+45ms3baTio2VVGyqZOWmHVRsrqRiYyUrN1Uy6+MNrNxUwe6qA5NVbrbo1n5/YqquTSXWrLq0beXJKsN4InLOxeKGc/tzy5T5BzTPFeRmc8O5/cnKEl3b5dO1XT5Dj0xefu9eY922XazcVEnFph1UbEpIWpsqmfvJRv6yoJJdVXsPKJeTFSWrfbWpjgV0r7HepW0rsj1ZNRpPRM65WFR3SDjUXnNZWaKwXSsK27ViSO8OSY8xM9Zv2/WpJFURkteC8k28+O4qdu45MFllZ4lu7VrRo2OoWbVPSFqhZlXYthU52T7DbTp4Z4V6eGcF55o3M2PD9t1UbNoRalf7a1grw7Ji0w4qd386WXVt12pfYurevoCeHfP3r3cooFu7lpusvLOCc86lSBKd2+TRuU0eg3rWXrPatGP3p5JU9frilVt4efGaT/UCzBIUtmtF9w4FoYPF/iRVvd6tfT65LTRZVYslEUl6HOgfVjsCG81sWNhXDDwAtAf2Ap8xs8oa5TsDjwN9gOXABDPbIKkPsAhYEg79p5ldE8q8AvQAqp+OnmNmq9N+c865ZkcSHVvn0bF1HgN7tE96jJmxuXLPgYlqY3i/uZKlq7fy2ntr2Larqsa5oUvbVgmJKrGjRdTJolv7fPJymm+yiiURmdmXqt9LuhfYFN7nAI8AV5pZqaQjgN1JTnEzMMPMJkq6OazfFPa9X53UkrjczLydzTmXdpLoUJBLh4JcBnSvPVlt2blnf21qY0LS2lzJB2u28caydWzZuedTZbu0bZXQXf3T3de7dWjVZAeZjbVpTtGv0yYAZ4VN5wDzzKwUwMzW1VJ0DDAyvH8YeIX9icg55zKSJNrn59I+P5d+3drVetyWyt2s2lzJio2Vn3pu9fG67bz1wTo2V346WR3RJm9frapHQlNg4np+bt3J6pk55Y0+7FLcz4hOA1aZ2dKw3g8wSdOBQuBPZnZXknLdzKwCwMwqJHVN2He0pDnAZuD7Zvb3hH2/l1QFTAbuMO+p4ZzLQO3yc2mXn0vfrrUnq62hZpWs+3rZhu28s3w9m3Z8ukGpc5u8T3VXr15fsGIT9/31PSpDL8LyjTu4Zcp8gAZNRg2WiCS9BHRPsutWM3s2vL8MeKxGPKcCnwG2AzNCz4sZKV62Aigys3WSRgDPSBpkZpuJmuXKJbUjSkRXAn+oJfargasBioqKUry0c841nratcujbtS19u7at9ZjtuxKaAWt0X1+xqZLZH29gw/ZkTz/227G7irunL2maicjMRtW1PzwPGgeMSNhcBrxqZmvDMdOA4UDNRLRKUo9QG+oBrA7X3AnsDO9nSXqfqJY108zKw/Ytkh4FTqSWRGRmk4BJEHXfTv2unXMuc7TOy+GYwrYcU1h7sqrcXbWv6e/Lv3sr6TErkoyAkU5xdsMYBSw2s7KEbdOBYkmtQ6I6A3g3SdnngKvC+6uAZwEkFUrKDu+PAY4DPpCUI6lL2J4LnA8saIB7cs65JiU/N5uju7ThlGO70KtjQdJjetayPV3iTESXcmCzHGa2AbgPeAeYC8w2s6kAkh6UVP3jqInAaElLgdFhHeB0YJ6kUuAp4BozWw+0AqZLmhfOWw78ruFuzTnnmp4bzu1PQY3ODNXDLjUkH1mhHj6ygnOuJUlXrzkfWcE559whGXtCr0afmLD5/lTXOedck+CJyDnnXKw8ETnnnIuVJyLnnHOx8kTknHMuVt59ux6S1gAfHWLxLsDaNIbjDp9/JpnJP5fMc7ifyVFmVpjKgZ6IGpCkman2o3eNwz+TzOSfS+ZpzM/Em+acc87FyhORc865WHkialiT4g7AfYp/JpnJP5fM02ifiT8jcs45FyuvETnnnIuVD3raQCR9CJxnZksStv1fYAXwNtEcSh8ArYFVwF1m9udw3HSgW8LpugKfmNlJjRN981bHZ1NBNGswQF+i6UJ2APOAT4AvJpymNXAs0MHMtjZC2M1WXd8VM7tL0qlE08O0BwT80sx+HY7z70oDqef/YU8Ai4Dqff80s2sknQQ8UONUA4DvmNlvar2YmfnSAAvwM+C/EtaziGagPQoYCfw5Yd8wYDlwdpLztAEWA6PjvqfmstT12SRsewUoqeMc/w+4I+57aQ5LPd+V7sDHwPCwrwswC7goyXn8u9J4n0sfYEEK5zgnfCZt6jrOm+YSSPqKpHmSSiX9UVI3SU+H9VJJpyQpc3s49m+Slkr6etj1GNHkf9VOB5ab2ad+HGtmc4EfAd9KEtYvgGlm9uJh32ATFtdnU0ssVxDVmG4/zNtqshrx8/gm8JCZzQYws7XAjcANScJq8d+VDPuedCGagPRyM9tW58FxZ91MWYBBRNXMLmG9M/A48J9hPZuoGaZmuduBUqCA6F9rnwA9w76FwNDw/rfAN8P7kSTUiMK2YcCiGtsuIppRtlXcf5+W8tkklH2FJDUion8JlgF94/67tITPA5gCjKlxng7AxhrbWvx3pZE/lz7ANmAO8CpwWpLzPgvclErsXiPa7yzgKYv+xYVFU4yfBfwmrFeZ2aZayj5rZjtC2ZeBE8P2x4BLJeUAY4An67i+DliRegG/BL5sZjsP8Z6ai7g/GwAkZQOPALeZ2bLDuaEmrjE/DwF1du3178o+jfm5VABFZnYCcB3wqKT21SeTdA3RM727UwncE9F+9f4HDyDpm5LmhqVn2FyzXPX6Y8AEYBQwz8xW13HqE4ge/iFJwMPARDN79yDuobmK+7Op9n2gwsx+n2LczVVjfh4LgZrDzIwAZoZr+Hdlv0b7XMxsp5mtC+9nAe8D/cL5BxB9V640s70pRR53dTJTFqJq7XvAEba/WvsnDqzWtq+lWjsXyAeOIHqw2jNh/9th/1cTto3kwM4KxcCHhM4KRO3ff07n/TXlpTE/m4R9r5DQNAecHD6jTnH/PeJeGvm70iMcNyysHwG8CZwR1v27Es/nUghkh/fHEPUw7QzkAbOBSw4q9rj/eJm0AFcBC4jaSx8i6hb6LDA/fBCfreVDnATMAJYCX6+x/ztAJQlts0SJaBNR++oS4O/ABQn7dxL9C2NuwvJy3H+flvDZJOx7hQMT0XSibvZzayzHxv23ae6fB9FD8rfDd2UncFnCPv+uxPC5ABcT1VZLiRLPBWH7ZcDuJN+T79QVt4+scJgk3Q5sNbN74o7FHcg/m8ySjs9D0jeBa4DTzWxDumJryTLhe+LPiJxzTYaZ/drMhngSal68RuSccy5WXiNyzjkXK09EzjnnYuWJyDnnXKw8ETnXwCT1kbQgDee5UNLN6Ygpybn/U1Lrhji3c/XxzgrONTBJfYh+dDk47lhqI2k50e+m1sYdi2t5vEbkXOPIkfRwGBn5qerah6QfSHpH0gJJk8KQNUj6D0nvhuP/FLZ9VdL94f0loUyppNdqXkxSD0mvhWFcFkg6LWw/R9KbkmZLelJSW0n/AfQEXpb0cmP9QZyr5jUi5xpYqBF9CJxqZq9L+l/gXTO7R1JniwanRNIfgSfM7HlJK4CjzWynpI5mtlHSV4lqLd+SNJ9o0rLy6v01rvldIN/MfhIGa20NtCIazfrzZrZN0k1Eo1X/yGtELk5eI3KucXxiZq+H948Ap4b3Z0p6KySWs4jGC4NoVtj/p2juoz1Jzvc68FCYOyY7yf53gH8Jv5ofYmZbiMbLOx54XdJcouFgjjrsO3PuMHkicq5xfGp0Y0n5wH8D481sCNEkYvlh/xeBXxONND0rDMO/v7DZNUQjHB8JzJV0RI39rxGN0VYO/FHSV4hGZ37RzIaF5Xgz+1pa79K5Q+CJyLnGUSTps+H9ZcA/2J901kpqC4wHkJQFHGlmLxPNRtoRaJt4MknHmtlbZvYDYC1RQkrcfxSw2sx+B/wPMBz4J/A5SX3DMa0l9QtFtgDt0ni/zqUsp/5DnHNpsAi4StIDRCMc/8bMtkv6HdHIyMuJmtMgamp7RFIHolrMz8MzosTz3S3puLB/BtEoyIlGAjdI2g1sBb5iZmvCc6bHJLUKx32faOqAScALkirM7Mw03rdz9fLOCs4552LlTXPOOedi5YnIOedcrDwROeeci5UnIuecc7HyROSccy5Wnoicc87FyhORc865WHkics45F6v/D/k2ASZXZW5BAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(my_basis_sets,energy)\n",
    "plt.plot(my_basis_sets,energy)\n",
    "plt.xlabel('basis set')\n",
    "plt.ylabel('total energy (AU)');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eed62018",
   "metadata": {},
   "source": [
    "# What is going on numerically?\n",
    "\n",
    "we can dig deeper into PySCF and extract the matrix representation of the Schrodinger equation and see numerically what some of the constituents are: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "3f8b4a88",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from pyscf import ao2mo\n",
    "\n",
    "mol_h2o = pyscf.gto.M(atom = 'O 0 0 0; H 0 1 0; H 0 0 1', basis = 'sto3G')\n",
    "\n",
    "# get the overlap matrix\n",
    "overlap = mol_h2o.intor('int1e_ovlp') \n",
    "\n",
    "# get the one-electron terms - kinetic energy and electron nuclear interaction\n",
    "hkinetic = mol_h2o.intor('int1e_kin') \n",
    "hnuc = mol_h2o.intor_symmetric('int1e_nuc')\n",
    "\n",
    "# get the two-electron terms - electron-electron repulsion elements\n",
    "eri_4fold_ao = mol_h2o.intor('int2e_sph')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "2659c828",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((7, 7), (7, 7, 7, 7))"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# check the shape of the matrices\n",
    "hkinetic.shape, eri_4fold_ao.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "0c7239c1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 2.90031999e+01 -1.68010939e-01  0.00000000e+00  0.00000000e+00\n",
      "   0.00000000e+00 -5.12015832e-03 -5.12015832e-03]\n",
      " [-1.68010939e-01  8.08127955e-01  0.00000000e+00  0.00000000e+00\n",
      "   0.00000000e+00  1.09023173e-01  1.09023173e-01]\n",
      " [ 0.00000000e+00  0.00000000e+00  2.52873120e+00  0.00000000e+00\n",
      "   0.00000000e+00  0.00000000e+00  0.00000000e+00]\n",
      " [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  2.52873120e+00\n",
      "   0.00000000e+00  2.52392544e-01  0.00000000e+00]\n",
      " [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00\n",
      "   2.52873120e+00  0.00000000e+00  2.52392544e-01]\n",
      " [-5.12015832e-03  1.09023173e-01  0.00000000e+00  2.52392544e-01\n",
      "   0.00000000e+00  7.60031884e-01  1.80967253e-02]\n",
      " [-5.12015832e-03  1.09023173e-01  0.00000000e+00  0.00000000e+00\n",
      "   2.52392544e-01  1.80967253e-02  7.60031884e-01]]\n"
     ]
    }
   ],
   "source": [
    "print(hkinetic)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "8a04322c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-6.16768093e+01 -7.43359610e+00  0.00000000e+00 -1.42231889e-02\n",
      "  -1.42231889e-02 -1.57330366e+00 -1.57330366e+00]\n",
      " [-7.43359610e+00 -1.00996038e+01  0.00000000e+00 -1.69721860e-01\n",
      "  -1.69721860e-01 -3.59672783e+00 -3.59672783e+00]\n",
      " [ 0.00000000e+00  0.00000000e+00 -9.95087124e+00  0.00000000e+00\n",
      "   0.00000000e+00  0.00000000e+00  0.00000000e+00]\n",
      " [-1.42231889e-02 -1.69721860e-01  0.00000000e+00 -1.00640216e+01\n",
      "   0.00000000e+00 -2.74859516e+00 -6.09772963e-02]\n",
      " [-1.42231889e-02 -1.69721860e-01  0.00000000e+00  0.00000000e+00\n",
      "  -1.00640216e+01 -6.09772963e-02 -2.74859516e+00]\n",
      " [-1.57330366e+00 -3.59672783e+00  0.00000000e+00 -2.74859516e+00\n",
      "  -6.09772963e-02 -5.70193488e+00 -1.77459817e+00]\n",
      " [-1.57330366e+00 -3.59672783e+00  0.00000000e+00 -6.09772963e-02\n",
      "  -2.74859516e+00 -1.77459817e+00 -5.70193488e+00]]\n"
     ]
    }
   ],
   "source": [
    "print(hnuc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "eb11386e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1.         0.23670394 0.         0.         0.         0.04875181\n",
      "  0.04875181]\n",
      " [0.23670394 1.         0.         0.         0.         0.44714688\n",
      "  0.44714688]\n",
      " [0.         0.         1.         0.         0.         0.\n",
      "  0.        ]\n",
      " [0.         0.         0.         1.         0.         0.37820346\n",
      "  0.        ]\n",
      " [0.         0.         0.         0.         1.         0.\n",
      "  0.37820346]\n",
      " [0.04875181 0.44714688 0.         0.37820346 0.         1.\n",
      "  0.28990152]\n",
      " [0.04875181 0.44714688 0.         0.         0.37820346 0.28990152\n",
      "  1.        ]]\n"
     ]
    }
   ],
   "source": [
    "print(overlap)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb6b0fc9",
   "metadata": {},
   "source": [
    "# Other systems"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "a12c3e91",
   "metadata": {},
   "outputs": [],
   "source": [
    "# here we create a H2O molecule and run a HF calculation on it using a minimal basis set (STO3G)\n",
    "atom_he = pyscf.gto.M(atom = 'He 0 0 0', basis = 'sto3g')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "8875441e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the one-electron terms - kinetic energy and electron nuclear interaction\n",
    "he_kin = atom_he.intor('int1e_kin') \n",
    "he_nuc = atom_he.intor_symmetric('int1e_nuc')\n",
    "\n",
    "# get the two-electron terms - electron-electron repulsion elements\n",
    "he_eri = atom_he.intor('int2e_sph')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "9cca4d73",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.]])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "he_overlap = atom_he.intor('int1e_ovlp') \n",
    "he_overlap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "cdc0a20d",
   "metadata": {},
   "outputs": [],
   "source": [
    "he_core = he_kin[0][0]+ he_nuc[0][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "61898885",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-52.56554983397122"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "he_core*27.2113845"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "aa26a1d3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "converged SCF energy = -2.80778395753997\n"
     ]
    }
   ],
   "source": [
    "rhf_he = pyscf.scf.RHF(atom_he)\n",
    "e_he = rhf_he.kernel()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "10329868",
   "metadata": {},
   "outputs": [],
   "source": [
    "atom_h = pyscf.gto.M(atom = 'H 0 0 0', basis = 'sto3g', spin = 1)\n",
    "h_kin = atom_h.intor('int1e_kin') \n",
    "h_nuc = atom_h.intor_symmetric('int1e_nuc')\n",
    "h_core = h_kin[0][0]+ h_nuc[0][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "ceaff8b2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-12.696338109024174"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "h_core*27.2113845"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "c814a04e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "converged SCF energy = -0.466581849557275  <S^2> = 0.75  2S+1 = 2\n"
     ]
    }
   ],
   "source": [
    "uhf_h = pyscf.scf.UHF(atom_h)\n",
    "e_h = uhf_h.kernel()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aebc8920",
   "metadata": {},
   "source": [
    "you can also set `charge` and other properties on your 'molecule'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8c12153",
   "metadata": {},
   "source": [
    "## Getting forces from PySCF\n",
    "\n",
    "as well as calculating the electronic structure we can calculate the forces on the nucelei from all the electrons and the other nuclei. These forces can then be used to move the atom to find favourable molecular geometries and / or to run molecular dynamics simulations. \n",
    "\n",
    "We can get the forces from a particular molecular configuration (like mol_h20 we created earlier) using the \n",
    "`mol.nuc_grad_method().grad_nuc()`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "31c43a5e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "converged SCF energy = -74.9611711378677\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[ 0.        ,  2.24022816,  2.24022816],\n",
       "       [ 0.        , -2.3392332 ,  0.09900503],\n",
       "       [ 0.        ,  0.09900503, -2.3392332 ]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grad = mol_h2o.nuc_grad_method().grad_nuc()\n",
    "grad"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e75bfad1",
   "metadata": {},
   "source": [
    "This is just an array of force vectors - just like we had in our MD code before. \n",
    "\n",
    "Once we have that we can update positions etc exactly as if we were using a classical forcefield (like LJ for instance)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}