Variational/Variational_Helium.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sympy import *\n",
    "from sympy.abc import r,x,y,z\n",
    "from sympy import init_printing\n",
    "from cubature import cubature\n",
    "import sympy.printing.llvmjitcode as jit\n",
    "from timeit import default_timer as timer\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "init_printing()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Energy of a Helium atom\n",
    "\n",
    "The Hamiltonian for this system is\n",
    "$$-\\frac{1}{2} (\\nabla_1^2 + \\nabla_2^2) - \\frac{2}{r_1} - \\frac{2}{r_2} + \\frac{1}{r_{12}}$$\n",
    "\n",
    "The first two terms are the kinetic energy, the next two are the electron-nuclear Coulomb interaction, and the final term is the electron-electron Coulomb interaction.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "r1 = Symbol('r1')    # magnitude of vector r1 (electron 1 to nucleus distance)\n",
    "r2 = Symbol('r2')    # magnitude of vector r2 (electron 2 to nucleus distance)\n",
    "r12 = Symbol('r12') # |r2-r1| (magnitude of vector r2-r1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simple trial wavefunction\n",
    "\n",
    "This trial wavefunction uses a Slater form for the electon-nuclear interactions.  The value in the exponent (-2) is fixed by the cusp condition.   The electron-electron interaction uses a simple Pade form.   There is one adjustable parameter, $\\beta$.   The integration occurs over the spatial coordinates of each electron, giving a 6 dimensional integral.\n",
    "\n",
    "In the QMCPACK distribution, the example input file in `examples/molecules/He/simple_he.xml` has the same wavefunction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAIcAAAAYCAYAAADQ1+6cAAAABHNCSVQICAgIfAhkiAAAA6xJREFU\naIHt2ltoHFUYwPGfSaMWtfHWmCqCl4oKKsXSaq1UlEqpRW1BW4qCq30QrFEL0oKIpCiK8YKgRXwQ\nH3ywIIqC4oO3aJ8UxYpPKohQ8FLFu4JaGx++2ewlO5ud2Wyyq/OHYc+cmT3nn92z5zvnm1BQkMIh\ncy1QkJlLsQ6/43jch6/n1KigaxjEM1hYVXcVTkjKK3EZbm23o752G5iGC7ANo3gdqzrcX156xRN+\nwbH4LjkfQkklCizFWzh71s0ycCQeqDrfiD9w0tzopNIrnmXOwl11daMYrjo/AvfMllAezsNBnJ6c\nL8CE+PC7iV7xbMao2sFxnZhd2qKTYeUTrMAXyfnJyevnWITbsQfPYyd+xPwO+qTRK55pDOFMsVCF\nrbjI1Nmlq3kWjyTlNejHXjHK+3Bqcq0Pj866XYVWPPswgu3J8Z8kz1Z2RGUKbsR7eK6u7iYRK3eI\nKRuOw2di1X0wqRvEFlyP83O4zZbnlWLA7MMLuB8ftun7v+QK8aHD4TglKW/EiynvGe+sUkOyeN6h\nMmOMYX3OPie6+ZiX849qlVUibr8qFkwX4lt8icvxZof7b5Wsnk/i0KR8Lh7L2W+/WCOU1zBjOdtp\nh740h0aDYwluEQuvAZGF24pfM3Z6Gl7BUXX1g8nrYrVbyCzMlCP5PP9MjpV4B1/l9FyHl1TC01Kd\nC08LVXIj1bTsUMLHavf423H1TFq2wHiTayXd4bgAdze5XjK9Z7PwdKKYwS7BbW26ln0a0VKIXIG/\nxa+hzHIx1Q7NgFwWxlPqu8nxZjEbDGB13bVWPQ9TmbFeEwOCmNHLWdrFuHcGfEsp9WkO+qtuehpH\ni+lyDdYmkjvxwwzItcJ88Uxgg9gZ7MWBLnOETXhcpNx3JOX9OTz/wV9iEE2ILweW4QMxY2zB27gx\nue9afITfxDOVn0T4OkOEiPdTnJeIz7OeNIdJBsSX8FRKw91ALziS3bNReFpeVR7BDaY+bBsSu6hh\nlbCzq66dc3BncuyuKm9TmwBtGiKHxajp5oROLziS3bNReFpbdX232Em93OC9oypp8+mep5QyOkyO\nnu/F075Gu5dFuKZJw7NFLziSzXMTHhL/j7Ef34jE3YBIh28W2+ZBkbRrxno8kcO3kcMUxvCu2qzp\nMpFOPiZHp52gFxxpz/PiFvsYEhnezWJ7vAsPN7m/1GK7k1TLz8ODYuu1T6xiPxWx88DUt84JveBI\ne56r8UYHnAbxcwfaLSgoKCgoKChoi38BgAMLZLfLj4EAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$e^{- 2 r_{1}} e^{- 2 r_{2}} e^{\\frac{r_{12}}{2 \\beta r_{12} + 2}}$$"
      ],
      "text/plain": [
       "                    r₁₂     \n",
       "               ─────────────\n",
       " -2⋅r₁  -2⋅r₂  2⋅(β⋅r₁₂ + 1)\n",
       "ℯ     ⋅ℯ     ⋅ℯ             "
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "beta = Symbol('beta')\n",
    "R1 = exp(-2*r1)\n",
    "R2 = exp(-2*r2)\n",
    "G = exp(r12/2/(1 + beta*r12))\n",
    "R1*R2*G"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "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)\n",
    "\n",
    "def grad_spherical(e, r, rhat):\n",
    "    \"\"\"Compute gradient for expresion e with respect to symbol r.\n",
    "       The symbol rhat is the unit vector in the radial direction.\"\"\"\n",
    "    t1 = diff(e, r)\n",
    "    t2 = simplify(t1)\n",
    "    return t2*rhat\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Wavefunction has the form F1 * F2 * G\n",
    "# grad (FG) = G grad(F) + F grad(G)\n",
    "\n",
    "def grad_combo(F, gradF, G, gradG):\n",
    "    t = F*gradG + G*gradF\n",
    "    return simplify(t)\n",
    "\n",
    "# del^2 (FG) = G del^2(F0 + F del^2(G) + 2 grad(F) * grad(G)\n",
    "def del_combo(F, gradF, delsqF, G, gradG, delsqG):\n",
    "    t = G * delsqF + F * delsqG + 2*gradF*gradG\n",
    "    t2 = t.subs(rhat1*rhat1, 1)\n",
    "    t3 = t2.subs(rhat2*rhat2, 1)\n",
    "    t4 = t3.subs(rhat12*rhat12, 1)\n",
    "    t5 = t4.subs(rhat1*rhat12, rhat1_12)\n",
    "    t6 = t5.subs(rhat2*rhat12, rhat2_12)\n",
    "    return simplify(t6)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Unit vectors\n",
    "rhat1 = Symbol('rhat1')\n",
    "rhat2 = Symbol('rhat2')\n",
    "rhat12 = Symbol('rhat12')\n",
    "rhat1_12 = Symbol('rhat1_12')\n",
    "rhat2_12 = Symbol('rhat2_12')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAETEAAAAqCAYAAAB19ahPAAAABHNCSVQICAgIfAhkiAAAIABJREFU\neJzt3XncbVVd+PHPvczIBRSZBO0BUUJRMWfGqyHwyyF/pBDiQOqv7OeE4IBkdc0cICNMTc3SUxqO\nmZY2W9c0p1Qcf5ZjpmKoiaJSoMLvj3XO6znPec7ee+157X0+79fred3n7LP32uusvfb+ruHc9YAk\nSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkqdBW4LLp77sAjwbOBh7X0Tm7PG+MBwMHN5DO4mec2Q+4\npIH0u5L1OVJjeUuSJEmSJEmSJEmSJEmSJEmSJEmSJEmSpCRt7TsDHdgPOB/YPn29Hfg48EZgX+CA\nDs7Z5HkPrJMx4CDgPGBLxeN3BU5j+WecOZtyi6TM0mxTVrkt+xzHAYdWSKsJZfI5U7a8of0yr1re\nRWUvSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIStAqLmHwXuAy4dvr6x8DFwJ7Tn2s6OGeT531Avazx\nDeATNY4/G3gvyz8jwBrw9YppzrsVcG/gFODJpXO5WVa5LfscHwPOqJBWE8rkE6qVN2wu81TKu6js\nJUmSJEmSJEmSJEmSJEmSJEmSJEmSJElSglZhEZNFO4EbgP8HXAncOPLzNu1Q4Lqc948DPlkzzV2B\no4APAF8DDiyZXhOuAw7q4bxlVSlv2FjmqZX3UMq+rDOBI/rOhNSC3YFH9Z0JSappf+B1hMUGpdQZ\neyWNwdhjr2MA42DMlTQGY4+50ljY7pA0BrY7JDVtX+BN03+lVWHfQNIYjLlv4Bxgt4yLksZgzHFR\nkqRVZV9F0hjE9FVi5upy91nFRUwOA94FvBi4HNgG3Bd4EfAc4GWEBR7Gct5FBwFHT89d1s2AH+e8\nfzfgowvbij7jsjTvRlhQ46nAJcCnCtJow78Bdyp5TNfXs0p5w+YyT628q5R96s4F7gV8qe+MSA3b\nCryGzc8iSRqa7wBvBl7dd0akAsZeSWMx5tjrGMA4GHMljcWYY640FrY7JI2F7Q5JTbsWeCXwdmC3\nnvMidcG+gaSxGGvfwDnAbhkXJY3FWOOiJEmryr6KpLGI6avEzNXl7rNrvTwmYRvwc0u2fx34myXb\nHws8H/gRoZC3A/9EWP3qEuCbc/s+GPggcHWN85U57wnAMcAdgJfOHXsscMb097sDt5z+/mPCoig3\n5px30TeAc0rsP29/QoXKcgxh8YlbAkcB9yEE5MWynS/XZWneBNwA/A6hvPYEDiD7+hwF7EGz5fb9\n6ecpk9ayz9pmPmPKe1mdWizz1Mp7seyH7h7AkwnXQkrdvYDjgf0Idfa5hFiV5WLgSsLiR025BHgH\n8J7E0krVScDhwN6EtsWrgX/sMZ0hKVvf205H/Xs7oZ3/aOCPes6LVkcKsReai5nG3v7SGoImY6bx\ndxzGGHsdAxgPY27/jLnVGXO1aIwxV0pZCn19x9nLsd1Rne0OLbLdoRnvaTXlXcCpwAuBC3vOi1RG\nleegfYP++R2cauwXaJmx9Q2cA+yecbF/jplVZ2zUorHFRVXnPS1J6XF+f3jsq1Tn/6PTopi+Ssxc\nnfN5wM7pv78O7DX9/e7AicAWQmHPOwh4K3BIA+eMPe+Tp/++LCfN8zK239TRz1NyPuPMGjCZ/r74\nGRfLde8laf6vud/fABxK/vWpU24zOxde3w04rWRaRfWojXxCfnkvO+dimadW3nllPzS7Ap8A7td3\nRqQI+wAvmHt9FnAdcFjG/kcCn6fZv3h0EmGBq+2JpZWybxEaigAPA35AuJZ9pTMUZet72+koHXcB\nriK0l6S2pRB7obmYaeztN63UNRkzjb/jMqbY6xjAeBhz02DMrcaYqyxjirlSylLo6zvOXp7tjmps\ndyiL7Q55T6tptwS+R/hDT9IQVHkO2jdIg9/BKc9+gfKMpW/gHGD3jItpcMysGmOjsowlLqo672lJ\nSo/z+8NkX6Ua/x+dssT0VWLm6pbus7Vu7gZgL+B84BjgAuAPgCcBjwJuC7wXOBr44MJx3yAMOjZx\nzr2Al0Sc93eBmxGCRVlbCItP3H16ni0N/PwsYdGKxe1ZnxFgd+CJhFWnT17yGRfL9bq5NCGUzW7A\nfYFzgJcTVmTKuz51yi3rc9wR+JeSaRXVozbyWVTey845X+YplneVsk/VucAPgX/oOyNShCOBZxKe\nCwB/TbhHs/56wMWEFQZ/2ND59yU0ej6TWFqp2w68Ze717myMq12nMxRl63vb6SgdHwe+ADy874xo\nJfQde6G5mGns7T+t1DUZM42/4zKm2OsYQLOeCDyhp3Mbc9OwHWNuFcZcZRlTzO3SKYS/HDIWuwMv\novmFurSu776+4+zVbMd2RxW2O5TFdoe8p9W0bwFXABf1nREpUpXnoH2DNGzH7+CUZb9AecbSN1jV\nOcAxzdUZF6vZjmNmVRgblWUscVHVeU9L7XN+X2U5vz9M27GvUoX/j05ZYvoqMXN1g5/P2wpcSFig\n4xHT39u2g7CIR5fOBW6R8/55Oe+dNP33WRXOe5+F1wcBb2Xz538EsEeF9OftWEh3Ps0TK6ZTp9yW\neUyDae2gvXzGWHbOWZmnWN55ZT80HwWe2uP596Z4MabbE4LDHwK/B/w2sK3lfLXpNoT6nedw4DnA\n84EXA28krPQ1FlnXfRfgcrJXRdsC3Iv1xvYdgZuAuy7Zd1/C6oJHRuQntrwvJDyXdlJ/1ccm02rD\n04iL1WXr6uuBX6udu+bSaUqV51RR2ZWp73maSqcLbdW7oYmJjU8iLDjYtcOBlwEvBV4F/AnD+Yty\nRfF3VetVX7EXuo+/qcdeaOc52GTMTCn+phx7m06rTcbeIOXY27S+xgCG3IbI8n8Jk32vAR7f8bmN\nufUZc+MZc5thzA2GEnOHNGZ9O+CSHs7bhLxyvitwaUQa9vWH2ddfpXF2sN1Rhu2OZtjuCIbS7mjC\nqsdDqBYTy97TMXP83n/jsXgti9pXMycCN5D/XaU2Dakv0yS/g9PNc9Dv4NTjd3DKabpvsIr9ArBt\nMrMqfYO8OcAxztPB+ObqVikugmNmZThm1gzjYrAqcREcMwPHzNS8qmNmbRraeNgqz+87huf8ftNp\ntSH1Mbym06or5b7KUPopYFtppqm+SsxcXd/zebW8irBCD4SVeZ7X8vkOIjx4zmn5PPOeQBjMfFHO\nPvtlbN8d+CnC4hG/Slg97jJC5Ylx6pJtO9i8iMkewOmRaS6zrFzn01yWj6J06pTbMnclf/GaMmm1\nmc8YWeeclXlq5V1U9kNye0IAvnNP578H8C/TPGS5A/AfrK80eQjwJcK9P0T7Aa8lf5XJEwgrH84v\nxvQKQsMuFXVWySy67rcH/jgyrdcSGrnLPAz4ckQaseX9IODY6e87qddhajKtRU2sYPoThI7rjoL9\nytTVewLPJkwc1hmoaiqdZaqWXZXnVJX7PK++l9FUOvNSrXd9afMZOXNP4EdU7zRVyeMBwNuAQ+e2\nHQH8G2FQNWVF8XfV61XXsRe6j79txl5I8znYZMxsK/6uSuxtOi1Is871aQixNwV9jQEMuQ2R5ZFs\nvHd+g27HZ4259Rhz4xlz06xzfRpKzG07nymNWb+Z5udsuhBTzi8hfKEgi339Yfb1hzLODmnGANsd\nG9nu2Mx2R5DyOHvTxhAPIY2YmHdPx8zxe/+loc1naUxd2hO4Hnh0A/koa2h9mab4HZxunoPgd3Dq\n8Ds45XTRN0i5XwDp1ru+DKFv0Le8OcAxztPB+ObqhhIXIc1nlGNmGzlmtplxMXDMbN0QrjWk31d0\nzGw4+h4zy7NKc/uwuvP7juE5v990WotSbTPaV9nI/0e3mW2loMm+SsxcXZ/zebX8JPDfhIyfCzyG\nsHqV1h1PWMgE4LcIA6APjzjuGMJCFk+f/nv83Hs7GM+CElLXngJcR1i1r0vHAO8EJsAHyA4wuxIm\nS86f23Zr4JuEvA/RSwiLOWU5kvXFsOa9AvjXVnJUzauAw0oeE3vdITzbH1mQ3mMIq3VuyXj/FcCf\nF6QRW963Ah4193on1TtMTaa1TJVrs+j3CddnR84+Vevq44GPAPtUzVzD6cyrUnZVnlNVyq6ovsdq\nKp1Fqde7rrX9jITQpr4eeEDZzE1VyePTWb744MWkvyp0Xvy1XgU76Cb2Qvfxt+3YC2k/B5uMmU3H\n31WIvU2nNZNynevDEGJvCvoaAxhyGyLLskmgLv+qiDG3HmNuHGNukHKd68NQYm6b+UxpzPp4wl8y\nGZIy9eH2wFtz3revP7y+/pDG2SHtGGC7w3ZHFtsd6Y+zN20M8RD6j4lF93TRHD94/6Wi7WfpDorb\nVx8g1MuuDLEv0yS/g9PNcxD8Dk4dfgcnXhd9g9T7BZB+vevaEPoGfcubAxzjPB2Ma65uSHER0n5G\nOWbmmFkW46JjZvOGcq0h/b6iY2bDkcKYWZZVmduH1Z7fdwzP+X3H8Oyr+P/ollv1tlJbfZWYubqu\n5/MacTbw4b4zkbgHA/cGHkKoYI8EDp6+dyvCAwHgFOAX2bwIzKkLrw8CXk+3qydLY/J64P0952FC\ndoB5DHADsH9nuWnXkcDfF+xzKZsnk7YAnwXe0EamKpoAazWPz2tY3BL4BNn/ue5nWI8Ze2bk5cPA\ncwvyEVvevwA8C7ho+nMVoaFYZYC2ybSWmVDv2pwJ/DzFjeDYsrs3cDXhLzhAWPTtJuChJfPVVDp5\nJpQvuyrPqbL3eUx9j9FUOstMaqbXdL3r24R2n5Eznyasplr1HGslj3k58GdLtl8AXFYxH10oir/W\nq6Cr2Avdx9+2Yy+k9RxsMma2HX8njDv2Np3WvEnNtIy9m49vO/amoK8xgKG2IVJmzK3OmBvPmBtM\naqZlzN18fBcxd0J7+UxpzPp1wH36zkQNE4rrw5XAzZdst6+/fvyQ+vpDGmeHtGKA7Y7NbHdsZrtj\n8/EpjrM3aSzxEPqNiUX3dMwcv/dfOia0+ywtal9BaJN8ukYe6pgwjL5MU/wOzsbj23oOzvgdnGr8\nDk45bfcNhtAvgPTqXd8mpN836FveHKDzdO1wzKw6x8ziOWYWTGqmZVzcfLxjZsO41pB2X9Exs2GZ\n0P+YWZYJqzG3D6s7v+8Y3sbjnd93DK8s+yr+PzqfkRvF9FVi5uo27LNrxIlT8Fngf+ZebwHOI77w\nVsFVbFzo5TML781WQX43cDhhVZx5X114/Q1cwESq4ydJa8WtRecAnwe+03dGGvIE4LU57x8FfBT4\n8cL2RxIWbVq2atpYfQv4GnBfNndYTwYOJay0dgjrDfJ/X9jv9sArc85Rprxfs7DP44ErCKs/ltVk\nWk3bh9BAfxxhcjNLmbL7EaFRd9X09ZHAD4GPl8xbU+k0rexzqux9HlvfizSVThvaqHer4jt0OyD+\nMcJqnVcQYto1wB7AIwj9nlTlxV/r1bouYi/0E39Tjr3Q/HOwyZiZYvwdSuxtOq0mGXur6zr2Nq2v\nMYChtiFSZsytxphbjjG3PmNudSnH3JTGrE9kfaJ9rD4EnAa8cWG7ff04qfX1V2WcHWx3lGW7oz7b\nHdWl3O4oYjyMlxUTY+7pojl+77/xiLmWee2rmW8T4usW0vq+YEp9mab4HZx4dZ6DM34Hpzy/g1Ne\nm32DVegXgG2TOobcN8ibA3Serh2OmVXjmFk5jpnVZ1ysbshx0TGzeI6ZKUZTY2Z9SG08bFXn9x3D\ni+f8fj9SHsNrOq2mDKWvkmo/BWwr1RHTV4mZq9uwz1AWMbkS+FPgKYQPsBehgqc0Idm3D+e8d2vC\nwiUAZwEHEG7G+YVMUl5sQRqiWwMf7DsTGbYQOmkfIjQaTiM8E9aA5xCeuTMPA04hfJ5zpq/vQHim\nvHiaRpEm0ijyM8DvFeTh8unvfwj8gPDZDyGUxZc7ymcq/hl4MBs7X0cC7wC2Ley738LrbdOfa3PS\njynvRYcDT57ucyGhTr5jLr0y1yYvrb48C3hBxH5lyu7DwKuBJwI3Tt9/IPC5hfSKyi4mna6VeU7N\nlCm7mPoeU3ax901fmqx3q/SMhND5OqjD800IneRzCINjTwfOAH6FsOrvTGpxOS/+xt6Tq1K32o69\n0G/8TTH2QvPxNzZmDjH+DiX2xqbVl6br3Ko8I6H72AvNlm9fYwATumtD5EmtjVKVMbc6Y248Y24z\njLnV9RFzY3Q9Zp3ndoQvsNywsH1s9eTjwHFsXsTEvn68FPv6RW0F2x22O2x3lOc4e3Vdtzsc++7P\nYkyMvaeL5vi9/4Ih530m9loua1/Nu4bw1/z2I53/IJFSX6bJc/gdnHKqPgfB7+BU5XdwymmzbxBb\n3/0Ozuo+I2HYfYO8OcAJ3c7TDWG+ry7HzKrrY8wsttxSi42OmTXDuFjdkOOiY2blOGbWriHnfaap\nMbOupTYetsrz+47hlTO0+X3H8Nodw4tNq0tD6auk3E8B20p1xPRVYubqNuwzlEVMAH6n7wwM2FeA\ni6a/v6nPjEgrZBvw3b4zkeEAQiA4iBA4nz3dvh34J8LqZ5+e7nN3QmPsc8CbgRcSGhofJQSUomDb\nRBpFDpt+li/k7LMX8N/T36+e5uujhMbDdsIqdW3nMyUfIzTK5n0R2Dfi2FkjM6/zVVTen1pyzFeB\nZ0x/5lWpQ1lp9eU44Pvk19GZsmX3urnfF9tKZcouL50+xD6n5pUpu6L6Hlt2sfdNH5qsd6v2jITQ\nmTq4w/NdD9yPULanE1ZFfiehfGdSi8tF8TfmnuyinZCKtmMv9Bt/U4u90F78LYqZQ42/Q4m9MWn1\npek6t0rPSOg+9jZdvn2NAXTVhsiTWhulDmNuNcbccoy59Rlz6+k65sbqcsy6yGHANxa2jbGefJtQ\nvvPs65eTYl8/r61gu8N2h+2O8hxnr6fLdodj3/1ajIkx93RROXv/BWOoS2Wu5bL21bxrpv/ejHQW\nMUmpL9PUOfwOTnlVnoMzfgenPL+DU16bfYOY+u53cFb7GQnD7RtA/hxgl/N0Q5jva4JjZtX0MWZW\nttxSio2OmdVnXKxnqHHRMbPyHDNrzxjqUpNjZl1LaTwMVnd+3zG88oY0v+8YXjdjeEVpdW0ofZVU\n+ylgW6mumL5KzFzdhn2GtIiJJA3JLsCP+s5Ehl2m/x4B/NHc9p2ECZdLCCvHbQfeD+xOWDHsCuA9\n09+vAd4Sca4m0iiyxuZO57y92dhRuHju9+sJK+B/voN8puS/CCuOVjGrPz/MeD+mvF9a4nzbGfa1\n2Qo8Ffg/Eftadutin1Mzlt1GTde7VXtGQijDruP4EwiLDz4QeAnwAMJg2WnAJ0kvLq+RHX9j78km\n85O6NmMv+BxcZPwtz9hbTxt1bjvDLpOyuo6922m2fPscA+iiDZGnifTbzmMsY255xtzyjLn1GHPr\n66O/G6PLMesiB7L5PyZ0cd6uXcPmv4Kyhn39Muzrd8t2R3m2O+pxnL2+Ltsd23Hsu09VYuIa2eXs\n/bduO8PNO5S7llBcl7ZM//1xnUw1LKW+TFPnWMPv4JTld3C6Y7+gGvsG9dg2qW+ofQMongPsap6u\nqXT6Sj+WY2bl9RUbtzPccjMu1mNcrG+ocXENx8zKcsysPdsZbt6h+TGzrqU0HgarO7+/hmN4ZQ1p\nfn87w742juFVY1+lHttK9cX0VWLm6jbs4yImktSO64D9I/e9CzBh/QFd5ErgFyrkaebb038/z/qK\nYTNfBe4P7EFYDeszhJXa9gDeON3nK4SVyeZtA14DXAD8x9z2JtK4F3A8ocN1AvBcwgpyMweT/xev\nTwDem/He94DbEFZZe3/L+Vz0R4Rrv+g2wD2BG5a891jgIzlpxvo2m7+gHuu66b9Zx8eU967ED8CW\nqUNNafLa/BLh/l52zCLLbl3sc+r66TbLbqOm692qPSMhrKr5/YJ9mszjU4GTCF9mgDDQ8JvAUwir\nu96FbuNyzPXKi7+x92QX7YR5Y4294HNw0SrFX2NvdanXuVV6RkJc7G2y79xE+c7rawygqzZEXl1q\non1RJo9F9brO/WDMTeP5l2rZGXOrS73OGXOX6zqfXY5ZF12vrUvy12R/el7XcwPzdmHzdbSvX459\n/WKpx4BUy852R3Up1znH2ZdrKo9N9iuHGA9heDExr5y9/9YNvS6VuZZQXJf2nv77vZx9xvz9m6Jr\n2VS/xe/glOd3cPKl3EZblbKzb7AuhXq3as9I6HYesOlxvLw5wK7m6cp8rrx0mpqrK8qrc3XFUn5G\nxZbdkMvNuLguhTpnXFwupXGFGcfMynPMbLOU6vaYxsxWeW4f6s3v5517mZTm9x3DK29I8/u2GYNV\na2/bV1mXQr1btWckxPVVYubqYvaRJNX0WcIKW32aADdlvPcNlge9d0+POXRu268DV+ec53HAjulx\naxn7VE1jH+AFc6/PIjT+D1vY9r6ctJ9F9qJdfwl8vaN8xpqQXY6xx2dd95k7snHFuDK2Eho9j814\nv0x5l1F0bbowody1OYTwFxYW3USoR4ssu43KPKcsu3Vt1rtVeUZC6JBeXuMcayX230JYifPYJe89\nnpDfY+a2tR2XY69XXvwte0/2XbcmDDv2gs/BecbfYIKxt6oJadW5VXhGQr3YW0cTcRX6GQPoqg0R\nW5eaaF8U5bEonbqfwZib1vMvVt9lN8GYW9WEtOqcMTf+PGs1j+9zzDrmev008FcZ6TYR77owobg+\nnAu8emGbff2Nx9vXb8+EtGJArL7LboLtjqompFPn+n4+wjDaHROq5bFuv3K2bSzxENKNiVnl7P23\n3BDzXvZaQnFdej7rX+Lt2oT0+zIx5/A7ONWO9zs47ZiQThutjCGWHdg3mJmQTr1blWck9DMP2NQ4\nXtYcYB/f9amTTlNzdUV5da6ungnpPKNiDbHcwLg4MyGdOmdcjD/HWoXjmohDjpltPt4xs/6vxxDz\n3saYWZYJ45/bh+rz+3nn7tqE8vP7juFtPn6M8/u2GVervQ32VWYmpFPvVuUZCXF9lZi5ug37zFYb\nu8kff/zxx59aP4u+ANx8yfZUfAA4cMn22Yps35zbdl/CKvFZ/oDsDnPdNI4Engncdvr6rwmrep0w\nt883yS/ro1i+atzBwKmsr4TWdj5TcnNCw7aKGwmrwmX9hYEy5V1G0bWBsHLdWwgrzM27F+EvH+wA\n/g44uWIeyjoV+EngbXM/75y+9/PT12fO7W/ZbVTmOWXZrWuz3q3KMxLgcOBzHZ3rQMIzddn5XkkY\nLNtzblvbcTn2euXF37L35CrUrTZjL/gcnGf8rc7YW03bdW4VnpHQbeyd10RchX7GALpqQ8TWpSba\nF0V5LEqn7mcw5pZjzK3OmFuNMbcZfcXcGF2MWcdcr68DB2Sk20S8S8UtgasWttnXL8e+/nhiQKwh\nlp3tjmocZ29GH+2Ouv1KMB5WUSUmZpWz999yQ8x72WsJxXVpDfi3RnPZjFT6MjHn8Ds47fA7ON2w\nX1CPfYNqbJs0I7W+QZkyzJoD7OO7PnXSaWquLu8cTZzfMbNyUoiNMeUG6ZWdcbEa42IzUouLEBeH\nHDMrzzGz9g0x722MmfUhpfGwqvP7eedO0eL8vmN45Q1xft8242q1t8G+SlW2lZoR01dZo3iuLmYf\nSVJNzwW+1nMeJixfYAXgHMLqXXvNbdsCfAd409y2PYH/AX454nw3sXxVrzppbCE0YLZMX99xus9d\n5/a5Pdkdib0JDbTdlrz3HOA/gYM6ymesCe2vjva/gX+scY63sPmvdEK58i4j5to8jvRXrmN6/E1s\nbqBadpvFPqcsu2Jr1K93q/SMvBmhQ3rPGudYK7H/FkJZn7jkvW3AVwmDDtBNXI69Xlnxt+w9mULd\nmjDc2As+B2OssXrxd4Kxdwx1blWekXVjb1VNxNWZPsYAumpDxNSlJtoXMXmMrdd17gdjbj1rGHNj\nGHODCenUOWNuufOs1Ty+zzHrmOu1Ffgym//iRlP96S5MKK4PrwUesrDNvv7G4+3rr0YMiJVC2U2w\n3TH0OpfC8xGG0e6YUD6PTfUrxxQPId2YmDfHv2iN1b7/hpz3RWssv5YzRXXpSuB3auahqgnp92Vi\nzlHnOQirc98tHu93cNaNsY1WxpDLzr5BMCGNerdKz8g+5gGbHMfLmgPsap5uXt/zfWXy6lzdsJ9R\nsWLrjGNmm1nnjIupjZnN5D3bHTPbfLxjZo6ZQRpjZlkmjH9uH6rP7+edu2sTys/vO4a3+fixze/b\nZlzN9rZ9lWBCGvVulZ6RsX2VmLm6PufzVMNW4EnAM6Y/TdsPuKSDY7q0DTgPOJewsvUuLR3Thyr1\noe06pI1+mv47Mm+Y5mHvJe9tBd4PPG1u21mEILw2t+2+0zSOiThf1udtIo2Z1wK/vbBtC2H1zFst\n2f/+hFXLfovwmWceSJhcukeH+Yw1KUi7SN51n3k+8MIa5zgf+PSS7WXKu4w61+bOhFUvZyvX7Tvd\n56wK+ZhQ/56+3fT8z1vYbtltFvucsuyKNVHvVukZuR34NpsHM2NNKJ/HswnP1fnj9gXeDDx0blvX\ncRmyr1dW/C17T6ZQtyYF6RbpM/aCz8EYqxh/Jxh7x1DnVuUZuZ16sbeqJsu3rzGALtsQM8vqUhPt\ni7JpZKUD9e4HY249xtw4xtxgQjp1zpgbb0J7+exyzHom63q9Cbjbwram+tNdiKkPnyLMic2zr7/O\nvv7qxIBYKZTdBNsdQ69zKTwfYRjtjgnl89hUv3JM8RDSjYl5c/yLVv3+G3LeF2Vdy5m8urSN8OXJ\n02vmoaqh9GWa6Lf4HZyN/A7O+NtoZQy57OwbBBPSqHer9IzcTvfzgE2O4+XNAXY9T5fafF/eOZyr\nG/YzKlbdOuOYmXVuxrgYb0J/4wqOmW3kmNlmjplVV2fMLM+E1Zjbh/rz+318729elfl9x/A2GuP8\nvm3G1Wxv21cJJqRR71bpGbmd4r5KzFzdpn26/s8Bqu4BwNuArwB/SmhcfWTJfrcDbgF8sGT6ZwMH\nd3BMl7YDxxIe2hcAdwI+1sIxfYitD3WP0bp7AgcC74zc/93AfwHHA/8OMM+XAAAPtklEQVTeUp6W\nOYgQ1A4jrNAFYVXJTwF/APzJdNuNhDpxGWHS5IeETsy92ZjfWxOeJ5+pkacm0gB4DKGT9cyF7TcB\nfwucBLxx4b17A5cS7u1XEVZwuxnwA8I98J8d5rNNsdd95kTq5e8dhIbTvsC1c9vLlHcZda7NJ4H7\nAF+cSwvgcxXzUtW+wNuBo6evLwTuR1gM7G1YdsvEPqcsu2xN1rtVekaeRnjO/airDBLi19eBy4Hr\nCbFtF+BlwM65/bqOy3nXKyv+lr0nh1q3Uom94HMwj/G3HGNvfU3XuVV5RvYRe6G58oX+xgC6bENA\ndl1qon1RNo9V63XRccbcaoy55Rhz6zPmBkOJuSmOWUP+9Xo1cCYb51Ca6k+3pUx9uAfwHuB7C2nY\n17evn2esMSDWUMvOdkd9jrMHQ2h3NNWvXPV4CN3ExLw5/hnvv2DIeZ8pupYzeXXpdMJfGvyHlvK4\nzBD7Mk30W/wOTjp9A9to2ewXlGffoD7bJsHY+gZFZZg3B9j1PF1K831tnH/GMbNq+oiNdeuMY2bW\nuRnjYruaKl/HzBwzy+OYWTVNjJk1bYjjYVB/fr8Pdef3HcNLZwwP0mtvj7HNWMaQy86+Sn22lYI2\n+ioxc3V9zOd17oi+M9CS84FnTH+/FHjIkn1uQ7UVwdaABxFW5GnzmK5tJazcswdwJeFh0sYxfYip\nD00co42eSmjgx7oE+KuW8pKiJlafzErjZwgBG2DPJfscT2icLXpBzfxkqZrPGJOKx8W6Lc0sYPQe\n1j/rTFvlXUbKK9flseyqs+yq8xm52RbgS4S4UtWEfldjnqkbl2Ou17L42+Y92VbdmlQ4pow2Yy/4\nHKxjzGU3KUi3jjGXGwy37Ib6jGwi9nYhJq6OZQygzbZcE2nEpFP3MxhzmzfmspsUpFvHmMsNhlt2\nqx5zJ6TR541V53q9A9irwjmbindt+mPgJzLes68fx77+asaAWLY7qrHObeQ4e7t5rDL3PMR4CGnH\nxKw5/ljef5sNOe9FdekvgOe2eP7UtHktY9LxOzjx/A7OarXRyrDsqrPsNhryMzL1ecDYMuxyDrCp\nvwDfdlsi7xxNnN8xs+aNuewmBenWMeZyg2GW3arHxQn9jCvMOGYWzzGzciYVj4s15LzXqUsT0pyr\nzlL3OlWd3887dyqy5vcdw4s39vl924zVWXbVWXYbDfkZGdtXiZmr27TP1ur5StIa6Q52lnEYcNbC\ntpdPfwDuBHxoyXHPA15U4XzHEVY6avuYrt0I7A/8CvB0wmpIbRzTh5j60MQx2ugVhLoR++y8nPBM\nukNrOVoNJwOHAn8JHAKcMX09733AboSOxczewP90kcGpmHzG+F3gaw3ma9EFwMUNpHMp8IS5112X\ndxWzleueVvH4tq6NZVedZVedz8jlHkJY6fN9NdJoO49diL1ei/G3j3uyibo11NgLPgfrGHvZWW7W\nORjGM7KJ2JuKMY8BNFGXmmoPVk2nzHHG3GaNvewsN+scrFbMXaU+728Av9pS2n06HfhXwl/XWMa+\nfhz7+saAqmx3VGOdCxxnb0aVueehxkNIOyYum+OP5f0Xbyh5z6tLdyT89bvLWjz/EHQ5vuZ3cOL5\nHRzbaFVYdtVZdsFQnpEpzwOWKcOxzAH2PXbpXF2xsTyjqnDMrBrrXLBKcbHveTrHzOI5ZlaOY2bZ\n6tSlvp8ZTShznarM7w9B3vy+Y3jxVnl+3zZjdZZddZZdMJRnZExfJWaurnCfI4AnAb9ICNp3LZvT\nBFxCWPVl0RnAIwkLg/wmsGuXmSqwDXgLcJuF7ecRLtqiE4CLlmy/K/CanPNklcHdpudeI6zIE6PK\nMXVdApyU8V7M9X0XcGDJc1Y5poxHAc8GHgdcmLFP0WfLqg95qhxTR91r14es+xLCCmAPLJHWE4HX\nNZGphJ1LWCDnJuANhM/cVBpHAtdOt8//7LskjaOAK1iPA/cn1LEmNZHPPp1MmFRryl8D953+3kZ5\nV9HFXxFommVXnWVXnc/IzXYBPgAc3XdGaqobl8ter/n429Y9OeS61WbsBZ+DdVh21Vhu1fmM3GwI\nsbdsXB3yGECbdamp+liUTpOfwZjbHMuuGsutOmPuZkOIuU1r6nr9PLA98pxDqAu7Ac9i+VzqPPv6\n+ezrr14MKGtoZWe5Vec4+2aptzuqzj0bD5drIiYuzvHH8v7baMh5h+K69JfAQzvKS9/avpZ+B6dZ\nfgcnWIU2WhWWXXWWXTD0Z2TKfYMqZdj2HGAT38HNS6fJepOXV+fqglV4RlUxtLKz3KozLm6WclyE\ncnHIMbNijpmlY8h5h+bHHVLV5HUqM7+fd+5UxMzvO4ZXbBXm920zVmfZVWfZBUN/Rsb2VWLm6gr3\n+a2537cQAtiQ3Bl48pLtZwC3nnv954RFONr2sIh9HgfsYPkNuwfwsoVt+xIWvVjmNwkVfpm8MngE\nYcGUpwHvBe5TmOtqxyyKKZ+Zk4CrWd6QzPts8420K4AzI85V5ZhFMZ/tPOCF09/XgBuAWy7sU1R3\n8+pDlirHzCtz3aD6tWtL3fsS4GcJgaWM1xHuG7XvRODx09+fTViQRsEuhDLZpcE0DwH+mfVnSwrl\nvezePRl4LCG/hxBWiKsSu9pi2VVn2VWXStml5AWUb+somMVf69VGbcdeSKfMfQ5WN7Sys9yqS6Xs\nUjLW2OsYwDgYc5tj2VVjuVWXStmlZKwxV+2yr7+cfX1jQIyhlZ3lVl0qZZeSsbU7jIfZmoyJ83P8\nsbwm41FUly4Afq277GiB38HJ5ndwbKMVseyqs+zGYWx9A3AOsGuOmTXHsqvGcqsulbJLydjiomNm\n2RwzU1PaGHfQeDmGl21V5vdtM1Zn2VVn2Y1DTF8lZq4uaj7vs8Cxc69fVXRAYp4N3Glh2xFsXBTh\nDsAXgb07yM+OEvtmLZbwh2xcWeeXCKuo7QacurDvB4BjlqQRWwZrwCQyv3WOmdkRud++hNWHdrJ5\nIYyiz3Yp8MvT398L/FTE+aocs2hHwfu7A98EfmJu25EL+8Rct7z6kKXKMfN2lNi3zrVry44S+2bd\nl0cCnyh53l2BV7L5OktjcEdCg6VvQ1+5rk+WXXWW3Xjsz7gmiTRuqcRe8DlYh2VXjeU2HmOOvY4B\njIcxdxwsu2ost/EYc8yVxsR2xzhYdtVYbuNhu0NS07ax/h0uaVWk0jewjVadZVedZTceY+0bOAfY\nvVTiIviMqsOyq8ZyG4+xxkVJklZZKn0V24zVWXbVWXbjEdNXiZmri57PewpwA2ERh4uBu0y3Hzp9\n7z3Am4HnANcAe2WkU3b/ptJ4O5tXxXro9N+LgcuAbxE+VxN5LLKjxL5ZiyU8CThj+vvZhJv4W4S8\nHruw71XALZakkVUG83YHXgR8mrDaERSX0bJjytgRud+FwB4sXwij6LMdCzycsHDHU+e25322rGPK\n2FHw/imEa34y4aF9GXDawj5Fny2rPuR9tqI6FGNHiX2rXLu2780dJfbNui/3JyxCI0mSJEmSJEmS\nJEmSJEmSJEmSJEmSJEmSljgYuAL4IvB14H7T7acTFgf5GGHBha3AEdP3thIWIJhXZv+thIU6njH9\nKUoja3+Av1vymWYLJZwDPBP4yPTfvDyWyWeeHZH7QfZiCQ+dnjfGDYRVnpelAZvLoEhRGdW1I2Kf\nB7G+0MZOshfCGNpnezjhmp80fb0P8B3g8Ll9hvrZZqpeuyaeH3l2RO4H2fflbsD1JdKRJEmSJEmS\nJEmSJEmSJEmSJEmSJEmSJGn0ZotebANeQVg04EbgWcBbgeOAvwEOAG4NvH76/peA/YDHsnlxgjL7\nPwB4G/AV4E+BuxEWM8hK40EZ+0NY+GDewcB/TH9//fTfHwCnAJdkpE/JfM6f63xgy9y2E4E9515/\nD3ge5XwXODpy32uBvaf/zucrqwyKZF0DCItHvAi4YO71E4C9pq8vXUirSvncCrg58BcZ+RvyZ/vu\n9N8PT//9PvDfwAMJ9+GQPxvUu3ZZ92ZK9+UBwH+VPEaSJEmSJEmSJEmSJEmSJEmSJEmSJEmSpJVw\nJuE/+897HnDW9PezCIuaLLNzybbY/c8HnjH9/VLgIQVp5O3/twv7/hxhkYZ5rwaeFpHHMvnMsiNi\nn5mbgLUl2x/M5uuS5ZNL0igqgyLLymg/wiIYH53b9iDCohOwvphEkR0F7/8CYTGdi6Y/VwG/T1i4\nAob92Q4jLM6xz9y2q4AnTn8f8meD+teuiedHlh0R+8xk3Zd3Aj5YIh1JkiRJkiRJkiRJkiRJkiRJ\nkiRJkiRJkkZvtpDA54HjFt7bhfX/qH9/4F0l0o3d/+XTHwgLA3yoII28/a8Gts29PhG4zdzrOwNH\nA79bMo9F523TLYD/jNz3fYS8zSsqgyLLyui7wGXAtXPbbgucM/39C6wvjFHHa4AXAC+c/twAXAG8\nc/r+kD/b1wiLcZwwfX0gcDPgz6avh/zZoP61a+L50aajgX/o6FySJEmSJEmSJEmSJEmSJEmSJEmS\nJEmSJA3CrtN/PwH8FfDbwFeBPYC/B748ff8owqIEsWL3v376cwLwbuCqgjTy9n83cA/WFxf4IHAm\n8D3CgiyHAKcTFlQok8ei89ZxLmFBBwiLPbwXeOnc+3cGXhyZ1juAU4C/mNtWVAZFYsvo5cDu09/v\nBFwemX6Mw4EnE/J+IbAP4bMO/bM9Gvh14A7AMcBDCIubwPA/20zVa9fE86OOovvyDOD3GzqXJEmS\nJEmSJEmSJEmSJEmSJEmSJEmSJEma2llz/32BZ5c4Pmv/WwDPm/5+JHCvkvlatDPyvFnOqXl+gD8u\nse9W4G+APaevmyiDPDuXbDsBuCjy+Drl42erbueSbV19Nqj/+XYuvO76vjwA+JOaaUiSJEmSJEmS\nJEmSJEmSJEmSJEmSJEmSpCV21tz/l4Ddpj+nRhyft/+FwIHAw6fv17GzxHnbcG/g3JLHnAQ8Zfp7\nE2WQZ+fC67KLSdThZ6tu58LrLj8b1P98Oxded31fPpewEIskSZIkSZIkSZIkSZIkSZIkSZIkSZIk\nSWrIXsD5wNXABdPXZfc/G7gW+BZwDXBsQRpF++8C/CLw2NgP0VI+69oNuAjYWuHYBwD3pF4ZxNi5\n8LrLxST8bNXtXHjd9SIgVT9fCvflPYBjWj6HJEmSJEmSJEmSJEmSJEmSJEmSJEmSJElSYw6heFGY\nvqSwmERb/GySJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmS\nJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmS\nJEmSJEmSJEmSJEmSJEmSJEmSJEmK8v8BwrdrSWSQ/EMAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\frac{e^{- 18 r_{1} - 18 r_{2} + \\frac{1}{\\beta r_{12} + 1} \\left(14 r_{1} \\left(\\beta r_{12} + 1\\right) + r_{12} + 14 r_{2} \\left(\\beta r_{12} + 1\\right)\\right)}}{8 r_{1} r_{12} r_{2} \\left(\\beta r_{12} + 1\\right)^{2} \\left(\\beta^{4} r_{12}^{4} + 4 \\beta^{3} r_{12}^{3} + 6 \\beta^{2} r_{12}^{2} + 4 \\beta r_{12} + 1\\right)} \\left(- 16 r_{1} r_{12} \\left(\\beta r_{12} + 1\\right)^{2} \\left(\\beta^{4} r_{12}^{4} + 4 \\beta^{3} r_{12}^{3} + 6 \\beta^{2} r_{12}^{2} + 4 \\beta r_{12} + 1\\right) + 8 r_{1} r_{2} \\left(\\beta r_{12} + 1\\right)^{2} \\left(\\beta^{4} r_{12}^{4} + 4 \\beta^{3} r_{12}^{3} + 6 \\beta^{2} r_{12}^{2} + 4 \\beta r_{12} + 1\\right) - r_{1} \\left(- 8 r_{12} r_{2} \\hat{r}_{2 12} \\left(\\beta^{4} r_{12}^{4} + 4 \\beta^{3} r_{12}^{3} + 6 \\beta^{2} r_{12}^{2} + 4 \\beta r_{12} + 1\\right) + 16 r_{12} \\left(r_{2} - 1\\right) \\left(\\beta r_{12} + 1\\right)^{2} \\left(\\beta^{4} r_{12}^{4} + 4 \\beta^{3} r_{12}^{3} + 6 \\beta^{2} r_{12}^{2} + 4 \\beta r_{12} + 1\\right) + r_{2} \\left(\\beta r_{12} + 1\\right)^{2} \\left(4 \\beta r_{12} + r_{12} + 4\\right)\\right) - 16 r_{12} r_{2} \\left(\\beta r_{12} + 1\\right)^{2} \\left(\\beta^{4} r_{12}^{4} + 4 \\beta^{3} r_{12}^{3} + 6 \\beta^{2} r_{12}^{2} + 4 \\beta r_{12} + 1\\right) - r_{2} \\left(8 r_{1} r_{12} \\hat{r}_{1 12} \\left(\\beta^{4} r_{12}^{4} + 4 \\beta^{3} r_{12}^{3} + 6 \\beta^{2} r_{12}^{2} + 4 \\beta r_{12} + 1\\right) + r_{1} \\left(\\beta r_{12} + 1\\right)^{2} \\left(4 \\beta r_{12} + r_{12} + 4\\right) + 16 r_{12} \\left(r_{1} - 1\\right) \\left(\\beta r_{12} + 1\\right)^{2} \\left(\\beta^{4} r_{12}^{4} + 4 \\beta^{3} r_{12}^{3} + 6 \\beta^{2} r_{12}^{2} + 4 \\beta r_{12} + 1\\right)\\right)\\right)$$"
      ],
      "text/plain": [
       "                                                                              \n",
       "                                                                              \n",
       "⎛                       2 ⎛ 4    4      3    3      2    2              ⎞     \n",
       "⎝- 16⋅r₁⋅r₁₂⋅(β⋅r₁₂ + 1) ⋅⎝β ⋅r₁₂  + 4⋅β ⋅r₁₂  + 6⋅β ⋅r₁₂  + 4⋅β⋅r₁₂ + 1⎠ + 8⋅\n",
       "──────────────────────────────────────────────────────────────────────────────\n",
       "                                                                              \n",
       "                                                                              \n",
       "\n",
       "                                                                              \n",
       "                                                                              \n",
       "                 2 ⎛ 4    4      3    3      2    2              ⎞      ⎛     \n",
       "r₁⋅r₂⋅(β⋅r₁₂ + 1) ⋅⎝β ⋅r₁₂  + 4⋅β ⋅r₁₂  + 6⋅β ⋅r₁₂  + 4⋅β⋅r₁₂ + 1⎠ - r₁⋅⎝- 8⋅r\n",
       "──────────────────────────────────────────────────────────────────────────────\n",
       "                                                                              \n",
       "                                                                              \n",
       "\n",
       "                                                                              \n",
       "                                                                              \n",
       "             ⎛ 4    4      3    3      2    2              ⎞                  \n",
       "₁₂⋅r₂⋅r̂₂ ₁₂⋅⎝β ⋅r₁₂  + 4⋅β ⋅r₁₂  + 6⋅β ⋅r₁₂  + 4⋅β⋅r₁₂ + 1⎠ + 16⋅r₁₂⋅(r₂ - 1)\n",
       "──────────────────────────────────────────────────────────────────────────────\n",
       "                                                                              \n",
       "                                                                              \n",
       "\n",
       "                                                                              \n",
       "                                                                              \n",
       "            2 ⎛ 4    4      3    3      2    2              ⎞                 \n",
       "⋅(β⋅r₁₂ + 1) ⋅⎝β ⋅r₁₂  + 4⋅β ⋅r₁₂  + 6⋅β ⋅r₁₂  + 4⋅β⋅r₁₂ + 1⎠ + r₂⋅(β⋅r₁₂ + 1)\n",
       "──────────────────────────────────────────────────────────────────────────────\n",
       "                                                                              \n",
       "                                                           8⋅r₁⋅r₁₂⋅r₂⋅(β⋅r₁₂ \n",
       "\n",
       "                                                                              \n",
       "                                                                              \n",
       "2                    ⎞                        2 ⎛ 4    4      3    3      2   \n",
       " ⋅(4⋅β⋅r₁₂ + r₁₂ + 4)⎠ - 16⋅r₁₂⋅r₂⋅(β⋅r₁₂ + 1) ⋅⎝β ⋅r₁₂  + 4⋅β ⋅r₁₂  + 6⋅β ⋅r₁\n",
       "──────────────────────────────────────────────────────────────────────────────\n",
       "    2 ⎛ 4    4      3    3      2    2              ⎞                         \n",
       "+ 1) ⋅⎝β ⋅r₁₂  + 4⋅β ⋅r₁₂  + 6⋅β ⋅r₁₂  + 4⋅β⋅r₁₂ + 1⎠                         \n",
       "\n",
       "                                                                              \n",
       "                                                                              \n",
       " 2              ⎞      ⎛                ⎛ 4    4      3    3      2    2      \n",
       "₂  + 4⋅β⋅r₁₂ + 1⎠ - r₂⋅⎝8⋅r₁⋅r₁₂⋅r̂₁ ₁₂⋅⎝β ⋅r₁₂  + 4⋅β ⋅r₁₂  + 6⋅β ⋅r₁₂  + 4⋅β\n",
       "──────────────────────────────────────────────────────────────────────────────\n",
       "                                                                              \n",
       "                                                                              \n",
       "\n",
       "                                                                              \n",
       "                                                                              \n",
       "        ⎞                 2                                                  2\n",
       "⋅r₁₂ + 1⎠ + r₁⋅(β⋅r₁₂ + 1) ⋅(4⋅β⋅r₁₂ + r₁₂ + 4) + 16⋅r₁₂⋅(r₁ - 1)⋅(β⋅r₁₂ + 1) \n",
       "──────────────────────────────────────────────────────────────────────────────\n",
       "                                                                              \n",
       "                                                                              \n",
       "\n",
       "                                                                     14⋅r₁⋅(β⋅\n",
       "                                                    -18⋅r₁ - 18⋅r₂ + ─────────\n",
       " ⎛ 4    4      3    3      2    2              ⎞⎞⎞                            \n",
       "⋅⎝β ⋅r₁₂  + 4⋅β ⋅r₁₂  + 6⋅β ⋅r₁₂  + 4⋅β⋅r₁₂ + 1⎠⎠⎠⋅ℯ                          \n",
       "──────────────────────────────────────────────────────────────────────────────\n",
       "                                                                              \n",
       "                                                                              \n",
       "\n",
       "r₁₂ + 1) + r₁₂ + 14⋅r₂⋅(β⋅r₁₂ + 1)\n",
       "──────────────────────────────────\n",
       "        β⋅r₁₂ + 1                 \n",
       "                                  \n",
       "──────────────────────────────────\n",
       "                                  \n",
       "                                  "
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# TODO: Explain this derivation\n",
    "grad_1 = grad_spherical(R1, r1, rhat1)\n",
    "del2_1 = del_spherical(R1, r1)\n",
    "\n",
    "grad_2 = grad_spherical(R2, r2, rhat2)\n",
    "del2_2 = del_spherical(R2, r2)\n",
    "\n",
    "gradG_1 = grad_spherical(G, r12, rhat12)\n",
    "del2G_1 = del_spherical(G, r12)\n",
    "\n",
    "grad_R1G = grad_combo(R1, grad_1, G, -gradG_1)\n",
    "\n",
    "delsq_R1G = del_combo(R1, grad_1, del2_1, G, -gradG_1, del2G_1)\n",
    "\n",
    "\n",
    "grad_R2G = grad_combo(R2, grad_2, G, gradG_1)\n",
    "delsq_R2G = del_combo(R2, grad_2, del2_2, G, gradG_1, del2G_1)\n",
    "\n",
    "# This is the integrand for the denominator\n",
    "psi2 = R1*R1*R2*R2*G*G\n",
    "\n",
    "h1 = -1/S(2) * (R2*R2) * R1 * G * delsq_R1G - 2*psi2/r1\n",
    "h2 = -1/S(2) * (R1*R1) * R2 * G * delsq_R2G - 2*psi2/r2\n",
    "h = h1 + h2 + psi2/r12\n",
    "h = simplify(h)\n",
    "# This is the integrand for the numerator\n",
    "h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Symbols for performing the integration\n",
    "x1,y1,z1 = symbols('x1 y1 z1')\n",
    "x2,y2,z2 = symbols('x2 y2 z2')\n",
    "# Add a very tiny softening to keep the evaluation of the integrand from blowing up at certain points.\n",
    "r1_c = sqrt(x1*x1 + y1*y1 + z1*z1 + 1e-20)\n",
    "r2_c = sqrt(x2*x2 + y2*y2 + z2*z2 + 1e-20)\n",
    "r12_c = sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) + (z2-z1)*(z2-z1) + 1e-20)\n",
    "rhat1_12_c = (x1*(x2-x1) + y1*(y2-y1) + z1*(z2-z1))/r1_c/r12_c\n",
    "rhat2_12_c = (x2*(x2-x1) + y2*(y2-y1) + z2*(z2-z1))/r2_c/r12_c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "For beta =  1.0  Energy =  -2.82611479297\n",
      "Time: 96.18 seconds\n"
     ]
    }
   ],
   "source": [
    "# Perform 6D integral. Using scipy.integrate.nquad will take too long.\n",
    "# Use Cubature and Sympy LLVM JIT \n",
    "#  TODO: insert link to explain different integrators\n",
    "\n",
    "beta_val = 1.0\n",
    "num = tmp.subs(beta, beta_val)\n",
    "denom = simplify(tmp2.subs(beta, beta_val))\n",
    "\n",
    "ndim = 6\n",
    "lim_min = [-30.0]*ndim\n",
    "lim_max = [30.0]*ndim\n",
    "\n",
    "repls = [(r1, r1_c), (r2, r2_c), (r12, r12_c), (rhat1_12, rhat1_12_c), (rhat2_12, rhat2_12_c)]\n",
    "repls_n = [a for a in repls]\n",
    "tmp_h = h.subs(beta, beta_val)\n",
    "tmp_d = psi2.subs(beta, beta_val)\n",
    "\n",
    "# Use CSE (Common Subexpression Elimination) to speed up the integrand evaluation\n",
    "cse_h = cse([tmp_h, tmp_d])\n",
    "repls.extend(cse_h[0])\n",
    "expr_h = [repls, cse_h[1]]\n",
    "\n",
    "jit_expr = jit.llvm_callable( [x1, y1, z1, x2, y2, z2], expr_h, callback_type='cubature')\n",
    "\n",
    "epsrel = 1e-4\n",
    "start = timer()\n",
    "res = cubature(jit_expr, ndim=ndim, fdim=2, xmin=lim_min, xmax=lim_max, relerr=epsrel)\n",
    "end = timer()\n",
    "E = res[0][0]/res[0][1]\n",
    "print(\"For beta = \",beta_val,\" Energy = \",E)\n",
    "print(\"Time: %.2f seconds\"%(end-start))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Collect computation of energy into a single function\n",
    "def compute_energy(beta_val):\n",
    "\n",
    "    repls = [(r1, r1_c), (r2, r2_c), (r12, r12_c), (rhat1_12, rhat1_12_c), (rhat2_12, rhat2_12_c)]\n",
    "    repls_n = [a for a in repls]\n",
    "    tmp_h = h.subs(beta, beta_val)\n",
    "    tmp_d = psi2.subs(beta, beta_val)\n",
    "\n",
    "    cse_h = cse([tmp_h,tmp_d])\n",
    "    repls.extend(cse_h[0])\n",
    "    expr_h = [repls, cse_h[1]]\n",
    "\n",
    "    jit_expr = jit.llvm_callable( [x1, y1, z1, x2, y2, z2], expr_h, callback_type='cubature')\n",
    "\n",
    "    epsrel = 1e-4\n",
    "    res = cubature(jit_expr, ndim=ndim, fdim=2, xmin=lim_min, xmax=lim_max, relerr=epsrel)\n",
    "    E = res[0][0]/res[0][1]\n",
    "\n",
    "    return E"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.01 -2.86001117585\n",
      "0.03 -2.86676650391\n",
      "0.05 -2.87135272931\n",
      "0.06999999999999999 -2.87442400954\n",
      "0.09 -2.87639789443\n",
      "0.11 -2.87755791558\n",
      "0.13 -2.8781022511\n",
      "0.15000000000000002 -2.87817447072\n",
      "0.17 -2.87788106408\n",
      "0.19 -2.87730296758\n"
     ]
    }
   ],
   "source": [
    "# Compute energies over a range of variational parameters to find the minimum.\n",
    "energies = []\n",
    "betas = []\n",
    "for i in range(10):\n",
    "    beta_val = i*.02 + 0.01\n",
    "    e = compute_energy(beta_val)\n",
    "    print(beta_val, e)\n",
    "    betas.append(beta_val)\n",
    "    energies.append(e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x7f38a6a972b0>"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAETCAYAAABwaNKCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xu0VWW9//H3FzYXwQQFRRAhlYuA3LaapKYbEMX7SFFR\nIwHFzNAUjyXDGuKvzslMy3IMMg3zWCAaXsILCCi7kxmSckc01JBzuKiYaGpsLvv7++OZOxbbvWXv\nvS7PunxeY6zBWvO2vnO64OMz5/PMae6OiIhIrjWLXYCIiJQmBZCIiEShABIRkSgUQCIiEoUCSERE\nolAAiYhIFNECyMxuN7M1ZrbMzB41s/3qWe56M1tlZivMbLqZtUyZd02yjZVmdlvK9MlmtjaZd2ou\n9kdERBonZgtoHtDP3QcBa4HJtRcwsy7ANUC5uw8AyoDRybyhwNlAf3fvD9yRTO8DXAj0AU4HppqZ\nZX93RESkMaIFkLsvcPfq5OMioGs9izYH2ppZGdAG2JhMvwq4zd13Jtvbkkw/F5jp7jvdfR0h3L6U\nhV0QEZE05Ms1oPHAnNoT3X0jcCewHtgAbHX3BcnsXsBJZrbIzBaa2dHJ9EOA/03ZzIZkmoiI5JGs\nBpCZzU+u3dS8ViZ/np2yzM3ADnefUcf67Qktmu5AF2BfM7skmV0G7O/uQ4DvAL/P5r6IiEhmlWVz\n4+4+4vPmm9lY4AxgWD2LnAK85e7/SJZ/DDgemAH8H/BY8j1/NbNdZtaB0OLplrKNrsm0ur5fN8IT\nEWkCd0/72nrMXnAjgRuBc9y9qp7F1gNDzKx10pFgOLAmmfcESXCZWS+gpbu/D8wGLjKzlmZ2GNAD\nWFxfHe6e1uvss51f/Sq9bRTL65ZbboleQzG9dDx1LPP1lSkxrwHdDewLzDezJWY2FcDMOpvZUwDu\nvhiYBSwFlgMG3Jusfz9wuJmtJLSIvp6s8yrwCPAq8AxwtWfyiNUyaRL87GdQXb33ZUVEZLesnoL7\nPO7es57pm4CzUj7fCtxax3I7gDH1bONHwI8yU+nnO/lkaN0a5s6FM87IxTeKiBSHfOkFV7DMdreC\nSl1FRUXsEoqKjmfm6FjmJ8vi2am8Z2YZOTu3fTscdhjMmQMDBmSgMBGRPGZmeCF3QigmLVvCt76l\nVpCISGOoBZSh/X//fejRA9asgYMPzsgmRUTyklpAeaZDB7j4Ypg6NXYlIiKFQS2gDO7/66/DV74C\nb78N++yTsc2KiOQVtYDyUO/ecNxx8Lvfxa5ERCT/qQWU4f1fuDB0SFi1Cpop3kWkCKkFlKcqKqBV\nK3j22diViIjkNwVQhpnB9derS7aIyN7oFFwW9n/7dvjiF0MrqH//jG9eRCQqnYLLYy1bwsSJagWJ\niHwetYCytP81A1Nfew06dcrKV4iIRKEWUJ7r0AFGj9bAVBGR+qgFlMX9f/11OOkkWLdOA1NFpHio\nBVQAeveGY4+F6dNjVyIikn/UAsry/j//fOiQsHp16KItIlLo1AIqEEOHhl5xGpgqIrInBVCW1QxM\n/elPY1ciIpJfdAouB/tfVRWemDpvHhx1VNa/TkQkq3QKroC0aqUnpoqI1KYWUI72f8sW6NlTA1NF\npPCpBVRgOnaEiy6CX/4ydiUiIvlBLaAc7v9rr8HJJ2tgqogUNrWACtCRR2pgqohIDQVQjtU8K6iE\nG54iIoACKOeGDYOystAlW0SklCmAcswMJk3SwFQREXVCiLD/VVXhianz52tgqogUHnVCKGA1A1Pv\nuit2JSIi8agFFGn/awamvv46HHRQlBJERJpELaAC17EjXHihBqaKSOlSCyji/q9ZEx7XsG4dtG4d\nrQwRkUZRC6gI9OkDRx+tgakiUpoUQJFNmqSBqSJSmhRAkQ0bBs2bhy7ZIiKlRAEUmZ6YKiKlSp0Q\n8mD/awamLlgA/frFrkZE5POpE0IRadUKrr5aA1NFpLSoBZQn+//ee9Crlwamikj+UwuoyBx4oAam\nikhpUQsoj/ZfA1NFpBCoBVSE+vSB8nKYMSN2JSIi2acAyjM1zwrKo4aZiEhWKIDyzPDh0KxZ6JIt\nIlLMFEB5RgNTRaRUqBNCHu5/zcDU556Dvn1jVyMisqeC74RgZreb2RozW2Zmj5rZfvUsd72ZrTKz\nFWY23cxapsy7JtnGSjO7LZnW3cw+NbMlyWtqrvYpUzQwVURKQbQWkJmdAjzv7tVJeLi7T661TBfg\nBeBId99uZg8DT7v7g2Y2FJgMnOHuO82so7tvMbPuwJPuPqABNeRlCwh2D0z929/CGCERkXxR8C0g\nd1/g7tXJx0VA13oWbQ60NbMyoA2wMZl+FXCbu+9MtrclZZ20D0xsBx4IF1yggakiUrzypRPCeGBO\n7YnuvhG4E1gPbAC2untN/7BewElmtsjMFprZMSmrfjE5/bbQzE7MdvHZct11IYC2bYtdiYhI5pVl\nc+NmNh/olDoJcOBmd38yWeZmYIe7f2b4pZm1B84FugMfArPM7JJk2TJgf3cfYmbHAo8AhwObgG7u\n/oGZlQNPmFlfd/+4rhqnTJny7/cVFRVUVFSkudeZ07cvDBoEDz0E48bFrkZESlVlZSWVlZUZ327U\nXnBmNhaYAAxz96o65o8CTnP3CcnnMcBx7j7RzOYQTsH9MZn3RjLv/VrbWAjc4O5L6th+3l4DqjF/\nPtxwAyxfHrpoi4jEVvDXgMxsJHAjcE5d4ZNYDwwxs9ZmZsBwYE0y7wlgWLKtXkALd3/fzDqaWbNk\n+uFAD+CtLO5KVp1ySrgrwnPPxa5ERCSzYl4DuhvYF5if2l3azDqb2VMA7r4YmAUsBZYTTuHdm6x/\nP3C4ma0EZgBfT6afBKwwsyWE03LfcPetOdqnjNPAVBEpVhqIWgD7v21bGJi6cGG4YamISEwFfwpO\nGq51aw1MFZHioxZQgez/e+9B797hiakamCoiMakFVGIOPBDOPx/uuSd2JSIimaEWUAHt/+rVoVfc\nunXhfnEiIjGoBVSC+vXbPTBVRKTQKYAKTE2X7AJquImI1EkBVGBGjIDqag1MFZHCpwAqMGYwaRL8\n7GexKxERSY86IRTg/mtgqojEpE4IJax1a/jmNzUwVUQKm1pABbr/774LRx4ZnpjasWPsakSklKgF\nVOIOOgjOO08DU0WkcKkFVMD7r4GpIhKDWkBCv34wcCDMnBm7EhGRxlMAFbhJkzQwVUQKkwKowI0Y\nEXrFTZsWuxIRkcbRNaAi2P+VK2HYMFi6FLp2jV2NiBQ7XQOSf+vfHyZODGODiiBPRaREKICKxOTJ\n8PbbMGNG7EpERBpGp+CKaP9ffhnOPBNWrIBOnWJXIyLFKlOn4BRARbb/3/0u/P3v8MgjsSsRkWKl\na0BSpylTYPlyeOyx2JWIiHw+tYCKcP9feAEuvBBWrYIDDohdjYgUG52Cy4BiDSCAa6+Fjz6CBx6I\nXYmIFBsFUAYUcwB9/HHonj11Kpx+euxqRKSY6BqQfK5994X77oOrrgotIRGRfKMWUJHv/xVXQIsW\n8Mtfxq5ERIqFTsFlQCkE0NatcNRR8LvfQUVF7GpEpBjoFJw0SPv2ofVzxRXw6aexqxER2U0toBLZ\n/0svhYMPhjvvjF2JiBQ6nYLLgFIKoC1bQq+4xx+HIUNiVyMihUyn4KRROnaEu+6C8eOhqip2NSIi\nCqCScuGF0Ls3/OAHsSsREdEpuJI5BVdj0yYYOBCefRYGD45djYgUIp2Ckybp3Bluvz2cituxI3Y1\nIlLKFEAl6LLLwvOCfvKT2JWISClr0Ck4M3sMmAbMcffqrFeVI6V4Cq7G22/DMcfA//wP9OkTuxoR\nKSS5PgU3FbgEWGtmt5lZ73S/WOLq3h1uvTWcitu1K3Y1IlKKGhRA7r7A3S8FyoF1wAIze9HMxplZ\ni2wWKNlz1VXQsiXcfXfsSkSkFDW4F5yZdQC+BowBNgLTgROB/u5eka0Cs6mUT8HVWLsWvvxleOkl\nOOKI2NWISCHI6Z0QzOxxoDfwW+ABd9+UMu9ldz8m3UJiUAAFd94JTz8NCxZAM3VLEZG9yHUADXX3\nhel+Wb5RAAW7dsHxx4frQd/4RuxqRCTf5TqAzqtj8ofASnd/N90iYlEA7bZ6dXhcw5IlcOihsasR\nkXyW6wB6GvgyUNMKqgBeAQ4D/p+7/zbdQmJQAO3pBz+Av/wlnI6ztH9aIlKsct0NuwXQx93Pd/fz\ngb6AA8cB3023CMkPN90EGzaEh9eJiGRbQwOoq7u/k/L5XeBQd/8HoBu6FIkWLeD+++GGG2Dz5tjV\niEixa2gAVZrZU2Z2mZldBvwhmdYW2Jq98iTXjj4aLr8cJk6MXYmIFLuGXgMy4DzCuB+APwOPpnMB\nxcxuB84GqoA3gXHu/lEdy10PXA5UAyuT5bab2UygV7LY/sAH7l6erDMZGA/sBL7t7vPqqUHXgOqw\nbRsMGgQ//CGMGhW7GhHJNznrhGBmzYEF7j403S+rtd1TgOfdvdrMbgPc3SfXWqYL8AJwZBI6DwNP\nu/uDtZa7A9jq7j80sz7ADOBYoCuwAOhZV9IogOr34otw/vmwahV06BC7GhHJJznrhODuu4BqM2uX\n7pfV2u6ClBubLiKERV2aA23NrAxoQ7gLQ20XEkIH4FxgprvvdPd1wFrgSxkrvEQcfzxcdBFcd13s\nSkSkWJU1cLmPgZVmNh/4pGaiu1+boTrGAzNrT3T3jWZ2J7Ae+BSY5+4LUpcxs68Am939rWTSIcBf\nUhbZkEyTRvrP/4T+/UO37DPPjF2NiBSbhgbQY8mrUZLA6pQ6idB9+2Z3fzJZ5mZgh7vPqGP99oQW\nTXfCwNdZZnZJrWUvBh5qbG01pkyZ8u/3FRUVVFRUNHVTRadtW7jvPhg7NpyKa5fRNrCIFIrKykoq\nKyszvt3G3Ix0H6Cbu7+esS83GwtMAIa5e1Ud80cBp7n7hOTzGOA4d5+YfG5OaOGUu/vGZNpNhOtJ\nP04+zwVucfeX6ti+rgE1wJVXhoGpv/pV7EpEJB/kdCCqmZ0NLAPmJp8HmdnsdL7YzEYCNwLn1BU+\nifXAEDNrnfTEGw6sSZk/AlhTEz6J2cBoM2tpZocBPYDF6dRa6n7yE3jmGXj++diViEgxaeg4oCmE\nC/lbAdx9GXB4mt99N7AvMN/MlpjZVAAz62xmTyXfsxiYBSwFlhNO4d2bso2LqHX6zd1fBR4BXgWe\nAa5WMyc97drBPffAFVfAJ5/sfXkRkYZo6DigRe4+xMyWuvvgZNoKdx+Q9QqzSKfgGudrX4OOHeGu\nu2JXIiIx5fpecKvN7BKguZn1NLO7gRfT/XIpLD//OTz8cBgjJCKSroYG0DVAP8JdCx4CPgI0QqTE\ndOgAv/hFuFXPtm2xqxGRQtfgXnDFSKfgGs893J6nd2/4r/+KXY2IxJDr5wH1Av4D+CIpY4fcfVi6\nBcSkAGqazZthwACYOxfKy2NXIyK5lusAWg7cQ3gI3a6a6e7+SroFxKQAaroHH4Sf/hT++tfwGAcR\nKR25DqBX3P3odL8s3yiAms4dzjgDTjgBvve92NWISC7lOoCmEB5C9zihIwIAyQPpCpYCKD3r14fn\nB1VWQr9+sasRkVzJdQD9vY7J7u7pDkaNSgGUvnvugd/8JnTNbt48djUikgs5DaBipQBKX3U1DB8O\nZ50VHuUtIsUvJwNRzew7Ke8vqDVPnXCFZs3g17+GH/0I1q6NXY2IFJK9DUQdnfJ+cq15IzNcixSo\nI46AW2+Fc86Bd96JXY2IFIq9PQ/I6nlf12cpYd/6Frz3HpxyCixcGO4ZJyLyefYWQF7P+7o+S4m7\n5RbYvh1GjIDnnoMDDohdkYjks8/thGBmuwiP4DZgH8JjsUk+t3b3gh6CqE4ImecON94YumYvWADt\n28euSEQyTb3gMkABlB3ucN118NJLMG8e7Ldf7IpEJJMUQBmgAMoe93BdaMWKcM+4ffeNXZGIZIoC\nKAMUQNlVXQ1XXglvvBEe6d2mTeyKRCQTFEAZoADKvupqGDcONm6E2bNhn31iVyQi6VIAZYACKDd2\n7YIxY+CDD+CJJ6BVq9gViUg6FEAZoADKnZ07YfTo0E171ixo2TJ2RSLSVDm5FY9IppSVwUMPhVv3\njB4NO3bErkhEYlMASc60aAEPPwxVVeGU3M6dsSsSkZgUQJJTrVrBo4/CP/4ROifs2rX3dUSkOCmA\nJOdatw6dETZsgAkTQk85ESk9CiCJok0bePLJMEbom98MA1dFpLQogCSatm3h6adh5Uq49lqFkEip\nUQBJVF/4AsyZA4sWhSeqKoRESocCSKJr1y7ctHThQpg8WSEkUir29jwgkZzYf3+YPx+GDQs95W69\nNXZFIpJtCiDJGx07hmcIVVSEMUPf+17sikQkmxRAklcOOig8TbWiItyu5zvfiV2RiGSLAkjyTufO\n8PzzcPLJIYSuuy52RSKSDQogyUuHHLI7hFq0CA+3E5HiogCSvNWtWwihmtNxEybErkhEMkkBJHnt\nsMPCNaGhQ0NLaOzY2BWJSKYogCTv9egRescNGxZC6NJLY1ckIpmgAJKC0Lt3GCd0yinhdNwFF8Su\nSETSpQCSgtG3L8ydC6eeGh5w99Wvxq5IRNKhAJKCMmAAPPMMnH56OB131lmxKxKRptK94KTglJfD\nU0/B+PHw7LOxqxGRplIASUE69tjwULsxY0IvOREpPAogKVjHHw+zZsHFF8Mf/xi7GhFpLAWQFLST\nToKZM2HUKPjzn2NXIyKNoQCSgjdsGEyfHnrFvfRS7GpEpKEUQFIUTj0VfvOb0CvugQf0UDuRQmBe\nwn9TzcxLef+L0bJlMG4cHHww3HsvHHpo7IpEio+Z4e6W7nbUApKiMmgQLF4MJ54Yumvfe69aQyL5\nKloLyMxuB84GqoA3gXHu/lEdy10PXA5UAyuT5bab2UygV7LY/sAH7l5uZt2BNcBrybxF7n51PTWo\nBVTEVq8OraEvfAF+/etwY1MRSV8xtIDmAf3cfRCwFphcewEz6wJcA5S7+wDCnRtGA7j7aHcvd/dy\n4FHgsZRV36iZV1/4SPHr1w9efBFGjgzjhu6+G6qrY1clIjWiBZC7L3D3mn8OFgFd61m0OdDWzMqA\nNsDGOpa5EHgo5XPaySzFoawMbrwxdNF++OHwgLu//S12VSIC+XMNaDwwp/ZEd98I3AmsBzYAW919\nQeoyZvYVYLO7v5ky+YtmtsTMFprZiVmsWwpE795hsOoFF4QBrHfcAbt2xa5KpLRl9WakZjYf6JQ6\nCXDgZnd/MlnmZmCHu8+oY/32wLlAd+BDYJaZXVJr2YvZs/WzEejm7h+YWTnwhJn1dfeP66pxypQp\n/35fUVFBRUVFo/dTCkPz5nDttXDmmXDFFfD734eu2337xq5MJL9VVlZSWVmZ8e1G7YZtZmOBCcAw\nd6+qY/4o4DR3n5B8HgMc5+4Tk8/NCS2j8qS1VNd3LARucPcldcxTJ4QSVV0desh9//tw/fXhNF2L\nFrGrEikMBd8JwcxGAjcC59QVPon1wBAza21mBgwn9HCrMQJYkxo+ZtbRzJol7w8HegBvZWMfpHA1\nawZXXQWvvBJOzR13HCxfHrsqkdIS8xrQ3cC+wPzkes1UADPrbGZPAbj7YmAWsBRYTjiFd2/KNi5i\nz9NvACcBK8xsCfAI8A1335rVPZGC1a1beMjdNdfAiBFwyy2wfXvsqkRKg+6EUML7L3vasCG0itat\nC9eGjjkmdkUi+angT8GJ5JtDDoHZs2Hy5NBR4aabYNu22FWJFC8FkEgKM7jkElixAt58M9za58UX\nY1clUpx0Cq6E91/2btascH3o4ovhhz+ENm1iVyQSn07BieTAqFGwciW8+y4MGABZGAohUrLUAirh\n/ZfGmT0brr4azj0Xbrst3ORUpBSpBSSSY+ecA6tWwb/+Bf37w7x5sSsSKWxqAZXw/kvTPfssXHll\nGDt0xx3Qvn3sikRyRy0gkYhOOy1cG2rRIrSGnn46dkUihUctoBLef8mM558PNzc94QT4+c/hgANi\nVySSXWoBieSJYcNCa+iAA+Coo+Cxx/a+joioBaQWkGTUCy/A5ZdDly6hVXTeebDPPrGrEsmsTLWA\nFEAlvP+SHVVV8Ic/wLRp8PLLMHp0CKXy8tiViWSGAigDFECSbW+/DQ88EG5uuv/+IYguvTS8FylU\nCqAMUABJrlRXw3PPhVbR3LlwxhkhjIYODc8mEikkCqAMUABJDO+/D9OnhzD65z9h3DgYOxYOPTR2\nZSINowDKAAWQxOQensg6bRo8/HB4Kuvll4c7LrRsGbs6kfopgDJAAST54tNP4dFHQxi9+ip87Wsh\njPr1i12ZyGcpgDJAAST56I034P774b//O5yWu/zy0JNONz+VfKEAygAFkOSznTtDh4Vp02DhQvjq\nV0MYnXBCeHCeSCwKoAxQAEmheOcdePDBEEYA48fD178OBx8cty4pTQqgDFAASaFxD48InzYNHn8c\nTj45tIpOPx3KymJXJ/nogw9g2TJYuhSWLAmf0715rgIoAxRAUsj++c/Qe+7++2HdOrjsstAy6tkz\ndmUSy6ZNIWSWLt0dOFu2wMCBMHhwuBtHeXn4nA4FUAYogKRYvPpqCKLf/haOPDK0ikaNgjZtYlcm\n2eAOb721Z9AsXRquG9YEzeDB4dWzZ+YHOyuAMkABJMVm+3Z46qlwiu5PfwrduAcMCP/HO2BAeO23\nX+wqpTF27oTXXtszaJYtC70iU4OmvBy6ds1NBxUFUAYogKSYbd0KK1bA8uXhzxUrwiPFDzpodyDV\n/HnEEbolUD7417/Coz1SWzarV4dgSQ2awYOhY8d4dSqAMkABJKVm1y54883doVTz55Yt4VlGqS2l\nAQOgXbvYFRevDz/cs3PA0qXhv02vXnsGzcCB+TcGTAGUAQogkWDr1vB/3qmhtGpV+L/sulpLzZvH\nrrgwVFXB5s2hc8CmTbBmze7WzebN4XHuqafRjjoKWrWKXfXeKYAyQAEkUr9du8KF7tRTeMuXw3vv\n1X1tqX372BXnzief7A6V+l4bN4aeip06QefO4dWz5+7WTa9ehRvkCqAMUACJNN6HH362tbRyJXTo\nsGcoDRwIPXoUzj+y7qEluLdg2bQJduzYHSqf9+rYsTivrSmAMkABJJIZ1dWhtZQaSsuXhzs4dOoU\n7u5d+9WqVXanp84DePfdPVsntUNl82Zo0WJ3eHTpUn+wtGtX2rdDUgBlgAJIJLs++ih0cNi+Pbyq\nqna/r/2qb15T1qk9r7o69P7bW4ulbdvYR6wwKIAyQAEkItJ4mQqgIjw7KSIihUABJCIiUSiAREQk\nCgWQiIhEoQASEZEoFEAiIhKFAkhERKJQAImISBQKIBERiUIBJCIiUSiAREQkCgWQiIhEoQASEZEo\nFEAiIhJFtAAys9vNbI2ZLTOzR81sv3qWu97MVpnZCjObbmYtk+kDzewvZrbUzBab2TEp60w2s7XJ\n9k/N1T6JiEjDxWwBzQP6ufsgYC0wufYCZtYFuAYod/cBQBkwOpl9O3CLuw8GbgF+kqzTF7gQ6AOc\nDkw1K+VnF+ZOZWVl7BKKio5n5uhY5qdoAeTuC9y9Ovm4COhaz6LNgbZmVga0ATYm06uBdsn79sCG\n5P05wEx33+nu6wjh9qUMly910F/yzNLxzBwdy/xUFruAxHhgZu2J7r7RzO4E1gOfAvPcfUEy+3rg\n2WS+Accn0w8B/pKymQ3JNBERySNZbQGZ2fzk2k3Na2Xy59kpy9wM7HD3GXWs3x44F+gOdAH2NbNL\nktnfBL7t7t0IYXR/NvdFREQyy9w93pebjQUmAMPcvaqO+aOA09x9QvJ5DHCcu080s63u3j5l2a3u\n3t7MbgLc3X+cTJ9LuFb0Uh3bj7fzIiIFzN3TvrYe7RScmY0EbgROqit8EuuBIWbWGqgChgOLk3kb\nzOxkd/+jmQ0nXOsBmA1MN7OfEU699UhZZw+ZOIAiItI0Ma8B3Q20BOYnndQWufvVZtYZuM/dz3L3\nxWY2C1gK7Ej+vC9ZfwLwCzNrDmwDrgRw91fN7BHg1WSdqz1mM09EROoU9RSciIiUrqK9E4KZjTSz\n18zsb2b23XqW+UUyYHWZmQ1qzLqlpAnHcnDK9HVmtrxmwHDuqs5fezueZtbbzF40s21mNqkx65ai\nNI+nfp+1NOB4XpIcs+Vm9oKZDWjoup/h7kX3IgTrG4Tecy2AZcCRtZY5HXg6eX8c4RRgg9YtpVc6\nxzL5/Bawf+z9yJdXA49nR+Bo4AfApMasW2qvdI5nMk+/z8YfzyFAu+T9yHT+7SzWFtCXgLXu/ra7\n7yCMMTq31jLnAg8CeOgh187MOjVw3VKSzrGEMEarWH9nTbHX4+nuW9z9FWBnY9ctQekcT9Dvs7aG\nHM9F7v5h8nERu8dZNvr3WawH/hDgf1M+/x+fHYxa3zINWbeUNOVYpg7+dUJHk7+a2YSsVVk40vl9\n6bf5WekeE/0+99TY43kFMKeJ6+bNnRDygbpkZ8cJ7r7JzA4k/EVf4+4vxC5KJKHfZxOZ2VBgHHBi\nU7dRrC2gDUC3lM9d2X2vuNRlDq1jmYasW0rSOZa4+6bkz/eAx9F9+dL5fem3+VlpHRP9Pj+jQccz\n6XhwL3COu3/QmHVTFWsA/RXoYWbdk8c3jCYMUE01G/g6gJkNAba6+zsNXLeUNPlYmlkbM9s3md4W\nOBVYlbvS81Jjf1+pLXP9Nj+rycdTv8867fV4mlk34FFgjLu/2Zh1ayvKU3DuvsvMJhIe+dAMmObu\na8zsG2G23+vuz5jZGWb2BvAJoSlZ77qRdiW6dI4l0Al4PLnlURkw3d3nxdiPfNGQ45l04HgZ+AJQ\nbWbfBvq6+8f6be4pneMJHIh+n3toyPEEvg8cwO5H3exw9y815d9ODUQVEZEoivUUnIiI5DkFkIiI\nRKEAEhGRKBRAIiIShQJIRESiUACJiEgUCiAREYlCASQiIlEU5Z0QRAqdmbUmPGulCujm7r+KXJJI\nxqkFJJKfzgb+4O5zCA8AEyk6CiCRPGNmXYE33N3N7HBgfeyaRLJBp+BE8s+A5Aav3wUGA/8RuyCR\nbFALSCTYfJTgAAAAmklEQVRPufuPgd8QbmsvUnQUQCL5J/XMxBHAP2IVIpJNCiCRPGJmHQhPkqwx\nEngsUjkiWaXnAYnkETMbDlQD+xMec/4nd18StyqR7FAnBJH80srdn4ldhEgu6BScSH7ZFbsAkVzR\nKTgREYlCLSAREYlCASQiIlEogEREJAoFkIiIRKEAEhGRKBRAIiIShQJIRESiUACJiEgU/x+ER7bG\n6lEt3gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f38a7550470>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.plot(betas, energies)\n",
    "fig[0].axes.set_xlabel(r\"$\\beta$\")\n",
    "fig[0].axes.set_ylabel(\"Energy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}