16_02_VMC.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Variational Monte Carlo\n",
    "\n",
    "The main ingredient is a trial wave function $|\\psi _{T}(\\overline{\\eta }\n",
    ")\\rangle $, that depends on a set of parameters $\\overline{\\eta }$. This\n",
    "wave function is represented in terms of a basis of orthogonal states $\n",
    "|x\\rangle$ \n",
    "$$\n",
    "|\\psi _{T}(\\overline{\\eta })\\rangle =\\sum_{x}\\langle x|\\psi _{T}(\\overline{\n",
    "\\eta })\\rangle |x\\rangle =\\sum_{x}C_{x}(\\overline{\\eta })|x\\rangle ,\n",
    "$$\n",
    "where the coefficients of the parametrization are known functions of $%\n",
    "\\overline{\\eta }$. We would like the wave function to be a good\n",
    "representation of the actual ground state of a model. Finding the best wave\n",
    "function means finding the right set of parameters $\\overline{\\eta }$ that\n",
    "maximize the overlap with the actual ground state. In practice this is\n",
    "impossible since we do not know the groud state \\textit{a priori}, and\n",
    "some physical insight is needed to derive a good analytical\n",
    "approximation. Then we apply the variational principle, that states\n",
    "that the variational energy of the trial state is always greater or equal to\n",
    "the exact energy of the ground state:\n",
    "\\begin{equation}\n",
    "\\langle E\\rangle _{T}=\\frac{\\langle \\psi _{T}\\mid H\\mid \\psi _{T}\\rangle }{%\n",
    "\\langle \\psi _{T}\\mid \\psi _{T}\\rangle }\\geq E_{0},\n",
    "\\label{energia variacional}\n",
    "\\end{equation}\n",
    "and we use the criterium of minimizing the variational energy. In order to\n",
    "do that we require to calculate this quantity for different sets of\n",
    "parameters $\\overline{\\eta }$, and once we found a proper wave function, we\n",
    "can calculate the physical quantities of interest.\n",
    "\n",
    "The expectation value of an arbitrary operator $O$ is\n",
    "\\begin{eqnarray}\n",
    "\\langle O\\rangle _{Var} &=&\\frac{\\langle \\psi _{T}\\mid O\\mid \\psi\n",
    "_{T}\\rangle }{\\langle \\psi _{T}\\mid \\psi _{T}\\rangle }  \\nonumber \\\\\n",
    "&=&\\frac{\\sum_{x}\\langle \\psi _{T}\\mid x\\rangle \\langle x\\mid O\\mid \\psi\n",
    "_{T}\\rangle }{\\sum_{x}\\langle \\psi _{T}\\mid x\\rangle \\langle x\\mid \\psi\n",
    "_{T}\\rangle } \\\\\n",
    "&=&\\frac{\\sum_{x}\\left| \\langle \\psi _{T}\\mid x\\rangle \\right|\n",
    "^{2}\\langle x\\mid O\\mid \\psi _{T}\\rangle /\\langle x\\mid\n",
    "\\psi _{T}\\rangle }{\\sum_{x}\\left| \\langle \\psi _{T}\\mid x\\rangle \\right| ^{2}%\n",
    "}  \\nonumber \\\\\n",
    "&=&\\sum_{x}P_{x}O_{x},  \\label{valor medio}\n",
    "\\end{eqnarray}\n",
    "with\n",
    "\\begin{equation}\n",
    "P_{x}=\\frac{\\left| \\langle \\psi _{T}\\mid x\\rangle \\right| ^{2}}{%\n",
    "\\sum_{x}\\left| \\langle \\psi _{T}\\mid x\\rangle \\right| ^{2}},  \\label{P}\n",
    "\\end{equation}\n",
    "and\n",
    "\\begin{equation}\n",
    "O_{x}=\\frac{\\langle x\\mid O\\mid \\psi _{T}\\rangle }{\\langle x\\mid \\psi\n",
    "_{T}\\rangle }.  \\label{o}\n",
    "\\end{equation}\n",
    "The equation (\\ref{valor medio}) has precisely the form of a mean value in\n",
    "statistical mechanics, with $P_{x}$ as the Boltzmann factor:\n",
    "$$\n",
    "P_{x}\\geq 0;\\,\\sum_{x}P_{x}=1.\n",
    "$$\n",
    "The first step in order to calculate it consists of generating a collection\n",
    "of configurations distributed according to this probability. For that\n",
    "purpose we employ the Metropolis algorithm: starting from a configuration $\n",
    "|x\\rangle $, we accept a new configuration $|x^{\\prime }\\rangle $ with\n",
    "probability $R=\\left| \\langle \\psi _{T}\\mid x^{\\prime }\\rangle \\right|\n",
    "^{2}/\\left| \\langle \\psi _{T}\\mid x\\rangle \\right| ^{2}$, or $1/(1+R)$ if we\n",
    "use the heat bath approach.\n",
    "\n",
    "The variational simulations are simple to perform and very stable.\n",
    "Since the probabilities do not deppend on the statistics of the particles\n",
    "involved, they do not suffer from the sign problem. However,\n",
    "the results deppend decisively of the\n",
    "quality of the variational wave function, because they are completely\n",
    "pre-determined by it, and the physical arguments that define it. In the\n",
    "particular situation in which the trial function coincides with the exact\n",
    "ground state, the matrix elements (\\ref{o}) for $O=H$ are all equal to $%\n",
    "E_{0}:$%\n",
    "$$\n",
    "E_{x}=\\frac{\\langle x\\mid H\\mid \\psi _{T}\\rangle }{\\langle x\\mid \\psi\n",
    "_{T}\\rangle }=E_{0}.\n",
    "$$\n",
    "This is the property called \"zero variance\": the more the wave function\n",
    "resembles the actual ground state, the more rapidly the variational enery\n",
    "converges with the number of iterations.\n",
    "\n",
    "In general, the computation results more complicated, or numerically more\n",
    "expensive, when the number of \\ variational parameters $\\overline{\\eta }$\n",
    "that define the trial state is large. Thus, we always try to keep the form\n",
    "of the wave function simple enough and with only a few variational degrees\n",
    "of freedom."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quantum Spin Chain\n",
    "\n",
    "We seek a good representaiton of the ground state wave function for $L$ spins interacting via a Heisenberg Hamiltonian:\n",
    "$${H}=\\sum_{i=1}^{N-1}\n",
    "{S}^z_i {S}^z_{i+1} -\n",
    "\\frac{1}{2}\\left[\n",
    "{S}^+_i {S}^-_{i+1} +\n",
    "{S}^-_i {S}^+_{i+1}\n",
    "\\right]\n",
    "$$\n",
    "Notice the minus sign in the off diagonal term. This results from a unitary transformation called \"Marshall rotation\". As a consequence, all coefficients in the wave-function are positive and real, and we don't have to worry about complex phases nor signs.\n",
    "As we learned in previous lectures, we can only solve the problem extraclty for a small number of spins (in reality, there is an exact analytical solution called \"Bethe Ansatz\"). We will now learn how to solve the problem using a variational approximation to the ground state. Leaving aside the justification for our choice of wave function, we introduce the so-called Jastrow wave-function:\n",
    "$$\n",
    "|\\Psi\\rangle = \\sum_{\\bar{\\sigma}} \\exp\\left({\\alpha\\sum_{i>j}\\frac{S_i^zS_j^z}{|i-j|}}\\right)|\\bar{\\sigma}\\rangle,\n",
    "$$\n",
    "\n",
    "where $\\sigma=\\pm 1/2$ represents the orientation of the spins along the $z$ axis. This expression enjoys an elegant simplicity, and depends on only one free variational parameter $\\alpha$. We start our simulation by chosing a pair of spins; if anti-parallel, we atempt a spin flip. The acceptance probability is determined by the ratio $|\\Psi(\\bar{\\sigma}')|^2/|\\Psi(\\bar{\\sigma})|^2$, where $\\Psi(\\bar{\\sigma})=\\exp\\left({\\alpha\\sum_{i>j}\\frac{\\sigma_i \\sigma_j}{|i-j|}}\\right)$. Notice that this probability reminds us of a Boltzmann factor with inverse temperature $\\alpha$ and classical Ising-like Hamiltonian $H'=\\sum_{i>j}\\frac{S_i^zS_j^z}{|i-j|}$. Therefore, in order to implement the code, we only have to introduce minor modifications to our previous Monte Carlo code for the Ising model. There are some important changes: (i) Since we are only inetrested in the sector with zero total spin $S^z=0$, we pick spins in pairs and we attempt spin flips that preserve the total magnetization; (ii) The interactions in the effective classical Hamiltonian are long-range, (iii) more importantly, the expression for the energy involves the actual quantum Hamiltonian and we need to account for all possible spin flips to evaluate the expression for the \"local energy\":\n",
    "$$E_\\bar{\\sigma} = \\frac{\\langle \\bar{\\sigma}|H|\\Psi\\rangle}{\\langle \\bar{\\sigma}|\\Psi\\rangle}=\\sum_i\\frac{\\langle \\bar{\\sigma}|S^z_iS^z_{i+1}|\\Psi\\rangle}{\\langle \\bar{\\sigma}|\\Psi\\rangle}-\\frac{1}{2}\\sum_i\\frac{\\langle \\bar{\\sigma}|S^+_iS^-_{i+1}+S^-_iS^+_{i+1}|\\Psi\\rangle}{\\langle \\bar{\\sigma}|\\Psi\\rangle}$$\n",
    "The last term involves all the allowed spin flips between nearest neighbors, and has at most $L$ non-zero contributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from matplotlib import pyplot\n",
    "import matplotlib.animation as animation\n",
    "\n",
    "class Direction:\n",
    "    RIGHT, LEFT = range(2)\n",
    "\n",
    "def distance(r,L):\n",
    "    return abs(abs(r)-L*int(2*abs(r)/L))\n",
    "    \n",
    "class VMC(object):\n",
    "    \n",
    "    def __init__ (self, L, alpha):\n",
    "      \n",
    "        self.L = L\n",
    "        self.alpha = alpha\n",
    "        # Table for probabilities\n",
    "        self.psi2 = np.zeros(shape=(3,int(L/2+1)))\n",
    "        \n",
    "        # Initialize neighbors table for boundary conditions\n",
    "        self.nn = np.zeros(shape=(self.L,2), dtype=np.int16)\n",
    "\n",
    "        # Periodic boundary conditions\n",
    "        for ix in range(L):\n",
    "            self.nn[ix,Direction.LEFT] = ix-1\n",
    "            self.nn[ix,Direction.RIGHT] = ix+1\n",
    "\n",
    "        self.nn[0,Direction.LEFT] = L-1\n",
    "        self.nn[L-1,Direction.RIGHT] = 0\n",
    "\n",
    "        # Initialize spins\n",
    "        self.spin = np.ones(self.L, dtype=np.int16)\n",
    "        pos = np.random.choice(range(L), int(L/2), replace=False) # we pick L/2 unique random positions and flip\n",
    "        for ix in pos:\n",
    "            self.spin[pos] = -1\n",
    "        \n",
    "        self.set_alpha(alpha)\n",
    "                                                             \n",
    "    def set_alpha(self, alpha):\n",
    "        self.alpha = alpha\n",
    "        # Lookup tables for transition probabilities\n",
    "        # Pair-wise contribution to the argument in the exponential\n",
    "        for i in range(-1,2,2):\n",
    "            for l in range(1,int(self.L/2)+1):\n",
    "                self.psi2[i+1,l] = self.alpha*i*0.25/l\n",
    "        \n",
    "    def metropolis(self):\n",
    "        nchanges = 0\n",
    "\n",
    "        # trial spin change\n",
    "        # pick a random pair of spins\n",
    "        i = 0\n",
    "        j = 0\n",
    "        while(self.spin[i] == self.spin[j]):\n",
    "            i = int(np.random.random()*self.L)\n",
    "            j = int(np.random.random()*self.L)\n",
    "\n",
    "        si = self.spin[i]\n",
    "        sj = self.spin[j]\n",
    "\n",
    "        arg = 0.\n",
    "        for k in range(self.L):\n",
    "            sk = self.spin[k]\n",
    "            if(i != k and k != j):\n",
    "                arg += self.psi2[si*sk+1,distance(i-k,self.L)]\n",
    "                arg += self.psi2[sj*sk+1,distance(j-k,self.L)]\n",
    "        p = np.exp(-4*arg)\n",
    "        if(p >= 1 or np.random.random() < p):\n",
    "            self.spin[i] *= -1\n",
    "            self.spin[j] *= -1\n",
    "            nchanges += 1\n",
    "        return nchanges\n",
    "    \n",
    "    def Psi(self, spins):\n",
    "        arg = 0.\n",
    "        for i in range(self.L):\n",
    "            for j in range(i+1,self.L):\n",
    "                si = spins[i]\n",
    "                sj = spins[j]\n",
    "                arg += self.psi2[si*sj+1,distance(i-j,self.L)]\n",
    "        return np.exp(arg)\n",
    "                \n",
    "    def Eloc(self):\n",
    "        Ediag = 0.\n",
    "        Eoff = 0.\n",
    "        for i in range(self.L):\n",
    "            j = self.nn[i,Direction.RIGHT]\n",
    "            si = self.spin[i]\n",
    "            sj = self.spin[j]\n",
    "\n",
    "            # Diagonal contribution\n",
    "            Ediag += si*sj*0.25\n",
    "            # Off diagonal contribution\n",
    "            if(si != sj):\n",
    "                arg_flip = 0.\n",
    "                for k in range(self.L):\n",
    "                    sk = self.spin[k]\n",
    "                    if(i != k and k != j):\n",
    "                        arg_flip += self.psi2[si*sk+1,distance(i-k,self.L)]\n",
    "                        arg_flip += self.psi2[sj*sk+1,distance(k-j,self.L)]\n",
    "                Eoff -= 0.5*np.exp(-2*arg_flip)\n",
    "        return Ediag+Eoff\n",
    "            \n",
    "        \n",
    "    def Eloc2(self):\n",
    "        Ediag = 0.\n",
    "        Eoff = 0.\n",
    "        for i in range(self.L):\n",
    "            j = self.nn[i,Direction.RIGHT]\n",
    "            si = self.spin[i]\n",
    "            sj = self.spin[j]\n",
    "\n",
    "            # Diagonal contribution\n",
    "            Ediag += si*sj*0.25\n",
    "            # Off diagonal contribution\n",
    "            if(si != sj):\n",
    "                aux_spin = np.array(self.spin)\n",
    "                aux_spin[i] *= -1\n",
    "                aux_spin[j] *= -1\n",
    "                Eoff -= 0.5*self.Psi(aux_spin)/self.Psi(self.spin)\n",
    "        return Ediag+Eoff\n",
    "            \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-1 -1  1  1 -1  1]\n"
     ]
    }
   ],
   "source": [
    "system = VMC(6,-2)\n",
    "for i in range(10):\n",
    "    system.metropolis()\n",
    "print(system.spin)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-2.6056655895852012 -2.6056655895852012\n"
     ]
    }
   ],
   "source": [
    "print(system.Eloc(),system.Eloc2())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-6.0 -2.2476454996991335\n",
      "-5.368421052631579 -2.3554867747015096\n",
      "-4.7368421052631575 -2.519872247494779\n",
      "-4.105263157894737 -2.6241387216420153\n",
      "-3.473684210526316 -2.7209015824737466\n",
      "-2.8421052631578947 -2.77675008036505\n",
      "-2.210526315789474 -2.7887071207778544\n",
      "-1.578947368421053 -2.7517243550783483\n",
      "-0.9473684210526319 -2.625733308156553\n",
      "-0.3157894736842106 -2.2854880822726833\n",
      "0.3157894736842106 -1.7900110468615151\n",
      "0.947368421052631 -1.1871430661525704\n",
      "1.5789473684210522 -0.586573018479168\n",
      "2.2105263157894726 -0.17318405340705434\n",
      "2.842105263157894 0.11417812064050045\n",
      "3.473684210526315 0.2702146123785564\n",
      "4.105263157894736 0.3600973675739084\n",
      "4.7368421052631575 0.4243142224346528\n",
      "5.368421052631579 0.44413201752312537\n",
      "6.0 0.4640753048997006\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "from matplotlib import pyplot\n",
    "\n",
    "alphas=np.linspace(-6,6,20)\n",
    "e0 = []\n",
    "for a in alphas:\n",
    "    Etot = 0.\n",
    "    Nmeas = 0\n",
    "    system.set_alpha(a)\n",
    "    for i in range(50000):\n",
    "        system.metropolis()\n",
    "        if(i%10 == 0):\n",
    "            Etot += system.Eloc()\n",
    "            Nmeas += 1\n",
    "    print(a,Etot/Nmeas)\n",
    "    e0.append(Etot/Nmeas)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEKCAYAAAAFJbKyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmUVOW57/HvwyDYigZFARm6HeMYh9shJs5xwqORE42KoAkhSGIc0JibqNzoiYboiYlDlKOCQ9S00aghMQmR4WgkOU40QiSIA0cZHQBBjczIc/94q1ZVN1XdVd1d9dbw+6zVq/fetbvqqVbq13u/k7k7IiIiueoUuwARESkvCg4REcmLgkNERPKi4BARkbwoOEREJC8KDhERyYuCQ0RE8qLgEBGRvCg4REQkL11iF1AIvXr18rq6uthliIiUjVmzZq10911yObcig6Ouro7GxsbYZYiIlA0zW5TrubpVJSIieVFwiIhIXqIGh5kNNrPXzWyBmV2Z4fERZrbCzOYkvkbFqFNERFKitXGYWWdgPHAisBSYaWZPuvurzU591N0vLnqBIiKSUcwrjkHAAnd/y903Ao8AQyLWIyIiOYgZHP2AJWn7SxPHmjvTzF4xs8fNbEBxShMRKR8NDVBXB506he8NDYV9vVJvHP8jUOfunwOmAQ9kO9HMRptZo5k1rlixomgFioi0V1s/+N1h4kQYNQoWLQr7ixbB6NGFDQ+LtXSsmX0R+A93PzmxfxWAu9+Q5fzOwCp337G1566vr3eN4xCRctDQED7o165NHevWDS66CA44AD74AFauDN/Tt1euhFWr4NNPMz9vbS0sXJh7HWY2y93rczk35gDAmcDeZrY7sAwYCgxLP8HM+rr7u4nd04H5xS1RRKR1DQ0wdiwsXgwDB8KPfwwnnrj1h37zD/8PPoDGRti8uenzbdgAN9/cvpoWL27fz7ckWnC4+2YzuxiYAnQG7nP3eWZ2HdDo7k8Cl5rZ6cBmYBUwIla9IiLNrVsH48bBz34GmzaFY4sWwYgRxavBLNyiam7gwMK9ZtQpR9x9MjC52bFr0ravAq4qdl0iIpl88gk89xw8+yzMmAEvvQQbNxbmtbp3h7PPhp13hl69wvf07V69YKed4PHHt77VVVMTAq1QKnKuKhGRjrB6Nfz97yEkZsyAWbOytylksv/+qQ/85h/8ya/nn4drrw1XL0k1NTBhAgwf3vprJM9Jv1U2blxuP9tWCg4RqXrpbRQ77wyHHQbvvw+vvJL5NlC6Ll22bqOA0Dg9b17rr33kkbDbbu374B8+vLBB0ZyCQ0Sq2i23wA9/mGqjWLkSpk7NfK4ZfO5zcPTRcMwxcNRRMG1a+28VFfuDv70UHCJSdRYtgt/+Fh59NNx+yqZz53D1ccwxISyOPBJ69mx6ToxbRbEpOESkKrzzDjz2WAiL559v/Xyz0MbRo0fr55bbFUN7KThEpGItXw5PPBHCYsaM1tsr0g0cmFtoVCMFh4iUvfTG7X79YPDgMGr66adhy5atz+/SBU44Ac45J7RtXHZZcbuzljsFh4iUteZTdixdCvfcs/V5nTrBscfC0KFwxhmh91RSTU11tVG0V7S5qgpJc1WJVI++feG997I/fuSRISzOPBP69CleXeWmXOaqEhFpsxUr4OqrWw6NJUugf//i1VQtFBwiUlY2b4a774b/9//gww+zn1dbq9AoFAWHiJSNv/0NLr44jOhO17lz06lA1LhdWKW+kJOICO+8Exqrjz66aWjsuSf88Y/wwAPhCsMsfM91nidpG11xiEjJ2rgRbr0Vrr8+zEyblOwF9b3vhVlkQUFRTAoOESlJU6bApZfCG280PX722fDzn8OAAXHqEgWHiJSYt98OVxK//33T4wccALffDscdF6cuSVEbh4hE19AQBt6ZwR57NA2NHXYIM9jOnq3QKBW64hCRqBoaYNQoWL9+68e++U244Qbo3bv4dUl2Cg4RieqHP8wcGn36wH33Fb8eaZ1uVYlINKtXw7JlmR97//3i1iK5U3CISBRr18JXvpL98YEDi1eL5EfBISJFt2kTnHUW/M//ZH5cI79LW9TgMLPBZva6mS0wsyszPN7NzB5NPP6imdUVv0oR6UhbtsDIkTB5curYeedp5Hc5idY4bmadgfHAicBSYKaZPenur6ad9i1gtbvvZWZDgf8Ezil+tSLSEdzhiivg179OHRs7Fn7yk3g1Sf5iXnEMAha4+1vuvhF4BBjS7JwhwAOJ7ceB483MilijiHSgG28MU4gkffvbYToRKS8xg6MfsCRtf2niWMZz3H0z8BGwMxmY2WgzazSzxhUrVhSgXBFpj4kTw/oZSV/7GowfH25PSXmpmMZxd5/g7vXuXr/LLrvELkdE0vzud/Cd76T2jz8+3K7q3DleTdJ2MYNjGZA+TVn/xLGM55hZF2BH4IOiVCciHeLpp+Hcc0OjOEB9PUyaBN26xa1L2i5mcMwE9jaz3c1sG2Ao8GSzc54EvpHY/hrwtFfiIukiFWrWLBgyJEyPDrDPPqE3VY8eceuS9onWq8rdN5vZxcAUoDNwn7vPM7PrgEZ3fxK4F3jIzBYAqwjhIiJl4I034JRTUuto9OsHU6eC7iSXv6hzVbn7ZGBys2PXpG2vB84qdl0i0j7LlsFJJ0Gyn0rPnmF9jdrauHVJx6iYxnERKQ2rVsHJJ8OiRWG/pibcnjrggLh1ScdRcIhIh1mzBk47DebNC/tdusATT8Dhh8etSzqWgkNE2q2hIdyG2n57eP75cMwMHnwQBg+OW5t0PK3HISLt0tAAo0eH2W7TnXde6IYrlUdXHCLSLmPHbh0aADNmFL8WKQ4Fh4i0y+LF+R2X8qfgEJE2c88+AlwLMVUuBYeItNlvfpN5vXAtxFTZFBwi0iYrV8KYMan9Hj20EFO1UK8qEWmTK64I4QEwYEAYu6E5qKqDrjhEJG/TpoUxGkl33qnQqCYKDhHJy9q1YeW+pHPOgVNPjVePFJ+CQ0Tycu218PbbYbtnT7jttrj1SPEpOEQkZy+/DDffnNr/xS+gd+949UgcCg4RycnmzTBqVGolvy9/GUaMiFqSRKLgEJGc3HorzJ4dtrt3h7vvDt1vpfooOESkVW+9Bddck9q/9lrYa6949UhcCg4RaZF76EW1bl3YP/jgMIZDqpeCQ0Ra9NBDMH162O7UCSZOhK5d49YkcSk4RCSr5cvh8stT+2PGwOc/H68eKQ0KDhHJ6vLLwxriEOaguu66uPVIaYgSHGa2k5lNM7M3E997ZjnvUzObk/h6sth1ilSzv/wFHn44tX/XXWFpWJFYVxxXAv/t7nsD/53Yz2Sdux+S+Dq9eOWJVLdPPoELL0ztDxumtcMlJVZwDAEeSGw/APx7pDpEJINrroFFi8L2zjuHMRwiSbGCo7e7v5vYfg/INmlBdzNrNLMXzEzhIlIEM2c2nX/q5pthl13i1SOlp2DrcZjZdKBPhofGpu+4u5uZZ3maWndfZmZ7AE+b2Vx3/98srzcaGA0wUGtWirTJpk1NpxU58UQ4//y4NUnpKVhwuPsJ2R4zs/fNrK+7v2tmfYHlWZ5jWeL7W2b2V+BQIGNwuPsEYAJAfX19tiASkRb84hfwyithe9ttQ4O4phWR5mLdqnoS+EZi+xvAH5qfYGY9zaxbYrsXcATwatEqFKkiDQ3Qrx9cdVXq2I9/DHvsEa8mKV2xguNG4EQzexM4IbGPmdWb2T2Jc/YDGs3sH8AzwI3uruAQ6WANDTB6NLzzTuqYGfTJdKNZBDD3yrurU19f742NjbHLECkLdXWpHlTpamth4cJiVyOxmNksd6/P5VyNHBepcosX53dcRMEhUuV23TXzcXVOlGwUHCJVbsCArY/V1MC4ccWvRcqDgkOkir39dlhHPMkstG1MmADDh8erS0pbwcZxiEjpu/32poP9pk6NW4+UB11xiFSpjz+Ge+5J7V92WbxapLwoOESq1K9+Bf/6V9j+7Gc1+63kTsEhUoU+/bTpRIZjxoRlYUVyof9VRKrQn/4Eb70Vtnv2hK9/PW49Ul4UHCJV6JZbUtujR8N228WrRcqPgkOkysyeDc8+G7Y7d4aLLopbj5QfBYdIlUlv2zjrrMwDAEVaouAQqSLvvQe/+U1qX11wpS0UHCJV5M47YePGsH344fCFL8StR8qTgkOkSqxfH4Ij6fLL49Ui5U3BIVIlHn4YVqwI2wMGwBlnxK1HypeCQ6QKuMOtt6b2L74YumimOmkjBYdIFXjmGZg7N2zX1MAFF8StR8qbgkOkCqQP+BsxIowWF2krBYdIhXvzzTDFSNKll8arRSpDTsFhZr8zs1PNTEEjUmZ++cvU9qmnhplwRdoj1yD4L2AY8KaZ3Whm+l9PpAx8+CHcf39qXwP+pCPkFBzuPt3dhwOHAQuB6Wb2nJl908y65vuiZnaWmc0zsy1mVt/CeYPN7HUzW2BmV+b7OiLV7p57YM2asH3ggXD88XHrkcqQ860nM9sZGAGMAmYDtxGCZFobXvefwBnAjBZerzMwHjgF2B8418z2b8NriVSlzZvD0rBJl10W1hQXaa+cenKb2STgs8BDwFfc/d3EQ4+aWWO+L+ru8xPP29Jpg4AF7v5W4txHgCHAq/m+nkg1mjQJFi8O2716wbBhceuRypHrEKBfuvszmR5w96y3mtqpH7AkbX8pkHVmHTMbDYwGGDhwYIFKEikf6QP+LrwQtt02Xi1SWXINjp5m1nyCgo+Aue6+PNMPmNl0oE+Gh8a6+x/yqDEn7j4BmABQX1/vHf38IuXkpZfguefCdteuIThEOkquwfEt4ItA8qrjWGAWsLuZXefuDzX/AXc/oZ21LQPSVwronzgmIq1Iv9oYOhT69o1Xi1SeXIOjK7Cfu78PYGa9gQcJt45mENo+OtpMYG8z250QGEMJXYJFpAVLl8Jjj6X21QVXOlquvar6J0MjYTkwwN1XAZvyfVEz+6qZLSVcxfzZzKYkju9mZpMB3H0zcDEwBZgP/Nbd5+X7WiLVZvz40KMK4Oij4bDD4tYjlSfXK46/mtmfgOTfMWcmjm0HfJjvi7r7JGBShuPvAP+Wtj8ZmJzv84tUq7Vr4e67U/u62pBCyDU4LiKMuzgysf8g8IS7O3BcIQoTkfw99BCsXh22d98dTj89bj1SmVoNjsRAvOnufhzwROFLEpG22LKlaaP4pZdC587x6pHK1Wobh7t/Cmwxsx2LUI+ItNHUqfDaa2G7Rw8YOTJuPVK5cr1V9Qkw18ymAWuSB91dEzSLlICGBhg1KrV/xBGwww7x6pHKlmtw/C7xJSIlJhka69enjj37bDg+fHi8uqRyWWjfzuFEs22Bge7+emFLar/6+npvbMx7Ci2RslRXB4sWbX28thYWLix2NVKuzGxWrlNI5bqQ01eAOcBTif1DzOzJtpcoIh0lOZFhrsdF2ivXAYD/QZit9kMAd58D7FGgmkQkDztm6baiuT6lUHINjk3u/lGzY1s6uhgRyc+mLPM21NTAuHHFrUWqR67BMc/MhgGdzWxvM7sdeK6AdYlIDh5/PCwPC9CpU1ioqbYWJkxQw7gUTq7BcQlwALAB+A3wMaDJDEQiSx/wd+21YRDgwoUKDSmsnLrjuvtaYGziS0RKwAsvhHU3ALbZBr7znbj1SPXIdenYfYDvA3XpP+PuXy5MWSLSmvSrjeHDYddd49Ui1SXXAYCPAXcB9wCfFq4cEcnFkiWhfSNpzJh4tUj1yTU4Nrv7nQWtRERyNn48fJr4E+644+Dgg+PWI9Ul18bxP5rZd82sr5ntlPwqaGUiktGaNaHXVJLW3JBiy/WK4xuJ7/837ZijQYAiRZe+5sYee8Cpp8atR6pPrr2qdi90ISLSui1b4LbbUvtac0NiaPFWlZn9IG37rGaP/bRQRYlIZs3X3PjmN+PWI9WptTaOoWnbVzV7bHAH1yIirUjvgvutb2nNDYmjteCwLNuZ9kWkgF59FaZMCdtmcMklceuR6tVacHiW7Uz7OTOzs8xsnpltMbOs87+b2UIzm2tmc8xMC2xIVfvlL1PbQ4aEhnGRGFprHD/YzD4mXF1sm9gmsd+9Ha/7T+AM4O4czj3O3Ve247VEyt4HH8CDD6b21QVXYmoxONy9IP013H0+gJnudonkYuJEWLcubB9yCBx9dNx6pLrlOgAwFgemmtksMxsduxiRGDZtgjvuSO1fdllo4xCJJdcBgHkzs+lAnwwPjXX3P+T4NEe6+zIz2xWYZmavufuMLK83GhgNMFBLn0kFeeIJWLYsbO+6Kwwd2vL5IoVWsOBw9xM64DmWJb4vN7NJhOVrMwaHu08AJgDU19e3ueFepNSkd8H97nehW7d4tYhACd+qMrPtzKxHchs4idCoLlI1XngBXnwxbGvNDSkVUYLDzL5qZkuBLwJ/NrMpieO7mdnkxGm9gb+b2T+Al4A/u/tTMeoViSV9epFhw6B373i1iCSZe+Xd1amvr/fGRg37kPK2dCnU1aWmT589O/SoEikEM5vl7lnH1aUr2VtVItUufc2NY49VaEjpUHCIlKC1a+HutOGxGvAnpUTBIVKCmq+5cdppcesRSafgECkxW7Y07YKrNTek1Cg4RErMtGlac0NKm4JDpMSkX22MHKk1N6T0KDhESsj8+fBUYrSS1tyQUqXgECkh6WtunH467LlnvFpEslFwiJSIVavggQdS++qCK6VKwSFSAhoawtVFcs2NAQPgmGPi1iSSjYJDJLKGBhg9Gj78MHXs/ffh4Yfj1STSEgWHSGRjx4aR4uk2bgzHRUqRgkMkssWL8zsuEpuCQySynj0zH9dCllKqFBwiEa1bF6YYaa6mBsaNK349IrlQcIhEdPvtqUbxTp3CoL/aWpgwAYYPj1ubSDYFW3NcRFq2ejXccENqf/x4LQ0r5UFXHCKR3HRT6mpjzz3hW9+KW49IrhQcIhG8+27TyQx/8hPo2jVePSL5UHCIRHD99alR4occAmefHbcekXwoOESKbMECmDgxtX/DDaFhXKRcRPnf1cxuMrPXzOwVM5tkZp/Jct5gM3vdzBaY2ZXFrlOkEH70I9i8OWwfeyycfHLUckTyFuvvnGnAge7+OeAN4KrmJ5hZZ2A8cAqwP3Cume1f1CpFOtjs2fDII6n9G24IXXBFykmU4HD3qe6e+JuLF4D+GU4bBCxw97fcfSPwCDCkWDWKFMLVV6e2//3f4fDD49Ui0lalcGd1JPCXDMf7AUvS9pcmjomUpb/+NbW6X6dOoSeVSDkq2ABAM5sO9Mnw0Fh3/0PinLHAZqChA15vNDAaYKAm+ZES4w5Xpd2Q/frX4YAD4tUj0h4FCw53P6Glx81sBHAacLy7e4ZTlgED0vb7J45le70JwASA+vr6TM8nEs2TT8ILL4TtbbaBH/84bj0i7RGrV9Vg4AfA6e6+NstpM4G9zWx3M9sGGAo8WawaRTrKp582bdu46CLNfCvlLVYbxx1AD2Camc0xs7sAzGw3M5sMkGg8vxiYAswHfuvu8yLVK9JmDz0Er74atnv0aHrLSqQcRZnk0N33ynL8HeDf0vYnA5OLVZdIR1u/Hq69NrX//e/DLrvEq0ekI5RCryqRinXXXamV/HbZBS6/PG49Ih1BwSFSIB9/3HQxph/9KNyqEil3Cg6RArn5Zli5MmzX1cHo0VHLEekwCg6RAli+HH7xi9T+dddBt27x6hHpSAoOkQIYNw4++SRsH3ggDBsWtx6RjqTgEOlgCxfCnXem9m+4ATp3jlaOSIdTcIh0sGuvhU2bwvYRR8Cpp8atR6SjKTgSGhpCA2anTuF7Q7tnz5Jq9M9/hgF/STfeqGnTpfJEGQBYahoaYNSoMFgLYNGiVA+Y4cPj1SXlZ+zYMKEhhCuNI4+MW49IIeiKg/CPPRkaSWvXwpVac1By1NAAffqEyQyTfvrTePWIFJKCg9TI3uaWLoWzz4Znnkn9FSnSXENDuEJ9//3Usc6dYe7ceDWJFJKCg5ZnKn3sMfjyl2G//eDWW2H16uLVJeVh7NhwhZru00/DcZFKpOAg9LmvqWl6rFOz38zrr4d5hvr1g5EjYebM4tUnpS3bFWu24yLlTsFBaACfMAFqa0MPmNpaePBBeOUV+O53m84vtG4d3H8/DBoE9fVwzz2wZk282iWudeugS5YuJlpzQyqVgiNh+PAwcGvLlvB9+HA46CAYPx6WLYO774ZDDmn6M7NmwQUXhKuQk0+G3XZTd95qc8klqTEb6Wpqmk5wKFJJFBw56NEjNH6+/DI8/3xYLzp93qGPPoKpU+Hdd0MjerI7r8Kjst17b/hK2mmn1BXrhAnqyi2VyzIv913e6uvrvbGxsaCv8cEH8MADYb2FN9/MfM6AAbrPXalefhm+9CXYsCHsn3deuL2pwX5SrsxslrvX53KurjjaaOed4Xvfg9dey/5hsWQJPP10ceuSwlu1Cs48MxUaBx0UbmUqNKRaKDjaqVOnlhtBjz8+3NpasaJ4NUnhbNkC558f2sEAdtgBnnhi6155IpVMwdEBMnXnTffQQ7DvvuF++JYtxatLOt5PfwqTJ6f2f/Ur2HvvaOWIRKHg6ACZuvPecQecc07qnFWrwnxYxx4Lr74arVRph2nT4JprUvs/+AF89avx6hGJRY3jBfbUU2EsyNtvp4517Ro+dMaOhW23jVeb5G7xYjjssNApAuCYY2D69OxjOETKTck3jpvZTWb2mpm9YmaTzOwzWc5baGZzzWyOmZVGEuRp8OAw1faVV6Y+ZDZtCre3Djoo/BUrpW3DBjjrrFRo9O0Ljzyi0JDqFetW1TTgQHf/HPAGcFUL5x7n7ofkmoSlqKYmrAKX7MKZ9L//CyedFG51pU+QJ6Xl8svhpZfCdpcuYf6yPn3i1iQSU5TgcPep7r45sfsC0D9GHcV20EHwt7+FrpufSbvGevhh2H330MVXI89Ly0MPNV0G9qabwqp+ItWsFBrHRwJ/yfKYA1PNbJaZjW7pScxstJk1mlnjihLu+9qpUxhV/tprMGxY6vi6daEBPTny/IILFB6xzZ0L3/52av/ss2HMmHj1iJSKgjWOm9l0INMF/Vh3/0PinLFAPXCGZyjEzPq5+zIz25Vwe+sSd5/R2muXUuN4a6ZODSvFbd689WPbbx+uUJrPkSWF99FHYRLLBQvC/r77httV6RNeilSSfBrHC9a85+4ntPS4mY0ATgOOzxQaiedYlvi+3MwmAYOAVoOjnJx0Uli7IZNPPoFDDw0fYBdcAOeeqw+uYnCHESNSobHddmGQn373IkGsXlWDgR8Ap7v72iznbGdmPZLbwEnAP4tXZfG0Nv12Y2O4ZdK3bwiQl17SioSFdNNN8Pvfp/bvvRf23z9ePSKlJlYbxx1AD2BaoqvtXQBmtpuZJcfl9gb+bmb/AF4C/uzuT8Upt7AyjTzv1i30wEqfhXfNmrD+xxe+EG5f3XEHfPhhcWutdH/9K1yV1sdvzJimAzlFRAMAS0ZDQxgQuHhxuAIZNy500121KvTsmTgR5s3b+ue23TaMMaitDbP1LlnS9OclNw0N8MMfhrVXkr70pbDe/DbbxKtLpFjyaeNQcJQJ97AWyMSJ8OijoRdWS2pqtCZErhoawi3A5r/T22+Hiy+OU5NIsSk4KjA40n30URj7MWECzJmT/bxevcIVjKY1yW716jCG5qOPtn6stjY1C65IpSv5KUekfXbcES68MIxEnzkz+3krV4YRzqNGwYwZmpk33fvvh2lgamszhwZoES6RbBQcZcwsdNWtrc1+zscfh15BxxwDe+4ZZnfNtmJhNVi8OKwTXlcH//mf8K9/ZT+3td5uItVKwVEBMvXK6toVevduemzhQrj+ethnn9Dwe9ddofG9Grz+OowcGcLzjjtg/frUY337bt0AXlMTfq8isjUFRwXItB7I/ffDu++GBvULL4SePZv+TPJ4377wta/BFVeEn6u0ubLmzAndaffbL/xO0kfoH3ZYGNi3dCncd1/T3586Fohkp8bxKrFhA/z5z/Dgg+F7pilO0nXvHj48zz+/OPV1tOeeC1cM6av1JR11VOj6fNJJWidcJEm9qhQcLVq5Mqwn8eCDLTeum8GRR8KgQamv5F/lpaahAa6+Ooxj2WabEJTNnXJKGNx31FHFr0+k1Ck4FBw5mz8/v+k0dt21aZB8/vOw007ZBzAWgjssXx7WM1mwACZNgj/+MfOcX2Zw5pkhVA49tDD1iFQCBYeCIy91dWEq97bq3TtcxaR/cHfvHhrizz03TBK4/fYtr5jXPHh+8hM47rgQDJm+Pvmk9bq22y7M87Xvvm1/byLVQsGh4MhLQ0NYI2Rt2nSTNTXws59B//5hUsXk18cft/11unVLhUjy+/bbh0F4c+dmnyW4rcw0dkUkVyUxrbqUj+QtpWy3moYMCd+3bIE33mgaJHPmhDXUc7FhQ/jqiC7APXrAXnuFrylTMgeaxmGIFIauOKRd1q8PU3a8997Wj3XpArvsEm4rrVnTtr/+Bw1KBcRee4VxGHvtFZ432Uif7YpJXWpFcqcrDima7t3h5z9v/YPbPYRMMkTSv597LmRa7be2Fl58sfUaWrtiEpGOpeCQdsvlg9ssTLa47bbhaiHdLbdkDp58Rm4PH66gECkWBYd0iPZ8cOuKQaS8KDikJOiKQaR8aK4qERHJi4JDRETyouAQEZG8KDhERCQvCg4REclLRY4cN7MVQFun7esFrOzAcmKqlPdSKe8D9F5KUaW8D2jfe6l1911aP61Cg6M9zKwx12H3pa5S3kulvA/QeylFlfI+oHjvRbeqREQkLwoOERHJi4JjaxNiF9CBKuW9VMr7AL2XUlQp7wOK9F7UxiEiInnRFYeIiORFwZGFmV1iZq+Z2Twz+1nsetrDzK4wMzezXrFraSszuynx3+MVM5tkZp+JXVM+zGywmb1uZgvM7MrY9bSVmQ0ws2fM7NXEv40xsWtqLzPrbGazzexPsWtpDzP7jJk9nvh3Mt/Mvlio11JwZGBmxwFDgIPd/QDg55FLajMzGwCcBCyOXUs7TQMOdPfPAW8AV0WuJ2dm1hkYD5wC7A+ca2b7x62qzTYDV7j7/sDhwEVl/F6SxgCx7FFfAAAD/ElEQVTzYxfRAW4DnnL3fYGDKeB7UnBkdiFwo7tvAHD35ZHraY9bgB8AZd2Y5e5T3X1zYvcFoH/MevI0CFjg7m+5+0bgEcIfJmXH3d9195cT2/8ifDj1i1tV25lZf+BU4J7YtbSHme0IHA3cC+DuG939w0K9noIjs32Ao8zsRTN71sw+H7ugtjCzIcAyd/9H7Fo62EjgL7GLyEM/YEna/lLK+MM2yczqgEOBHBb4LVm3Ev6w2hK7kHbaHVgB3J+47XaPmW1XqBer2oWczGw60CfDQ2MJv5edCJfinwd+a2Z7eAl2QWvlfVxNuE1VFlp6L+7+h8Q5Ywm3SxqKWZs0ZWbbA08Al7n7x7HraQszOw1Y7u6zzOzY2PW0UxfgMOASd3/RzG4DrgR+VKgXq0rufkK2x8zsQuB3iaB4ycy2EOaAWVGs+nKV7X2Y2UGEv0L+YWYQbu28bGaD3P29IpaYs5b+mwCY2QjgNOD4UgzxFiwDBqTt908cK0tm1pUQGg3u/rvY9bTDEcDpZvZvQHdgBzP7tbufF7mutlgKLHX35NXf44TgKAjdqsrs98BxAGa2D7ANZTYJmrvPdfdd3b3O3esI/2MdVqqh0RozG0y4pXC6u6+NXU+eZgJ7m9nuZrYNMBR4MnJNbWLhr5B7gfnufnPsetrD3a9y9/6Jfx9DgafLNDRI/LteYmafTRw6Hni1UK9XtVccrbgPuM/M/glsBL5RZn/hVqI7gG7AtMQV1Avu/p24JeXG3Teb2cXAFKAzcJ+7z4tcVlsdAZwPzDWzOYljV7v75Ig1SXAJ0JD44+Qt4JuFeiGNHBcRkbzoVpWIiORFwSEiInlRcIiISF4UHCIikhcFh4iI5EXBIdLBzGxhazMR53KOSKlScIiISF4UHCLtYGa/N7NZibUpRjd7rC6xNkJDYn2Ex82sJu2US8zsZTOba2b7Jn5mkJk9n5io7rm0kcAiJUPBIdI+I939/wD1wKVmtnOzxz8L/Je77wd8DHw37bGV7n4YcCfw/cSx14Cj3P1Q4BrgpwWtXqQNFBwi7XOpmf2DsEbIAGDvZo8vcff/SWz/Gjgy7bHkBIGzgLrE9o7AY4npbm4BDihE0SLtoeAQaaPEVNwnAF9094OB2YRZVtM1n9MnfX9D4vunpOaNux54xt0PBL6S4flEolNwiLTdjsBqd1+baKM4PMM5A9PWfh4G/D2H50xOuT6iQ6oU6WAKDpG2ewroYmbzgRsJt6uae52wLvd8oCehPaMlPwNuMLPZaPZqKVGaHVekQBJLq/4pcdtJpGLoikNERPKiKw4REcmLrjhERCQvCg4REcmLgkNERPKi4BARkbwoOEREJC8KDhERycv/B6AaB48JhHVgAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pyplot.xlabel(\"alpha\")\n",
    "pyplot.ylabel(\"Energy\")\n",
    "pyplot.plot(alphas,e0,'bo-',lw=3);\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}