Variational/Variational_Hydrogen.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sympy import *\n",
    "from sympy.abc import r,x,y,z\n",
    "from scipy.integrate import quad, nquad\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "init_printing()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Energy of the Hydrogen Atom\n",
    "The variational principle states a trial wavefunction will have an energy greater than or equal to the ground state energy.\n",
    "$$\\frac{\\int \\psi H \\psi}{ \\int \\psi^2}  \\ge E_0$$\n",
    "\n",
    "First consider the hydogen atom. Let us use a trial wavefunction that is not the exact ground state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAEMAAAAbCAYAAAAnFzLpAAAABHNCSVQICAgIfAhkiAAAAfBJREFU\nWIXt10+ITWEYx/GPa65JkjJI/oVbFpKIpoRZTWHhbygmZaUYykpKiaUp9krYWYiFsrBQIoWS2Cg1\nGwurqamx8idZvO81Z7j3Duc67xm53zqd55zzvuf3O2/v+7zPocMPppVtoAkVDGJmvB4q0Uvp7MTS\nGN/BhhSilRQiOajhUIyHjQ/Mf0k3Zsf4ARalEO1KIRLZjvn4jI04F+NGfIrHVjzChwT+krHDxKl+\nD+sn6TMH5wtzVBIrTUyAq/EU1Un6DcY2VfQXYy09B+L5LK5gBMuxBKfxBLdxEaPCdnoYY7HtKNb8\nJS+tNJNQH4wjQp54icvYJuxmrzAQ4xUFe2mpmafoOiVsfc14jlsxXohleJF5fhJ9OIi5eIcF+Fag\njyztaLbFfr/WMtdxLPP8bkpDrTSLLrq2CDOjzlph3V6L1/14WLCHn2mq2ajOWIcTQmKpYp6Q2T/m\nEH6GfbFvBT3YZXx6rpL+v+O3NY/iNRZn7p3B7hyiNfTm6Dcl2IQv2Jy514v7QrL5UwakrXDbJmv2\ngrA09mAvpuOt8MM0luPd3fjapr9SqArGr5ZtpEzqu0mPMBOGS/RSOvXBGBGWQqM1XsPxZI6mCEN4\nbGJV2oebmFWGodRkP7wLl4Si6D1m4A1uSFy2dujQoUOHf4zv2otPzxZAFD8AAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$e^{- \\beta r^{2} - r}$$"
      ],
      "text/plain": [
       "      2    \n",
       " - β⋅r  - r\n",
       "ℯ          "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "beta = Symbol('beta')\n",
    "R_T = exp(-r - beta*r*r)\n",
    "R_T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Hamiltonian for this system is\n",
    "$$-\\frac{1}{2} \\nabla^2 - \\frac{1}{r}$$\n",
    "The first term is the kinetic energy of the electron, and the second term is the Coulomb attraction between the electron and proton.\n",
    "\n",
    "The first step is to compute the derivative of the trial wavefunction in spherical coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def del_spherical(e, r):\n",
    "    \"\"\"Compute Laplacian for expression e with respect to symbol r.\n",
    "        Currently works only with radial dependence\"\"\"\n",
    "    t1 = r*r*diff(e, r)\n",
    "    t2 = diff(t1, r)/(r*r)\n",
    "    return simplify(t2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "del1 = del_spherical(R_T, r)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construct $\\psi H \\psi$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "H = -1/S(2) * R_T * del1  - R_T*R_T/r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUMAAAAcCAYAAAD7nN2RAAAABHNCSVQICAgIfAhkiAAACRRJREFU\neJztnHuwVVUdxz9X7hVUBEcEgcJAyJRRQAhCFAuGEh9RGSBImFL2VDM1LaHmFFiW5sTUqDnJHKNC\nw/H9fpuVSFOmWc2USYWVaXpVRMRS+uO795x19tlrv85+3Mtdn5k75+611/7t9d6/9fv99gaHw9HX\nGVN1ASLYF9i9jBvtUsZNHA5Hj2U0MKPqQkTwPHB21YUAWAEcUnUh2uRdwOeBGnA3cGTFchw7LyOB\nS4HOCu6ddXx+E+iwnJsLLAVOAC4Cdm2viLHsCVwL7BdIn+aVozK+AHykygLkwEDgQuN4IfAq8JaK\n5Dh2fuYB3y/5nlnH5wTgDMu5o4FRxvFNwKEZynZawnwfRwv5DqStBrkq5X2PT5nfyqHAHXkJq5AJ\nwJvAWO94EGrshRXJcfQN1gLzS7xf1vFp2/ntD0wxjscDvwC6MpStljK/bTG8BBiXUMaBwA0p72tV\n579N+krkyTmo4b8Rk28/4FSgH7AHehJ+BnjWO/874DDgKe/Yf9L9OWV58pKTJ3m1UdUcgMbaNmC7\n91sDtkRcU3Wddgc2oEUojIvRZLweeKOE8kSNz7eihfnDwDPAH5A2OBKYSvj4mQKsB74IDANOAt6J\nnBk2WdtyrlOQx70yPJkg74nAukDaNFSf3VCbbAS+AvwjStDEhDcsircBW4lfjGcC5wP9jbTLgZ9G\nXLMWLfTtkpecrBTZRmUyHvg7DQP+cGAT0fWquk5TgV8hDSaKnwEfKr44oZjj8yjkKH0UWOL973uP\n77Zcv8D7XQosB37tyYuSZaOWsuw2zXAeUgCS8Bu06PlMBu4C9vKOB6L+eda8V5hmuBB7I5XBcuJd\n6W9HT8Kvh5yzOXw+BvwLOC970XKV0w5FtVFWJiAt4X8prulEmtMlwC+9tC40ULst15RZpyAHIY3v\nOZJpe/cCx6I6ZuV0GlvfMB6hVQMKjs87gb2RtrgObac3eef6hcgcDvzN+3+t9/sScshEyQJpjmfS\n7JA5AhhgHG8BLoiok41tJHPgTAN+T7OmuhLtHF70jl8BPocWzQuBRTZhjyBjZhUcjwq2g+gnyrdo\nDQvqAP4EXBmS/zg0SEAdMzpj+fKS0w5FtVE71EnfFsuA12k8rZNQZp2iqBOvGb6X5oWiDGzjcz5w\nXUj+u0LS5tPaxmuAT8TIslFLkRfsmuF84FMJrl+NnD8mr6AdyLBAejfwH/8gWOlOYBJaMctmIHAM\ncHVMvgPRfv/NQPpSpBoHNbZ3oyfWreipNxcYkaF8eclph6LaqAoWI3PMi3EZPXpDnUweRZO6rIiD\nqPE5B2mqQf6NwllMjqA5tGUCsrH9IEZW0YygYRO10Q+YRevOdhNqmz0C6dsxttPBbfIopIpWYVz/\nEvHOAJDx1reHrEFbsymowyZjrPTIK3YLWkRMBgfkzUJ1X4zsJeM9eavRBEwipwzybqMkdS+CDjTp\nNqLt1/tQ244GvooWkiBp+r0n0I20nJGEG+knoa1bNzIP7AN8lmjHkY248XkA0qqDPIhsoPcZaRvQ\n7mMLUpaGIHud/xCyyWqXJWhMgLauPwe+Z5yfjPo9itnAQ7Saa6ajRf8ZI20kWiAf8BOCi+FQ7/fl\nmJvmzSSkyv4lQd4BwGve/88jQ/rjwDuQGv1lI+9TtD75TPqjYNXTkOdtPeqIW5B23I0mbJycMsi7\njZLWvQiGeGUchhbfFV76e5Bhezqy+5ik6feewBtoQRkacu5kFCB9DI2F8lw0mW/McK+48Tnbkn4d\nervDXwzHerKidh42WTaSepl/7P19OuTcADT2t8bIOJGGBmuyNeTaM9ACv9wmbCZ6moUZVotiFxRQ\naRpHbfawQcBZFjmX0bzyJ+Eo5PHbFU20VV76KOC3aHL2BIpoozzrXiedzXBfVP7XaPb6ATyNFmST\nvPu9XerE2wwBNiONx+Qw4L/A4UbaNLS9Ddq0yuBsGgv2Eqp5eyaOZcDBMXkGoHFre5vGZBxaXFeZ\nicGKb/d+9yTaljMRDYgkNwZte06xnPukJ+v1BHJmIs0hjC1oknWS3Kv5EtoWTEeahh+esRlpYu3Q\n09uoyLrH8YL3+yStmsPTyPnQn8Z4zNrvefZBFsw6+NSQ1v1B9DDqB/wRmSnK3pEBfAc5Xa5A5U0T\nEVAGo9BO4ImYfMcCtxP/kOoP/ATVd0VUxjGesLK+YjEc+G5Iuk3rWYH94xK30QgJSEsN2ZySTpoy\nKbqNaiSv+1Xo6Rv8ewGF1oSdmxIqSXbph0LSH0R1M51TRfV7Vuok0wy30axhd6HFpuzX9foC1xKv\nPXaghfBrYSeDmuFmpH0MoZywgDm0vjrjv/KzCGkoP6Thyh9HqzcRtGDMQXFgWZhFYxL2NIpuozR1\n/6glvY4W1b8mkOGzAcUNBvG1qeeMtKL6vUiGoq2baeMdgjTBJHZfR3IGIw94nPa4EmnhK420k9D8\nCWUj1cUZgmxPYVrPYDRBwmwaF6A4ojQxaz67ocl3eoZrq2I0+bRRXnWvkz7OcDH6mIBpM+xA5hnz\nbZKi+r0d6sQ/PA6n1ZbpmyfOD8k/lnDngSOeZehVuyhOIVwjvML/J2zrcQ/Vft+sK/DrMxN9RmgV\nzVu6D6BYs+NIHrNmMgM5ER7IcG1V5NVGVdb9GuAxFE7iswA5F8410orq93bw36iIegtoBppLJv4W\neS7NdTkSecOtGoojkkW0voljMhuNof2BHxl/VxMTJTIRbZHLtp8NAu4H/knD0/gwMjSDVNsuFF1+\nJbKjrUHxcGHhC0lZirThnmgvDJJ3G+VV9zrZ3sbZ27t2PbLlrKP19bOi+j0tw9DraE+gtt+BNNb7\nafUY46XPC0nvRPGS16AJuho5MNyHlrMxHJl5ovBjPsP+VkZcByju6P3tlTF3Lqq6AL2AqtqoTnGv\nJvbGfh+DYjfdAlc8Z6Lg9baxddZZ3k16CnuRLTK/L1FlG71MIyA6T3prv9fQx5HDnD6OfFlACV8s\nOoecVtwcmIe8hg47O2Mb9cY6zUaf/ncUzzgUrJ4LUWr8xcgmMzGvm7XBVGQbc9jZGduot9VpBLLf\n2j6l78iXsI+4OhwOR5/jZlo/TuFwOBwOh8PhcDgcDofD4XA4HLnwfzLyQErtnCoHAAAAAElFTkSu\nQmCC\n",
      "text/latex": [
       "$$\\frac{1}{2} \\left(- 4 \\beta^{2} r^{2} - 4 \\beta r + 6 \\beta - 1\\right) e^{- 2 r \\left(\\beta r + 1\\right)}$$"
      ],
      "text/plain": [
       "⎛     2  2                  ⎞  -2⋅r⋅(β⋅r + 1)\n",
       "⎝- 4⋅β ⋅r  - 4⋅β⋅r + 6⋅β - 1⎠⋅ℯ              \n",
       "─────────────────────────────────────────────\n",
       "                      2                      "
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "simplify(H)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The integration occurs in 3D over the electron coordinates.  Because the integrand only has a dependence on $r$, it can be converted to spherical coordinates, and reduced to a 1D integral over $r$. (There should be an additional factor of $4 \\pi$, but it will cancel since it occurs in the numerator and denominator)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "h1 = simplify(r*r*H)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Substitute a concrete value for $\\beta$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "h2 = h1.subs(beta, 0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perform the integral"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMQAAAASCAYAAAAEwXFZAAAABHNCSVQICAgIfAhkiAAABZ5JREFU\naIHt2lusXVUVBuDvYGspPT2ACikEApSbJijUEJCUFEQSQ+QBGwgvDeCLQAwiERQxXAKt3BIIF6kK\nSIxGJQgpARKg4RJekKpcIiF4jQpKNaaAUWmBtj6MuTzzLNY6Z6591mH3Yf/Jzt5rzH+POec/1p5r\nzDE3I4wwwv8x1qOvQ3E1/ozt+Aguwt/nyNdJ+A5ewFvYgm1Z+y/w7Rb/F2J+6iPHPvhG6nMBdsE1\n+HWNtxSXYzO2YhG+jo0D8o7CxViYxrABl+GvNV6pxsOYx1zEo1SXUt6ncHYa30Khy7fSmHvFrngV\nqzLbJXgRH5wjX18VAW97ndTifz/8B1fU7B/GOuyV2Q7Ab3BwzfZPrMxsq4So8wbgfRKPYrd0PY6n\n8A/sn/FKdRnWPPqOR6kupbxleBA7Z7a1+BeOaBnbwFiTBpAL9CG8g3PnyNda7CtWlp0y+3LcOo3/\n74kAXVGzX4TzGviX4Nrseh021fpcKFbPswbgPYSDan0uS2P8aWYr1WVY8+g7HqW6lPJuTLbTM9vJ\nyXZzZcgHPhuchmfwbmbbhJdT21z42opXxA1RPZrHcSm+1uJ7JR5vaVuKExrsm0WQiZX4c/i9qenA\nWyKNObUjD45LY9ozsz2HN3BiZivVZVjz6DsepbqU8p4TT4NNmW08vf+3ZQwDYbH4lTXlh4/gzffR\n11oc3dI2jjvS56YV6exk/zF2T7YFeBafSNdLEuepBv8bxArehUfk9VtEapJjo0gl6KbLsObRhNnE\no0SXLrwmXCcWmMNn4HXCYWJC9Q0R3GdyYzfXvpaLx28b1uDA9LkpAAvExm87XhN58o9MzX0/IETe\n0OD/1fTdeR14xAZ1SY2zd+I8ka676DKsedQx23iU6NKFV8cBohjxxdzYR8o0kd7fbmirfqG7NbT1\n7etmzTcMsWn6N/4wTd9bRKrxiBD4h6mvZzPOVtwt0pK8QrdEBIHI60t5xLzqVZ0vi7Tjm+m6iy7D\nmkcds41HiS5deBVOxlW4HzeZfEr1hqM1/8KJTc12k+LNla/P4KUWnzvhB6ZWYtr6uBi3i7z5jyZX\n2Y9nnD3wW5GaECvkanHDbRdVni68Og4SN8vqzNZVl2HPo6945GjSZTa8eViPn4vy9XtwuNh4PF/4\nuit9b6n2CT2Q2hbPMLgKg/q6F99v8XkuPl2zNfVxgahYVFgkKhPbvLdOvbuoy98ghN8PvxSbzbEB\neBUWiPTkhpq9iy47wjz6iEeONl0G5VVYkfr+WSG/CIuE2Dc1tD2O1+fY13yxIqxpaFuCWxrs9QCM\nJd+HNXDPSfyPtQ064RU8NgNnOt6Y2Ahf2dBWqsuOMI8+4pFjOl268D5qsqhQYSL1vc1kxakX/Eps\n7ur4i+k3Nn34Wi4mdUFD2yrxWFyXvR5K/JfT9UpRsmvbsI+JKs6yaca8R/r+OdNwZuKtFiXKHGdk\nn0t02RHm0Uc8csykSwlvQjzN3jW5kScWmurgcNcGnwPjSnFMnj8+D0wdfanGPVgc7PThC85MbaUH\ngPtrfkJsxLEN/MWiolLdZF8R5cZ9Ms6FiZOfgpby4AuaV7a8SlOiy7DnQT/xqFCiSwlvZ3E+8jtT\niwBHpr6fqQxtJbOuWCt29qtEVYM4LX1JbO4qrMCTYpX47Cx9VagOZJoqME2YX3snRDkf3xUb0T8l\n+wTuFDfFlmQbFwc51QHZMnE6/Hlx+KUj7wRcj4dFebTCvDSuCiW6DHMeFfqIB+W6lPA2J85rpqbd\n54v0ruu/KYpwhHj83ShKWfeJo/wch4ra7209+KpwijiVPHIGnxMitfibEGoznk7fr7BCPLbvFtWb\ne3B8zc9CIe5dYjO2Hsc09FfKe137/3+uqnFLdRnGPCr0FY9SXbrodxZ+IhaU9UKfQ2YY5wgjjDDC\nCCOMMMIIU/E/53G1N9/WAb8AAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$-0.0748992089974223$$"
      ],
      "text/plain": [
       "-0.0748992089974223"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "num = integrate(h2, (r, 0, oo)).evalf()\n",
    "num"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also construct and integrate the denominator (the normalization)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "norm1 = r*r*R_T*R_T\n",
    "norm2 = norm1.subs(beta, 0.1)\n",
    "norm3 = simplify(norm2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAKsAAAASCAYAAAAkLS6kAAAABHNCSVQICAgIfAhkiAAABjJJREFU\naIHt2muMXVUVB/DftB1LoS1QWjo2kBYL4iPRlhTQoo1vQ8IH0w+KseLrA6IxEoT4COjgABJMwBCw\npjE6FCNFEEmQRMQIqQlUfKFWCajRGlSwxipapRYcP6x9nX337HPn3NvbD8b5Jzf3nrX3Wf+119l7\n7bXXucxhDv8jGCmuT8WnsAdTWI5L8OQhcByJXXhJjz7Pxzj+iQPpexx/69O2s/E5/DjT9e+s/Xu4\nMf1+Hj6Bp/EsjsKH8URh26OYwH3Yjw34CD6ARwbgfRnOT/0WCf9cle4dxC9nJHsW4QQ8hI/jdxV9\nOS7GqPBpjmHrazvevniPxuPYksk+ht14ziyGNuF08aCmevR5EX6Ljel6DL8WD6Vf2z6UuJo+Z6d+\nJ+FP2Jzdu0U4cEFhX6njX3h/0act73p8HUdk927FU1hX6Gzjl9PwTRyTrhdjJ/6INZqxWiy88UI+\nbH1tx9s375WpMX9Yy3AQF/QwtIYX4m5MiqjaNFkXiMh1YSY7EXvxwQFs25ruH8W8TH4Wbsiu78Sf\niz6LRJR9Z2HjHmzD13BNGluJtrzXCV+8JZOdk2TXZ7K2frkbJxe2rE/6dlTs7GBb6jNeyIetr+14\n++Z9DHdV5D/Ft3sYOhsmNU/Wd4tIdUxDewdtbbuh0mcxviG2HyISHxTbTIlHRSTIcf8strXlhfPw\nV7w+k50r/HN1Jmvrl7+L6Ht8Id8ndo4aNmec44dZX9vx9sW7JCm4sWzAPYlwUExqnqz34uez3H+o\ntm3Fmdn1WNK3s9L3IRHBc9w/i/62vE24Bs/gpZmsjV+IxXpApDU5nhDbconF+Hz6XZtcw9ZXQ228\nrXg72+rq9P1URfl+LMXCpHBYGMErxATZhDeIwa/B5fjREGw7C/Px3Uy2F//QnUd1sAorhF+eSbKF\nuAzHiYi8VhwEHusxthpvDSfhHXif6QNHW78Qh5clug+Fq7BSfZF91MwDUI5h6ytRG2/fvBvFyri8\nQnBzalvZh1E5JtUj6/IkfwTvzeSvEqfdFw/Bth+YuVrhC2J7yashY+IEP6V7O/qlyBc72CIqEGMN\nnL14OzhHVBh+Ig6KeZ7b1i9NuFpUNzYW8nVicnXQNhIOQ1+v8fbL68weZDtS26oWBDVMqk/WlUn+\ntDjc5HjcdO44qG2v1byVrhCR8fx0vQBX4IdJ33FZ39Kx80Vkvl4dvXhLLBBb/i4xSWnvlxpOFvnf\nFYV8Hm7SXTlpM1mHra823n54ETXHJrK7UtuSWQxpwqT6ZB1N8t2Vtl1iW194CLZ9VUTQJhwr6qzX\nCqesxvdFPbCsP5fYkz41zMZbYpMYw+3puq1fSiwUqcO1lbYL8OpCNtvkGra+DsrxtubtRI0nk4Jj\nKzcfhb/oLkQPAwdF/riv0nZArNplA9o2ijfiDz3494nU4iJcKibfSjxgenHtxHcq985XTz1m432B\nmS9HHk7fm0Vu2tYvOUbwRVF9uKhoGxM12/sabKphWPrajLct738PWPtF4n5i2UGE5Icr8mFgF06p\nyDsHpr3ioNOvbWeIidxUbqlhhXhzcmUmWy9y1hLL1SNrL96lYhyj4m3cr5L82fQ9IhYB7fySY0Lk\nuBOZ7Dxsx+vEpLkzaxtN3+eK3HM77hiyvm/1Md42vF34pHi1lW+Ba0WUKd/YnGJmPtWESc2lq7eK\n/C/XNSKi5VcGtI04cU5pfplxoShRnZDJLhY5YV4luE13/sp0sfqyPnmPEFHzF7oj44Z0T145aOsX\neJfwT4ltFVkHazRv28PS1894++Z9rnDG2zPZZ/Az3Yn0JnFqvqeH8Tk6h6AjK23z8KCYKB28WUSO\nNQPY1sElifM9DTZdit+YPtGvF+lGefI8XeRWHdtHxCp/QD1vnI33KvGfgnzR3SzSmNMyWVu/vEZE\n8S8Vnx24pcEGIthM6d5FDoe+tuMdiHedePV1nSj23mHm9nuqeLCf7WH88WIy7zb9fnyvyHPeVvRd\nJqLvbfhyMm7tgLZ18CYxuTc0tC/Cp0V+dLs4ob68oe8rk23bxdY3oV6jbcNLvM69RTy0e3Gr+MNK\niTZ+2af5/wgTZmKpeAa/N11xeDDZfTj0tR1vv7xzmMMc5jCHOczh/wv/AdXUUbu5/Mk5AAAAAElF\nTkSuQmCC\n",
      "text/latex": [
       "$$0.160795736242432$$"
      ],
      "text/plain": [
       "0.160795736242432"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "denom = integrate(norm3, (r, 0, oo)).evalf()\n",
    "simplify(denom).evalf()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAALsAAAASCAYAAAADg69MAAAABHNCSVQICAgIfAhkiAAABtdJREFU\naIHt2m2MXkUVB/Bf2y19ZYWWQiEQWrBgTZAWC1Wxje+KEsFGRWIVE18aUSNGqoCgVSy+JYBErYoJ\nazEWrKn9oB9ME00wEdREC6KAoAhUpCwpAlZaylI/nLns7N377DP32f1i3H/y5Oaec+658585c+bM\n3IdJTOL/BFMm0NfJ+BIewEEcgfXYPQG+L8b05L+Ok7ABT2N/um7AU5nNPbgSv8RerMAl+BjuyuyO\nxaWp/TMwG1/GH2vvLOX6MqxLbZqV/F2F23vkW8qjwmzchpd0eM8J+Bz2YQhz8Gk80iOPXmNgrPHt\nZncWvp3aUsXAc5n+d/hmF7+t8ALswtpMdhnuxCHj9H28GNgNDboX40G8It0vxP0Ntgdrv2fwkZrN\nfGzH0ZlssQiwJZmslOty/BQzM9kmPIllDVwqjMW3hEeF08VAH+ygX4zHsCaTrRVB05fJSnn0GgNj\n8S2x+6TR/ZL/zuritzU24lEjO2keDuDD4/T9XdHoDTV5nwjEizLZcRjEx2u2DyQ/P8FXsbThPetF\nhqzjMnwluy/lek1q93mZ7Owku67hPRU68aWMx1L8DAMiq3cK9u3Yg6mZbJbI8u/rgUevMTAW3xK7\nTWLcpxvJ5Ux8IzfMlePBO/AbPJvJ9uDupOsVa/CLDrr3iuw0kMkewgJ8vWZ7Pz6Et+FTmpf8E/Ca\nBvk+0ZEVSrn+QWS/PZlsbrr+p+E9jM2XMh534S0iYO/u4OeQZHOfkUv+02JCvT2TlfLoJQa68S2x\nGxLjfsAwl7m4QvTR85iIYD9ULPMPNugexkt79DsXb8ZNHfTni8H6V4/+69iJc/FDHJ5kM8SyPJDu\n23DdLJb2HZnsNDE4Wxqe78Z3IjFPZOB9DboncEZ2X8Kjlxgo5dvN7qMNsq+JvciIpNLXYNgWx6fr\nkw26vegXQbO/pd9Ldd6wTMEr8VusxhtEpyzC50U2yjFDzPT5IgOcKDZ2f8lsBvABMYleLcqaN+Ez\nuCPZjIfrYlyACzVvUMfi24ZHCQZFIMxs0B0jVsc+I7N0hSYevfRLCd82dhXOxDSxyozARAR7f7o+\n06Dbm66HaXcqswz/xl876OeLgTpSbFIvT/JX4RZxevCnzH6BCOaH0v1a/AqnGj552C/KmK14I24U\nte/vMz+9cD0bK3GOKK++1/BsN75teJRgCDfjrSJxVHX9QhHsRPZ/NHtmLB5t+6WUb6ldjuuMLMMm\nFCt13mDclHTHNOg6YSq+b+QOvu7/qCTbJzZVOXaJ04O6zxzTRGarbxQvwfWinv1besc/cUrSj4dr\nnygFbhNHcnnbuvFty6PCgM4b1AViRViXte+LYnIfFAmlCU082vRLKd82/VLhtfhzJ2XeeaeK5X9n\n4e+G9NzgGC+fk65PjWFTxzoxSE1ZokK1WbpPbKpy7MLrxbJZ4bmazZBo9zmZ7BNYhQ+KjH4KrhUT\n6wfJZjxcnxVn5CvFuXCFEr4VSniUYjC1ZSGuFkF0fXrHPiM3pDmaeLTpl1K+bfqlwoViEjYiL2Nu\nF2eqbbFbzLjDG3RzxAayNNgXirJkUxe7A6KDH2/Q7RfZYJ7IyreIpXpVzW6a4cw0BZ+t2ewVE+Ce\n1J6lYgNWyvVFqR13ZDY703WN2GPMVcZXIY+2eFzscXIchV8bXhFKeJTGQOn4ltrlmC7Kz/pJ3POY\niJp9r1gRjmvQvdBwx5TgdaJzt2ey6tjvXaKG24xtYgYvMRrVRqjKNsvFClDHEeKYjVjSD8O9DXbf\nEefsM5Vz7U9208VXxarmHErXKSJI2/At4TFeLBBfkTem+1IeTyjrl1K+swvttmX6M8TEeqwry3Hi\nC/iHkX8/OFHM9voXviVG19ljYZHmWu18Ua/mvqaILPKjTLbV6PpzefJ5RfbcI+KEp45DRWlUlUUl\nXGeK1edescJUWJHsRp0UZFikmW8JjzoGdK7ZLxIb0GMz2cWCa3VK04ZHmxjIsUjZR6Vudhck/Xg/\nYnbF0SLI3pPJrhUnIvkGY7WoCX/ewvcSQWJjTT4Vt4oBqvBOkdEXZbLT8WORLYjB2CyW6ryuPy+1\nN3+2XwRZvrsv5XqV+CKbD/6NYjk/TWd04lvKI0e1OZzdoLscfxclAzFxdhv+60VbHqX9Ukcnvm3t\n1if9+7v4mRAsExu7a8Sx1Dajl7WTRYd+q8Bfv/jD08OGT15uFR9+KswT2Wur+Bi0RWSTOlYlm81i\nabxS8xnz6qS/WQTKVnGcWUcJV+Ir5hYRHDuS35Ma7CjjW8LjSJFM7jT8/5DB5Pvdmd0s8fHlBjGJ\nduDlHdpWyqO0X0r5trE7V0y2FR3eN4lJTGISk5jEJCbxv4r/AmY6RtXOL3TYAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$-0.465803451930447$$"
      ],
      "text/plain": [
       "-0.465803451930447"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E = num/denom\n",
    "E"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And, as expected, energy is greater than the exact ground state energy of -0.5 Hartree.\n",
    "\n",
    "## Find the minimum energy\n",
    "Collect all the steps for computing the energy into a single function.   Even though this particular integral could be done symbolically, use numerical integration instead."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def compute_energy(R_T, beta_val):\n",
    "    \"\"\"Energy given a value for beta\"\"\"\n",
    "    \n",
    "    # Normalization integrand (denominator)\n",
    "    norm1 = r*r*R_T*R_T\n",
    "    norm2 = norm1.subs(beta, beta_val)\n",
    "    norm3 = simplify(norm2)\n",
    "\n",
    "    # Integrand for the numerator\n",
    "    del1 = del_spherical(R_T, r)\n",
    "    # Construct psi * H * psi\n",
    "    H = -1/S(2) * R_T * del1  - R_T*R_T/r\n",
    "    h1 = simplify(r*r*H)\n",
    "    h2 = h1.subs(beta, beta_val)\n",
    "    \n",
    "    lim = 20.0\n",
    "    \n",
    "    denom_func = lambdify([r], norm3)\n",
    "    denom_res = quad(denom_func, 0.0, lim)\n",
    "    \n",
    "    num_func = lambdify([r], h2)\n",
    "    num_res = quad(num_func, 0.0, lim)\n",
    "\n",
    "    e = num_res[0]/denom_res[0]\n",
    "\n",
    "    return e\n",
    "   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the energy can be computed vs. $\\beta$, and we can find the minimum energy.  In this case, the minimum occurs at $\\beta = 0$, which we know is the exact wavefunction for the hydrogen atom."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fa4bc2c4910>]"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEACAYAAAByG0uxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYFNW5x/Hvi0ZFRcQoiI6guAUTNSJucaEVcUEFtyAm\ncc01GJeQXMVAghdMclUQxSUuj0B0XCISBB0BZRFacc9lQDRwCWpcgjBRFBUIFwbe+8epkWboWZjq\nmeru+X2eZ56pqT7V/U4zMz/OOVWnzN0RERFpiBZJFyAiIoVLISIiIg2mEBERkQZTiIiISIMpRERE\npMEUIiIi0mCxQsTM2pjZNDNbZGZTzax1LW1bmFm5mZVl7Bsb7Ss3s3+YWXnGY4PMbLGZLTSzU+LU\nKSIijSNuT2QgMMPdDwRmAoNqadsfWJC5w937unsXd+8CPAVMADCzzkAfoDNwOnCfmVnMWkVEJMfi\nhkhvoDTaLgXOztbIzEqAnsDoWp6rD/DnjOcd6+6V7v4BsBg4MmatIiKSY3FDpK27VwC4+zKgbQ3t\nRgIDgKyXx5vZ8cAyd38/2rUn8HFGkyXRPhERySNb19XAzKYD7TJ3EcJgcJbmm4WEmZ0BVLj7PDNL\nRcdXdyHwRH0KFhGR/FFniLh7j5oeM7MKM2vn7hVmtjvwryzNjgV6mVlPoCXQyswecfeLo+fYCjgX\n6JJxzBJgr4yvS6J92WrQ4l8iIg3g7rHnmuMOZ5UBl0bblwDPVG/g7r9x9w7u3gnoC8ysCpBID2Ch\nu39S7Xn7mtk2ZrYPsB/wZk1FuHvefQwZMiTxGlSTamqOdamm+n3kStwQGQb0MLNFQHfgVgAza29m\nk+r5HBdQbSjL3RcA4whnc00BrvJcftciIpITdQ5n1cbdPwdOzrJ/KXBmlv0vAi9W23dZDc99C3BL\nnPpERKRx6Yr1RpJKpZIuYTOqqX5UU/3lY12qqWlZoY8SmZlGukREtpCZ4XkwsS4iIs2YQkREpJlZ\nty53z6UQERFpRior4Uc/yt3zKURERJqJykq46CJYuTJ3z6kQERFpBtavh0svheXLYeLE3D1vrOtE\nREQk/61fD5dfDkuXwqRJsN12uXtuhYiISBHbsAGuuAI++ggmT4aWLXP7/AoREZEitWED9OsH774L\nzz0H22+f+9dQiIiIFCF3uPpqWLgwBMgOOzTO6yhERESKjDv84hcwbx5MnQqtWjXeaylERESKiDv8\n6lfwxhswfTrstFPjvp5CRESkSLjD9dfDyy/DjBnQunXjv6ZCRESkCLjDwIEwa1YIkJ13bprXVYiI\niBQ4dxg8GJ5/HmbOhF12abrXVoiIiBS4oUOhrCwEyLe/3bSvrRARESlgv/sdjB8fhrF2263pX18h\nIiJSoG6+GZ54AtJpaNs2mRoUIiIiBWj4cCgtDQHSrl1ydShEREQKzB13wKhRIUDat0+2FoWIiEgB\nuesuuPfeECB77pl0NQoREZGCce+9IUTSadhrr6SrCRQiIiIF4IEH4LbbQoB06JB0NRspRERE8tzo\n0eFMrFmzYO+9k65mUwoREZE89tBDcNNNIUD23TfpajanEBERyVOPPAI33hiuRN9vv6SryU4hIiKS\nhx5/HAYNghdegAMOSLqamrWIc7CZtTGzaWa2yMymmlmNCw+bWQszKzezsox9Y6N95Wb2DzMrj/Z3\nNLPVGY/dF6dOEZFC8uSTYUn3adPgO99Jupraxe2JDARmuPtwM/s1MCjal01/YAHwzS1S3L1v1baZ\njQBWZLR/1927xKxPRKSgjB8Pv/xlCJDvfjfpauoWqycC9AZKo+1S4OxsjcysBOgJjK7lufoAT2Qe\nFrM2EZGCMnEiXHNNWNL94IOTrqZ+4oZIW3evAHD3ZUBNS4CNBAYAnu1BMzseWObu72Xs3jsayppl\nZsfFrFNEJK+VlcGVV8Jzz8GhhyZdTf3VOZxlZtOBzOW9jBAGg7M03ywkzOwMoMLd55lZiuw9jAvZ\ntBfyCdDB3b8wsy7A02Z2kLuvrKteEZFCM2kSXHEFTJ4Mhx2WdDVbps4QcfceNT1mZhVm1s7dK8xs\nd+BfWZodC/Qys55AS6CVmT3i7hdHz7EVcC7wzfyHu68Dvoi2y83sPeAAoDxbHUOHDv1mO5VKkUql\n6vq2RETywvPPw+WXhyDp2rXxXiedTpNOp3P+vOaedYSpfgebDQM+d/dh0cR6G3evaWIdM+sGXOfu\nvTL2nQb82t1PzNi3a/S8G8ysE/AicLC7r8jynB7nexARScr06fDjH8Mzz8AxxzTta5sZ7h577jnu\nnMgwoIeZLQK6A7dGxbU3s0n1fI4L2HQoC+AEYH50yu84oF+2ABERKVQzZ4YAmTix6QMkl2L1RPKB\neiIiUmjSaejTJ5zOe8IJydSQLz0RERHZArNnhwAZNy65AMklhYiISBN55RU477xwX/RiOf9HISIi\n0gRefx3OOQceewy6d0+6mtxRiIiINLI334TevaG0FE45JelqckshIiLSiObMgbPOgjFj4PTTk64m\n9xQiIiKN5JVXoGdPGDUKzjwz6Woah0JERKQRTJ8e5kAefRR69aq7faFSiIiI5Ngzz4QLCSdMKL45\nkOoUIiIiOfTnP29cjfe4ZrD+uEJERCRHHnwQbrgBZsyAww9PupqmoXusi4jkwO23wx//CC++CPvu\nm3Q1TUchIiISgzvcdBOMHRuWNCkpSbqipqUQERFpIHe4/np44QV46SVoW9O9XYuYQkREpAHWrw8T\n6G+/DbNmQZs2SVeUDIWIiMgWWrcOLrkEli0L14O0apV0RclRiIiIbIE1a+CCC0JPZPJkaNky6YqS\npVN8RUTqaeXKsHzJdtuFCwmbe4CAQkREpF5WrIBTT4WOHcMFhdtsk3RF+UEhIiJSh08/hRNPhK5d\nw2KKW22VdEX5QyEiIlKLJUvCbWzPPBPuvBNa6K/mJvR2iIjU4P334fjj4bLL4Pe/B7OkK8o/OjtL\nRCSLhQvDCry/+Q38/OdJV5O/FCIiItXMnRtuJjV8OFx0UdLV5DeFiIhIhldfDTeTuv9+OPfcpKvJ\nfwoREZHIjBlw4YXw2GPhdF6pmybWRUQIdyP80Y/gqacUIFtCISIizd4TT0C/fmEZkxNOSLqawqIQ\nEZFmbdSosJz7jBlwxBFJV1N4NCciIs3WHXfA3XdDOg377590NYUpVk/EzNqY2TQzW2RmU82sdS1t\nW5hZuZmVZew7wszeNLO50eeuGY8NMrPFZrbQzE6JU6eISKaquxE+8EC4mZQCpOHiDmcNBGa4+4HA\nTGBQLW37Awuq7RsODHb3w4AhwG0AZnYQ0AfoDJwO3Gema0VFJD53GDAgTKDPng0dOiRdUWGLGyK9\ngdJouxQ4O1sjMysBegKjqz20FKjqvewMLIm2ewFj3b3S3T8AFgNHxqxVRJq59evDBPrs2WEIq127\npCsqfHHnRNq6ewWAuy8zs5ruMDwSGMDGwKgyEHjFzG4HDPhBtH9P4LWMdkuifSIiDbJuHVx6KXzy\nSZhEb853I8ylOkPEzKYDmXltgAODszT3LMefAVS4+zwzS0XHVxkDXOvuT5vZ+cCfgB71Lz8YOnTo\nN9upVIpUKrWlTyEiRWzNGujbF9auhSlTmufNpNLpNOl0OufPa+6b/d2v/8FmC4GUu1eY2e7ALHfv\nXK3NzcBPgEqgJdAKmODuF5vZV+6+U0bbFe6+s5kNBNzdh0X7nweGuPsbWWrwON+DiBS3Vavg7LNh\n553h8cd1M6kqZoa7x55rjjsnUgZcGm1fAjxTvYG7/8bdO7h7J6AvMNPdL44eXmxm3QDMrDth7qPq\nefua2TZmtg+wH/BmzFpFpJlZsSKsxFtSEi4oVIDkXtw5kWHAODO7HPiQcEYVZtYeGOXuZ9ZxfD/g\nXjPbBlgD/AzA3ReY2TjC2VzrgKvU3RCRLfHpp2H5kmOPhbvu0s2kGkus4ax8oOEsEaluyRLo0SMM\nY/33f+tmUtnky3CWiEhe+cc/wt0IL7oIbr5ZAdLYFCIiUjQWLgwLKP7nf8Kg2i59lpxRiIhIUXj1\nVTjxxHAv9GuuSbqa5kMhIiIFb/x46N0bHnooXFAoTUer+IpIwXKHkSPh9tth6lTo0iXpipofhYiI\nFKT16+GXv4RZs+C117SQYlIUIiJScFavDrey/eorePnlcDW6JENzIiJSUCoqIJWCnXaC559XgCRN\nISIiBWPRIjjmGDjtNCgt1TIm+UDDWSJSEGbPhvPPh1tugcsvT7oaqaIQEZG89+STcO218NhjYUFF\nyR8KERHJW+4wfDj88Y8wfTocemjSFUl1ChERyUuVlaH38eqr4RTekpKkK5JsFCIikndWrtx4J8LZ\ns8OZWJKfdHaWiOSVpUuhWzdo2xYmT1aA5DuFiIjkjQULwim8Z58NY8bAt76VdEVSFw1niUheSKfh\nggtgxIhwLxApDOqJiEjiHnsM+vQJ90FXgBQW9UREJDHu4e6DDz4YFlL87neTrki2lEJERBKxbh1c\ndRXMmRNO4d1jj6QrkoZQiIhIk/v6a/jhD6FFC3jpJdhxx6QrkobSnIiINKklS8J90Dt2hLIyBUih\nU4iISJN5+234wQ/CWVgPPABbayyk4OmfUESaxAsvwIUXwl13hc9SHNQTEZFGV1oa7kQ4frwCpNio\nJyIijcYdfve7ECLpNHTunHRFkmsKERFpFGvXQr9+8M474RTedu2Srkgag0JERHLuyy/DXQhbtgw9\nkB12SLoiaSyx5kTMrI2ZTTOzRWY21cxa19K2hZmVm1lZxr4jzOxNM5sbfe4a7e9oZquj9uVmdl+c\nOkWk6Xz8MRx/PBx4IEycqAApdnEn1gcCM9z9QGAmMKiWtv2BBdX2DQcGu/thwBDgtozH3nX3LtHH\nVTHrFJEmMG9eOIX3kkvgnntgq62SrkgaW9wQ6Q2URtulwNnZGplZCdATGF3toaVAVe9lZ2BJ5mEx\naxORJjR1arj/+R13wHXXgek3uFmIOyfS1t0rANx9mZm1raHdSGAAGwOjykDgFTO7nRAaP8h4bG8z\nKwe+BG5095dj1ioijWTMGPjtb8Pw1bHHJl2NNKU6Q8TMpgOZ51UY4MDgLM09y/FnABXuPs/MUmza\nwxgDXOvuT5vZ+cCfgB6EHkoHd//CzLoAT5vZQe6+MluNQ4cO/WY7lUqRSqXq+rZEJAfc4cYbYezY\nsAbWAQckXZHUJJ1Ok06nc/685r7Z3/36H2y2EEi5e4WZ7Q7McvfO1drcDPwEqARaAq2ACe5+sZl9\n5e47ZbT90t03m5w3s1nAde5enuUxj/M9iEjDrF0LP/0pLF4Mzz4Lu+2WdEWyJcwMd4896Bh3TqQM\nuDTavgR4pnoDd/+Nu3dw905AX2Cmu18cPbzYzLoBmFl34O/R9q5m1iLa7gTsB7wfs1YRyZHPPoNT\nT4VVq2DmTAVIcxY3RIYBPcxsEdAduBXAzNqb2aR6HN8PGG5mc4E/AD+L9p8AzI/mRMYB/dx9Rcxa\nRSQH5syBrl3hqKPgL3+B7bdPuiJJUqzhrHyg4SyRpvPQQ/DrX8P998N55yVdjcSRq+EsXbEuInVa\nuxb69w+3sH3xRa2BJRspRESkVkuWhCVMdt8d3nwTdtqp7mOk+dBS8CJSo5degiOOgF694KmnFCCy\nOfVERGQz7nD33XDzzfDoo+FKdJFsFCIisolVq+BnP4MFC+D112GffZKuSPKZhrNE5BvvvQfHHBPu\nff7qqwoQqZtCREQAmDIlrMB75ZXw8MPhXiAiddFwlkgzt2ED/OEP8OCDMGGCFlCULaMQEWnGVqyA\niy4Kn//6V2jfPumKpNBoOEukmXrnnXD67j77wAsvKECkYRQiIs3Qk0/CiSfCkCHhVN5ttkm6IilU\nGs4SaUYqK8PaVxMnwvTp8P3vJ12RFDqFiEgz8a9/wQUXwLbbwv/8D+yyS9IVSTHQcJZIM/DGG2H5\n9uOOg8mTFSCSO+qJiBS5UaPC/c9HjYLevZOuRoqNQkSkSK1ZA9deG648nz0bDjww6YqkGGk4S6QI\nffwxnHBCuP7j9dcVINJ4FCIiRWbmTDjySOjTB8aNg1atkq5IipmGs0SKhDvcfnv4ePxxOOmkpCuS\n5kAhIlIEVq6Eyy+Hf/wjnInVoUPSFUlzoeEskQL397/DUUeFYavZsxUg0rQUIiIFrKwsXPvRvz+M\nHg3bbZd0RdLcaDhLpACtXw9Dh0JpKTz7bOiJiCRBISJSYD7/HH78Y/j3v8Py7e3aJV2RNGcazhIp\nIG+9FZZvP+igsICiAkSSpp6ISIF47DH41a/gnnugb9+kqxEJFCIieW7tWhgwICycOHMmHHxw0hWJ\nbKThLJE89s47YdL8gw/C8u0KEMk3sULEzNqY2TQzW2RmU82sdS1tW5hZuZmVZew7xMxeNbO3zOwZ\nM9sx47FBZrbYzBaa2Slx6hQpNBs2wB13hLsPXn01PP007Lxz0lWJbC5uT2QgMMPdDwRmAoNqadsf\nWFBt32jgBnc/FJgI3ABgZgcBfYDOwOnAfWZmMWsVKQgffgjdu8OECeHq8//4D9BPv+SruCHSGyiN\ntkuBs7M1MrMSoCchNDLt7+4vR9szgPOi7V7AWHevdPcPgMXAkTFrFclr7vDII+HmUaeeCi++CJ06\nJV2VSO3iTqy3dfcKAHdfZmZta2g3EhgAVB/u+puZ9XL3MkLPoyTavyfwWka7JdE+kaL02WfQr19Y\nwkT3PpdCUmeImNl0IPNsdAMcGJyluWc5/gygwt3nmVkqOr7K5cA9ZnYjUAasrX/pGw0dOvSb7VQq\nRSqVasjTiCRiyhS44gq48MKw+q6WLpHGkE6nSafTOX9ec9/s7379DzZbCKTcvcLMdgdmuXvnam1u\nBn4CVAItgVbABHe/uFq7/YFH3f1oMxsIuLsPix57Hhji7m9kqcHjfA8iSVm5Eq6/Hp5/Pixf0q1b\n0hVJc2JmuHvs2ba4cyJlwKXR9iXAM9UbuPtv3L2Du3cC+gIzqwLEzHaLPrcg9GweyHjevma2jZnt\nA+wHvBmzVpG88dprYchqzZpwFboCRApV3BAZBvQws0VAd+BWADNrb2aT6nH8hdGxC4Al7v4wgLsv\nAMZF+6cAV6m7IcVg7VoYPBjOOQduuw0efhha13hivEj+izWclQ80nCWFYsEC+MlPYI89wrLtu++e\ndEXSnOXLcJaI1GHDBrjzzjBk9fOfh6XbFSBSLLR2lkgj+ugjuOyyMPfx+uuw775JVySSW+qJiDQC\n97DqbteucPLJ8NJLChApTuqJiOTY8uVw5ZVhDmTqVDjssKQrEmk86omI5NBzz8Ehh0CHDjBnjgJE\nip96IiI5sGpVuHBwypQwjHXiiUlXJNI01BMRien110OPY/VqmD9fASLNi3oiIg20bh38/vfw4INw\n771w3nl1HyNSbBQiIg2wcCFcdBG0awdz50L79klXJJIMDWeJbIENG+Duu+GEE8LKu5MmKUCkeVNP\nRKSePv44XDi4alVYQHG//ZKuSCR56omI1MEd/vxnOPzwMGk+e7YCRKSKeiIitfj887De1dtvh2tA\nDj886YpE8ot6IiI1mDo1XDi4xx7hwkEFiMjm1BMRqWbVKrjhhrDabmkpdO+edEUi+Us9EZEMs2dD\nly7w1VfhwkEFiEjt1BMRIZx5dcMN8MorMHKkLhwUqS/1RKRZ+/e/4Q9/CPc733//cBGhAkSk/tQT\nkWbJHSZOhOuuCxPmc+bA3nsnXZVI4VGISLPzzjvQvz9UVMCYMXDSSUlXJFK4NJwlzcbnn8O114bQ\nOOccmDdPASISl0JEit769XD//dC5c9hesACuuQa2Vj9cJDb9GklRe+kl+MUvoHVrmDYNDj006YpE\niotCRIrSRx/BgAHhhlG33QY//CGYJV2VSPHRcJYUlX//G266KdxpsHPncMpunz4KEJHGop6IFAV3\neOqpcJ/zI46A8nLo2DHpqkSKn0JECt78+eGU3eXL4aGHdI9zkaak4SwpWMuXw9VXw8knhzmP8nIF\niEhTixUiZtbGzKaZ2SIzm2pmrWtp28LMys2sLGPfIWb2qpm9ZWbPmNmO0f6OZrY6al9uZvfFqVOK\nS2Ul3HdfmPOAMO9x1VU6ZVckCXF7IgOBGe5+IDATGFRL2/7Agmr7RgM3uPuhwETghozH3nX3LtHH\nVTHrlCKRTodVdv/yF5gxA+69F7797aSrEmm+4oZIb6A02i4Fzs7WyMxKgJ6E0Mi0v7u/HG3PADKX\nvtP5NPKNDz8MQ1aXXgr/9V8wc2a4YZSIJCtuiLR19woAd18GtK2h3UhgAODV9v/NzHpF232AkozH\n9o6GsmaZ2XEx65QCtXo1DBkSeh8HHxyGrs4/X6fsiuSLOkeRzWw60C5zFyEMBmdpXj0kMLMzgAp3\nn2dmKTbtYVwO3GNmNwJlwNpo/1Kgg7t/YWZdgKfN7CB3X1mP70mKgHsYshowAI4+GubOhQ4dkq5K\nRKqrM0TcvUdNj5lZhZm1c/cKM9sd+FeWZscCvcysJ9ASaGVmj7j7xe7+d+DU6Ln2B86IXnMtUaC4\ne7mZvQccAJRnq2Po0KHfbKdSKVKpVF3fluSxt94KS5V8+SU88gh065Z0RSKFL51Ok06nc/685r5Z\n56H+B5sNAz5392Fm9mugjbsPrKV9N+A6d+8Vfb2bu39qZi2Ah4BZ7v6wme0aPe8GM+sEvAgc7O4r\nsjynx/keJH989hnceCNMmBCuOr/iCthqq6SrEilOZoa7xx4YjjsnMgzoYWaLgO7ArVFx7c1sUj2O\nvzA6dgGwxN0fjvafAMw3s3JgHNAvW4BIcaishHvugYMOCqfpLlwIV16pABEpBLF6IvlAPZHCNnNm\nuNq8bVu46y743veSrkikechVT0SXZ0mTcw9LtA8fHu7tcfvt4SZROuNKpPAoRKTJVFaG+Y4RI2DF\ninB/8/HjoWXLpCsTkYZSiEijW7kS/vQnGDkSSkrgt7+Fs86CFlq5TaTgKUSk0SxdGibMH3wQUil4\n4olwzYeIFA/9X1BybsEC+OlPw9lWX30Fb7wRhq0UICLFRz0RyQn3sDjiiBEwZw5ccw28+64WRxQp\ndgoRiaWyMvQyRowIcx/XXRfuMLjddklXJiJNQdeJSIN8/TWMGQN33hluQ3v99XDGGZosFykUuk5E\nEvHJJ2GyfNQoOOkkePJJOOqopKsSkaTo/41SL++8A5ddFq4oX7UK3nwTxo1TgIg0d+qJSI3cw7Ik\nI0bAvHlw7bVhsnyXXZKuTETyhUJENrNuXbiXx4gRsGZNmCyfOFGT5SKyOU2syze++gpGjw6T5fvu\nGybLTz9dk+UixUgT65IzS5aEFXTHjIEePcL6Vl27Jl2ViBQC/R+zGZs/Hy65JNy7fO3acJHg2LEK\nEBGpP4VIM+MO06fDqafCaadB587w3nthCGvvvZOuTkQKjYazmol168I1HSNGhO3rr4eyMth226Qr\nE5FCphApYitXhl7HpEnh43vfg1tuCT0Q3QBKRHJBZ2cVmY8+CoHx7LPwyivhYsAzzwwf++6bdHUi\nki9ydXaWQqTAbdgAf/1rCI1nnw3LkvTsGULj1FNhp52SrlBE8pFCJNIcQ6RqmOrZZ2HyZNhttxAa\nZ50V7tmx1VZJVygi+U4hEmkuIfLhh5sOUx19dAiNM8+ETp2Srk5ECo1CJFKsIbJhQ1jksGqYaunS\nMEx11llwyikaphKReBQikWIKkZUrYdq0EBpTpoRhqrPOCh9HHaVhKhHJHYVIpNBD5MMPQ2hMmgSv\nvrrpMNU++yRdnYgUK4VIpNBCZP36MExVNb+xbNmmw1StWiVdoYg0BwqRSCGEyNdfbzybasoUaNt2\n4zDVkUdqmEpEmp5CJJKvIVI1TPXss/Daa3DMMRuHqbRGlYgkLS9CxMzaAE8CHYEPgD7u/mWWdh8A\nXwIbgHXufmRdx5vZIOByoBLo7+7TaqghkRBZvTosof7xx/DPf276sXgxLF8OZ5wRgqNHDw1TiUh+\nyZcQGQYsd/fhZvZroI27D8zS7n3gcHf/oj7Hm9lBwOPAEUAJMAPYP1taNEaIrFq1MRCqh0TV16tW\nQUnJph977RU+d+wIy5en6d49ldO64kqn06RSqaTL2IRqqp98rAnysy7VVD/5clOq3kC3aLsUSAOb\nhQhgZF92vqbjewFj3b0S+MDMFgNHAm/ErJevv84eCplfr1mzaSiUlMAhh4SeRdXXu+5a+yKGQ4cq\nROpDNdVPPtYE+VmXampacUOkrbtXALj7MjNrW0M7B6ab2XrgQXcfVcfxewKvZRy/JNpXqy+/rDsg\n1q3bGA5Vn7t0gV69Nn69yy5a5VZEpD7qDBEzmw60y9xFCIXBWZrXNK50rLsvNbPdCGGy0N1f3oLj\na3XQQSEkNmzYtPew117h7Kdzz9349c47KyBERHIl7pzIQiDl7hVmtjswy90713HMEOBrd7+jpuPN\nbCDg7j4sOuZ5YIi7bzacZWb5d2qWiEgByIc5kTLgUmAYcAnwTPUGZrY90MLdV5rZDsApwE11HF8G\nPG5mIwnDWPsBb2YrIBdvgoiINEzcnsguwDhgL+BDwim6K8ysPTDK3c80s32AiYShqq2Bx9391tqO\njx4bBPwUWEctp/iKiEhyCv5iQxERSU62027zhpmdZmb/a2Z/j64jydbmbjNbbGbzzOz7W3JsE9V0\nWMb+MWZWYWbzc1VPA2v6frSvxMxmmtnfzOxtM/tFHtS0rZm9YWZzo7puzlVNcerKeKyFmZWbWVmC\nNWX+TH1gZm9F71fWId8mqinzd6+1mf3FzBZG/4ZHJVmTmR0QvT/l0ecvc/WzHvN9GhS9P/PN7HEz\n2yYXNeWgrv7R34P6/U1w97z8IATcu4Sr2b8FzAO+U63N6cDkaPso4PX6HtvUNUVfHwd8H5ifJ+/T\n7sD3o+0dgUV58j5tH33eCnidcHZfou9VxuO/Ah4DyvKhJuB9wkW6efG7F339MHBZtL01sFPSNVV7\nnk+AvZKsKTrmfWCb6OsngYuT/vcDvgvMB7aNfv+mAZ1qe7187okcCSx29w/dfR0wlnBxYqbewCMA\nHs7cam2mBX2KAAADPUlEQVRm7ep5bFPXhIfTmr8gtxpck7svc/d50f6VwELqcT1OY9YUfb06arMt\n4RciV+9ZrLrMrAToCYzOUT2xa6LmC3kTqcnMdgKOd/eHoscq3f2rJGuq1uZk4D13/zjhmr4C1gI7\nmNnWwPaEcMuFOHV1Bt5w9/9z9/XAS8C5tb1YPofInkDmP/Q/2fwPXE1t6nNsU9VUrwslk67JzPYm\n9JJirwoQt6ZoyGgusAxIu/uCHNQUuy5gJDCABl7P1Eg1VV3I+1czuyIPatoH+MzMHoqGjx40s5YJ\n15TpAuCJHNQTqyYPS0DdDnwU7Vvh7jOSrgt4BzjezNpYOLO2J+HEpxrlc4g0hE73bQAz2xEYTzgL\nbmXS9bj7Bnc/jLBu2glm1q2uYxqbmZ0BVEQ9NyN/ftaOdfcuhF/2q83suITr2RroAtwb1bWa7Esh\nNTkz+xZhSaW/5EEtnQhDox2BPYAdzexHyVYF7v6/hEsupgNTgLnA+tqOyecQWQJ0yPi6JNpXvc1e\nWdrU59imrqmxxKop6kqPBx51982u80mipirRMMhkoGse1HUs0MvCYqJPACea2SMJ14S7L40+f0o4\nlf7IhGv6J/Cxu/9PtH88IVSSrKnK6cCc6L3KhTg1dQVecffPo2GjCcAP8qAu3P0hd+/q7ilgBfD3\nWl8tFxM5jfFBmNSpmhzahjA51Llam55snBw6mo2TQ3Ue29Q1ZTy+N/B2PrxP0dePAHfk0b/drkDr\naLslYUy2e9J1VWvTjdxNrMd5r7YHdoy2dwBeAU5J+n0CXgQOiLaHAMOSrina9wRwSZ78nB8KvA1s\nR+jVPgxcnXRd0de7RZ87AAuo48SInLyZjfUBnEY4Y2gxMDDa1w/4WUabP0Zv2FtAl9qOzYOa/kyY\nPPs/wljoZQnVdFi071hCV3UeodtaDpyW5PsEHBzVMTfaf32+/ExlPJ6zEIn5Xu2T8W/3dh79nB8K\n/DWqbQLRfwoSrml74FOgVb78PBHm1/5GOBuqFPhWntT1EmFuZC5hWapaX0sXG4qISIPl85yIiIjk\nOYWIiIg0mEJEREQaTCEiIiINphAREZEGU4iIiEiDKURERKTBFCIiItJg/w/5hdnjarXo7AAAAABJ\nRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa4bc37f650>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "energies = []\n",
    "betas = []\n",
    "for i in range(10):\n",
    "    beta_val = i*.01\n",
    "    e = compute_energy(R_T, beta_val)\n",
    "    betas.append(beta_val)\n",
    "    energies.append(e)\n",
    "\n",
    "plt.plot(betas, energies)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}