Session6_PhysicalChem.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Programming for Chemists: Physical Chemistry in Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Importance for scientists:**\n",
    "* The purpose of this session is to learn a variety of tools applicable in many different scientific applications on the job, including:\n",
    "    * Symbolic computation.\n",
    "    * Finding maxima of functions.\n",
    "    * Automatically varying fitting parameters to minimize or maximize certain values.\n",
    "\n",
    "* This hands on session will explore the hydrogen atom in intricate detail, adding more tools to your programming toolbox. \n",
    "* First we discuss a very useful Python library called [SymPy](https://www.sympy.org/en/index.html). SymPy is a Python library for symbolic mathematics; aiming to become a full-featured **computer algebra system (CAS)**. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computer Algebra in SymPy \n",
    "\n",
    "* Computer algebra refers to software packages that are capable of doing **symbolic mathematical computations** such as:\n",
    "\n",
    "$$\n",
    "    (x + 1)^2 \\hspace{0.5cm} \\rightarrow \\hspace{0.5cm} x^2 + 2x + 1,\n",
    "$$\n",
    "* where the word \"algebraic\" refers to operations with symbolic objects originating from algebra. \n",
    "* Two well known software packages allowing for this functionality are *Maple* and *Mathematica*, but Python has a library called **SymPy** which provides us the ability to do symbolic computation in Python. \n",
    "* Scientific computing is usually based on **numerical** computation using **approximate floating point numbers**, while symbolic computation emphasises **exact** computation with expressions containing variables that have no given value and are manipulated as symbols. \n",
    "\n",
    "\n",
    "| Symbolic Pros                        | Symbolic Cons                                                      |\n",
    "|:-------------------------------------|:-------------------------------------------------------------------|\n",
    "| Exact Arithmetic                     | **Much slower** than numerical computation                         |\n",
    "| Non-linear expressions with symbols  | Poor hardware support                                              |\n",
    "| Inequalities, differential equations | Difficult to parallelize                                           |\n",
    "| Results can be formulae              | Numerical problems better handled with \"orthodox\" computer methods |\n",
    "| It \"knows\" how to integrate          | Number of expressions can grow exponentially                       |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Analytical vs. Numerical Integration\n",
    "\n",
    "Consider the [known integral](https://functions.wolfram.com/Constants/Pi/07/01/01/0002/), one of the first applications of calculus by Isaac Newton (follow [this link](https://www.youtube.com/watch?v=gMlf1ELvRzc) for an interesting video on the subject)\n",
    "\n",
    "$$\n",
    "    4\\int\\limits_{0}^{1}\\sqrt{1-t^2}\\text{d}t = \\pi.\n",
    "$$\n",
    "\n",
    "* We pass this integral to SymPy and it produces \\\\(\\pi\\\\) **exactly**, not a floating point representation (we will discuss SymPy terminology throughout this tutorial):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sy\n",
    "\n",
    "# define t as symbol not requiring an explicit value\n",
    "t = sy.symbols('t')\n",
    "\n",
    "# define the integrand\n",
    "integrand = sy.sqrt(1 - t**2)\n",
    "\n",
    "# calculate the integral\n",
    "4*sy.integrate(integrand, (t, 0, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Native programming languages; Python, C++, Fortran, Java etc... **cannot do this** without the help of customised libraries, but can solve the integral **approximately** using **numerical** methods. \n",
    "* I programmed the following animation using matplotlib simulating a numerical technique called a **trapezoid Riemann sum** of the same integral above; approximating the area under the function using an increasing number of trapezoids:\n",
    "\n",
    "<center><img src=\"https://raw.githubusercontent.com/adambaskerville/ProgrammingForChemists/master/images/riemann_sum.gif\" width=\"auto\" height=\"auto\" /></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving the Hydrogen Atom "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* In this section we will conduct a detailed investigation on how to solve the hydrogen atom; the building block of a lot of computational and quantum chemical methods. \n",
    "* It can be described using the **time independent, non-relativistic Schr&ouml;dinger equation**\n",
    "\n",
    "$$\n",
    " \\hat{H}\\psi = E\\psi,\n",
    " \\tag{1}\n",
    "$$\n",
    "\n",
    "* where: \n",
    "    * \\\\(\\hat{H} \\\\) is the Hamiltonian, the total energy operator.\n",
    "    * \\\\(\\psi\\\\) is the wavefunction, containing all the key physical behaviour of the system.\n",
    "    * \\\\(E\\\\) is the energy eigenvalue. \n",
    "\n",
    "* This is the textbook system studied in quantum mechanics as it is one of the most simple, consisting of a single proton and electron. \n",
    "    * The radial part of the hydrogen atom is **exactly solvable** using a series solution method but lets pretend we are unaware of this and explore how a computer can solve it. \n",
    "    * The aim will be to get the computer to do *everything*, requiring no preliminary mathematics by us humans. \n",
    "    * We are going to use a very powerful mathematical technique called the **variational method**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Variational Method \n",
    "\n",
    "$$\n",
    "\\frac{\\int\\psi^* \\hat{H} \\psi \\text{d}\\nu}{\\int \\psi^* \\psi \\text{d}\\nu} \\ge E_g.\n",
    "\\tag{2}\n",
    "$$\n",
    "\n",
    "* If we do not know the exact wavefunction our best bet is to have an educated guess using some physical intuition; unlikely to be the exact wavefunction, but by labelling it as \\\\(\\psi_1\\\\) it will calculate a value of the energy, \\\\(E_1\\\\). \n",
    "\n",
    "* We can now guess a second wavefunction, \\\\(\\psi_2\\\\), resulting in an energy value \\\\(E_2\\\\) and if desired we could guess even more wavefunctions.\n",
    "* The variational principle tells us that if \\\\(E_g\\\\) is the exact ground state energy then both \\\\(E_1\\\\) and \\\\(E_2\\\\) will always be greater than \\\\(E_g\\\\) unless we somehow guess the exact wavefunction; which would then be equal to \\\\(E_g\\\\).\n",
    "* Due to the inability of solving many-body problems exactly, computational and quantum chemistry is built upon approximate methods and a number of methods utilize the variational principle, including **Hartree Fock theory**. \n",
    "    * These work by making repeated iterative changes to the trial wavefunctions or density matrices, calculating the energy as they go until the lowest value has been found; safe in the knowledge it can not drop below the exact ground state solution. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gaussian Wavefunction\n",
    "\n",
    "* For our first trial wavefunction, \\\\(\\psi_1\\\\), we will use a single [Gaussian function](https://mathworld.wolfram.com/GaussianFunction.html) \n",
    "\n",
    "$$\n",
    "\\psi_1 = e^{-cr^2},\n",
    "$$\n",
    "\n",
    "* where \\\\(r\\\\) is the nucleus-electron distance and \\\\(c\\\\) is a **variational parameter** which our program will change throughout the calculation in order to find the lowest energy, \\\\(E_g\\\\). \n",
    "\n",
    "* <font color='red'>**Exercise:** Plot the Gaussian function when $c=1$ using matplotlib:</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "r = np.arange(-2, 2, 0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* The Hamiltonian for the hydrogen atom in atomic units is given by\\\\(^{[1]}\\\\)\n",
    "\n",
    "$$\n",
    "\\hat{H} = -\\overbrace{\\frac{1}{2}\\nabla^2}^{K.E} - \\overbrace{\\frac{1}{r}}^{P.E}.\n",
    "\\tag{5}\n",
    "$$\n",
    "\n",
    "* The Laplacian operator for the **spherically symmetric 1s state** of the hydrogen atom is\n",
    "\n",
    "$$\n",
    "    \\nabla^2 = \\frac{\\partial^2}{\\partial r^2} + \\frac{2}{r}\\frac{\\partial}{\\partial r}.\n",
    "$$\n",
    "\n",
    "* The first thing we do is import NumPy, SciPy, SymPy and a useful module within the SciPy library called `Optimize` containing a lot of functionality relevant for non-linear optimization problems. We want to import the `minimize` function which we do using\n",
    "\n",
    "    `from scipy.optimize import minimize`\n",
    "\n",
    "* As we are going to solve this symbolically using computer algebra we need to tell Python to treat the coordinate `r` and parameter `c` as a symbol not requiring an explicit value. We do this using `sy.symbols('r')` from SymPy:\n",
    "\n",
    "    `r, c = sy.symbols('r c')`\n",
    "\n",
    "* which tells the program to treat `r` as a generic variable not needing an explicit value.\n",
    "\n",
    "* <font color='red'>**Exercise:** Import numpy, scipy, sympy and the minimize function from scipy:</font>\n",
    "* <font color='red'>**Exercise:** Tell Python to treat r and c as a symbol not a numerical value using sympy:</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Define wavefunction:**\n",
    "\n",
    "* The next step is to tell the program what our Gaussian wavefunction looks like:\n",
    "\n",
    "    `psi_1 = sy.exp(-c*r**2)`\n",
    "\n",
    "* <font color='red'>**Exercise:** program this wavefunction into the following code box:</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* where we call the **SymPy** implementation of the exponential function. \n",
    "\n",
    "**Define Hamiltonian:**\n",
    "\n",
    "* We next define the Hamiltonian operator and get SymPy to evaluate the differential operators for us, **symbolically!** \n",
    "    * Remember an operator, our Hamiltonian, only makes sense if it **is acting on something** so we will act it on the wavefunction on its right hand side of equation (2):\n",
    "\n",
    "        `hamiltonian = (-1/2)*(sy.diff(psi_1, r, r) + (2/r)*sy.diff(psi_1, r)) - (1/r)*psi_1`\n",
    "\n",
    "* There is a lot to unpack here:\n",
    "\n",
    "    1. `sy.diff` calls the differential operator from SymPy. The first differential in our \\\\(\\nabla^2\\\\) operator is a second order derivative with respect to \\\\(r\\\\) which can be programmed in two equivalent ways:\n",
    "\n",
    "        `sy.diff(psi_1, r, r)` or `sy.diff(psi_1, r, 2)`\n",
    "\n",
    "    2. the next derivative is first order:\n",
    "\n",
    "        `sy.diff(psi_1, r)`\n",
    "\n",
    "    3. we also need to multiply the potential energy term at the end by `psi_1` as it is part of the Hamiltonian operator.\n",
    "    \n",
    "* <font color='red'>**Exercise:** program the Hamiltonian into the following code box and print it. What has sympy done?</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Volume of integration:**\n",
    "\n",
    "* Quantum mechanics tells us the wavefunction gives the probability amplitude for finding the electron in a particular region of space d\\\\(\\nu\\\\)\n",
    "    * For the spherically symmetric 1s state of the hydrogen atom this region of space is a sphere surrounding the proton. \n",
    "    * The surface area of a sphere is \\\\(4\\pi r^2\\\\) which we will program as our volume element, `dv`, remember to use `sy.pi`.\n",
    "    \n",
    "* <font color='red'>**Exercise:** program the volume element into the following code box:</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Integration:**\n",
    "\n",
    "* We have now defined our wavefunction, Hamiltonian operator and volume element, so we can now conduct the integration. \n",
    "* We will implement two expressions, one for the numerator (`num`) of equation (2) and one for the denominator (`den`):\n",
    "\n",
    "    `num = sy.N(sy.integrate(psi_1*hamiltonian*dv, (r, 0, sy.oo)))`\n",
    "\n",
    "* Lets break this down:\n",
    "\n",
    "    1. To symbolically integrate in SymPy we call:\n",
    "\n",
    "        `sy.integrate(integrand, (variable to integrate over, start, end))` \n",
    "\n",
    "    2. We want to integrate over \\\\(r\\\\) from 0 to \\\\(\\infty\\\\) as the theory of quantum mechanics tells us probability densities are asymptotic at long-range; but the probability is **effectively zero** past a radius of several Bohr radii.\n",
    "\n",
    "        `(r, 0, sy.oo)`\n",
    "\n",
    "    infinity is represented as `oo` in Sympy. \n",
    "\n",
    "    3. As SymPy evaluates it symbolically, it will leave it in a general form (print out the result to check for yourself!) and we want a numeric answer so call `sy.N` which numerically evaluates the expressions for us.\n",
    "    \n",
    "* <font color='red'>**Exercise:** Program the denominator in the following code box and divide with the numerator to get the energy expression:</font>. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* **We are finished!** \n",
    "*  <font color='red'>Move your code segments into the following function and run the cell</font>. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def min_energy(c):\n",
    "    '''\n",
    "    This function minimises the energy with respect to a single variational parameter c\n",
    "    '''\n",
    "    # The scipy minimize function passes the parameters in an array. \n",
    "    # To use it within our mathematical formulae we need to extract it as a separate number, done using the index [0]\n",
    "    c = c[0]\n",
    "    \n",
    "    # define the Gaussian wavefunction\n",
    "\n",
    "    # define the Hamiltonian\n",
    "\n",
    "    # define the volume of integration\n",
    "\n",
    "    # define the numerator of the variational principle\n",
    "\n",
    "    # define the denominator of the variational principle\n",
    "\n",
    "    # calculate the energy\n",
    "    E = num/den\n",
    "    \n",
    "    # return E to the minimize function\n",
    "    return E"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* We are going to minimize this function using the **minimize** function from **scipy**:\n",
    "    1. Set an initial value for c.\n",
    "    2. Call the minimize function using the syntax:\n",
    "    \n",
    "        `minimize(term to minimize, variables to minimize, tol=convergence tolerance, bounds=boundaries for your parameters)`\n",
    "    3. Minimize will pass the variables inside a numpy array to our function which we need to extract when inside the function.\n",
    "    4. We apply boundaries in this example as otherwise scipy will make c negative which will seize up the function whilst sympy tries to solve the difficult integration.\n",
    "* <font color='red'>**Exercise:** Add the variable to be varied into the following code block and run:</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set an initial value for the variational parameter 'c'\n",
    "c = 3\n",
    "\n",
    "bnds = ((0.01,10),)\n",
    "\n",
    "# call the minimize function, tell it the parameter to vary, c, and tell it the tolerance for convergence, tol\n",
    "hydrogen_gs = minimize(min_energy, , tol=1e-6, bounds=bnds)\n",
    "\n",
    "# print the final result to terminal\n",
    "print(hydrogen_gs)\n",
    "\n",
    "# extract out the value of the variational parameter and the energy value\n",
    "print(\"\\nEnergy = {:.10f} hartrees when c = {:.6f}\".format(hydrogen_gs.fun, hydrogen_gs.x[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* We have a final energy of \\\\(E_g = -0.4244131816\\\\) hartrees when \\\\(c = 0.28294212\\\\).\n",
    "* We know the exact ground state energy for a one-electron system\n",
    "\n",
    "$$\n",
    "E_g = -\\frac{Z^2}{2},\n",
    "$$\n",
    "\n",
    "* where \\\\(Z\\\\) is the nuclear charge, so for the hydrogen atom \\\\(E = -0.5\\\\) hartrees. Using a Gaussian wavefunction produces an **Energy value with a \\\\(\\approx 15\\%\\\\) error!**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Slater-type orbital Wavefunction\n",
    "\n",
    "* Let us consider a better approximation for our wavefunction, this time using a **Slater-type orbital**:\n",
    "\n",
    "$$\n",
    "\\psi_2 = e^{-cr}.\n",
    "$$\n",
    "\n",
    " <font color='red'>**Exercise:** Add the slater function to the plot below:</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# create the figure (fig) and axes (ax) objects\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "# create r range from 0 -> 5 with a step of 0.05\n",
    "r = np.arange(0, 5, 0.05)   # start, stop, step\n",
    "\n",
    "# define gaussian function\n",
    "gaussian = np.exp(-r**2)\n",
    "\n",
    "# define slater function\n",
    "\n",
    "\n",
    "# set the range of the x-axis data\n",
    "ax.set_xlim([0, 5])\n",
    "\n",
    "# set the range of the y-axis data\n",
    "ax.set_ylim([0, 1])\n",
    "\n",
    "# set the axis labels\n",
    "ax.set_xlabel('r')\n",
    "ax.set_ylabel(r'$\\psi(r)$')\n",
    "\n",
    "# plot the Gaussian function\n",
    "ax.plot(r, gaussian, label = 'Gaussian')\n",
    "\n",
    "\n",
    "# plot the legend\n",
    "ax.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* These two forms look quite different, and the key reason why the Slater form is more **physically realistic** than the Gaussian form is that it reaches a **cusp** at the origin, \\\\(r=0\\\\), which is a rule particles have to obey when they coalesce; [rigorously proven](https://onlinelibrary.wiley.com/doi/abs/10.1002/cpa.3160100201) by Tosio Kato in 1957.\n",
    "\n",
    "* <font color='red'>**Exercise:** Return to your program, or the complete version at the bottom of the script, and change the Gaussian function for a Slater-type function and see what the new result is.</font>\n",
    "\n",
    "    * We end up with \\\\(E_g = -0.5\\\\) hartrees when \\\\(c=1\\\\), matching the exact energy to the specified precision. \n",
    "* This also reveals the nature of the exact wavefunction of the hydrogen atom but also any one-electron system\n",
    "\n",
    "$$\n",
    "\\psi = e^{-Zr},\n",
    "$$\n",
    "\n",
    "* where \\\\(Z\\\\) is the nuclear charge. \n",
    "\n",
    "* <font color='red'>**Exercise:** Complete the for loop below to plot the exact wavefunctions (r against slater) with differing nuclear charges:</font>\n",
    "* <font color='red'>**Exercise:** add the legend using `plt.legend()`:</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "import seaborn as sns\n",
    "import numpy as np\n",
    "\n",
    "# set a dark grid style for the plot (looks nice!)\n",
    "sns.set_style(\"darkgrid\")\n",
    "\n",
    "# create the figure (fig) and axes (ax) objects\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "# create r range from 0 -> 5 with a step of 0.05\n",
    "r = np.arange(0, 5, 0.05)   # start, stop, step\n",
    "\n",
    "# set the range of the x-axis data\n",
    "ax.set_xlim([0, 5])\n",
    "\n",
    "# set the range of the y-axis data\n",
    "ax.set_ylim([0, 1])\n",
    "\n",
    "# set the axis labels\n",
    "ax.set_xlabel('r')\n",
    "ax.set_ylabel(r'$\\psi(r)$')\n",
    "\n",
    "for  in range(1,5):\n",
    "    # define slater function\n",
    "    slater = \n",
    "    \n",
    "    # call the plot command\n",
    "    ax.plot(, slater, label = \"Z = {}\".format(Z))\n",
    "\n",
    "    # plot the legend\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Radial Distribution Functions\n",
    "We now have our optimised Gaussian and Slater type wavefunctions for the hydrogen atom\n",
    "\n",
    "**Gaussian:**\n",
    "\n",
    "$$\n",
    "\\psi_1(r) = e^{-0.28294212 r}\n",
    "$$\n",
    "\n",
    "**Slater:**\n",
    "\n",
    "$$\n",
    "\\psi_2(r) = e^{-1r}\n",
    "$$\n",
    "\n",
    "* We can use these wavefunctions to calculate the **radial Distribution Function (RDF)** of the hydrogen atom; i.e. the electron density around the proton\n",
    "\n",
    "$$\n",
    "    \\text{RDF} = 4\\pi r^2\\psi(r)\\psi(r)|_{r=r'}\n",
    "    \\tag{7}\n",
    "$$\n",
    "\n",
    "* where we fix the coordinate \\\\(r\\\\) to be values in the range \\\\(0 \\le r' \\le \\infty\\\\). \n",
    "* which we need only calculate **once** using SymPy, and then divide the RDF function by the normalisation value. \n",
    "* We are going to make use of the [`lambdify`](https://docs.sympy.org/latest/modules/utilities/lambdify.html) function from SymPy which transforms SymPy expressions to a numpy recognisable format meaning we can calculate numerical values very fast. \n",
    "\n",
    "* <font color='red'>**Exercise:** Add the code to calculate and plot the RDF for the slater wavefunction:</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "import sympy as sy\n",
    "from sympy import lambdify\n",
    "\n",
    "r = sy.symbols('r')\n",
    "\n",
    "# create the figure (fig) and axes (ax) objects\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "# define gaussian function\n",
    "gaussian = sy.exp(-r**2)\n",
    "\n",
    "# define slater function\n",
    "slater = sy.exp(-r)\n",
    "\n",
    "# define volume element\n",
    "dv = 4*sy.pi*r**2\n",
    "\n",
    "# calculate the inner product <psi|psi> of the gaussian wavefunction\n",
    "gaussian_norm = sy.N(sy.integrate(gaussian*gaussian*dv, (r, 0, sy.oo)))\n",
    "\n",
    "# calculate the inner product <psi|psi> of the slater-type orbital wavefunction\n",
    "slater_norm = sy.N(sy.integrate(slater*slater*dv, (r, 0, sy.oo)))\n",
    "\n",
    "# calculate the RDF using the Gaussian wavefunction\n",
    "RDF_gaussian = gaussian*gaussian*dv / gaussian_norm\n",
    "lam_gaussian = lambdify(r, RDF_gaussian) # use lambdify to then treat as a function\n",
    "\n",
    "# calculate the RDF using the slater-type orbital wavefunction\n",
    "\n",
    "# set the range of the x-axis data\n",
    "ax.set_xlim([0, 5])\n",
    "\n",
    "# set the range of the y-axis data\n",
    "ax.set_ylim([0, 1.5])\n",
    "\n",
    "# set the axis labels\n",
    "ax.set_xlabel(\"r'\")\n",
    "ax.set_ylabel(r'$4\\pi r^2\\psi(r)\\psi(r)$') # This notation uses LaTeX math font!\n",
    "\n",
    "# generate 500 r' values (rp) between 0 and 5\n",
    "rp_vals = np.linspace(0, 5, 500)  # start, stop, no. points \n",
    "\n",
    "# calculate the y-values of the graph using the RDF functions\n",
    "y_vals_gaussian = lam_gaussian(rp_vals)\n",
    "\n",
    "# call the plot command\n",
    "ax.plot(rp_vals, y_vals_gaussian, label = 'RDF Gaussian')\n",
    "\n",
    "# plot the legend\n",
    "ax.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* We can check they are normalised correctly by integrating the RDF functions to calculate the area underneath the curves.\n",
    "* <font color='red'>**Exercise:** Write the code for the slater-type orbital wavefunction:</font>:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check normalisation of the RDF using a gaussian wave function\n",
    "print(sy.N(sy.integrate(RDF_gaussian, (r, 0, sy.oo))))\n",
    "\n",
    "# check normalisation of the RDF using a gaussian wave function\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Most Probable Value of r\n",
    "\n",
    "* If we look at the RDF plots above, they are not symmetric about a point like a standard Gaussian function is; meaning the **most probable** and **average** value of \\\\(r\\\\) are not equal! \n",
    "* The most probable value is given by the position of the maximum of the curve. To find and classify extreme points of a function take the first derivative and set equal to 0 and solve. \n",
    "* As we know the function has a maximum by visual inspection we do not need to take the second derivative in order to classify whether it is a maximum, minimum or point of inspection, "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_extremes(func, arg):\n",
    "    '''\n",
    "    This function find the extreme points of a sympy function, returning the maximum value\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    func : sympy.core.mul.Mul\n",
    "           This is the symbolic sympy expression which we want to find the maximum of\n",
    "        \n",
    "    arg : sympy.core.symbol.Symbol\n",
    "          This is the variable to differentiate wirht respect to.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    max(extremes) : The maximum value in the extreme points list\n",
    "    '''\n",
    "    # calculate first derivative\n",
    "    dy = sy.diff(func, arg)\n",
    "    \n",
    "    # calculate extreme points\n",
    "    extremes = sy.solve(dy, arg)\n",
    "\n",
    "    return max(extremes)\n",
    "\n",
    "print(\"Most probable r for Gaussian = {:.3f}\".format(find_extremes(RDF_gaussian, r)))\n",
    "\n",
    "print(\"Most probable r for Slater = {:.3f}\".format(find_extremes(RDF_slater, r)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* The most probable value of \\\\(r\\\\) for hydrogen is exactly \\\\(r = 1\\\\) a.u. so this shows that the Gaussian wavefunction does not just underestimate the ground state energy, but also the physical properties of the system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Review"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Today we covered:\n",
    "\n",
    "* Differences between symbolic and numerical computation.\n",
    "* How sympy allows for symbolic mathematics to be used in solving a problem.\n",
    "* How to symbolically differentiate and integrate in sympy.\n",
    "* How to program a minimisation problem to solve the hydrogen atom.\n",
    "* How to find maxima of functions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the expectation value of the nucleus-electron distance using sympy for both the Gaussian and Slater wavefunctions. Expectation values represent the average value of a given quantity, which in the case of the proton-electron distance, \\\\(r\\\\), in the hydrogen atom is\n",
    "\n",
    "$$\n",
    "    \\langle \\psi(r)| r |\\psi(r) \\rangle = \\frac{\\int\\limits_{0}^{\\infty} \\psi(r) r \\psi(r) \\text{d}\\nu}{\\int\\limits_{0}^{\\infty}\\psi(r)\\psi(r)\\text{d}\\nu}\n",
    "    \\tag{9}\n",
    "$$\n",
    "\n",
    "**If you get stuck the answer code is at the bottom of the notebook:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Complete Hydrogen Atom Program"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import sympy as sy\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "r = sy.symbols('r')\n",
    "\n",
    "def min_energy(c):\n",
    "    '''\n",
    "    This function minimises the energy with respect to a single variational parameter c\n",
    "    \n",
    "    Parameters:\n",
    "    -----------\n",
    "    c : np.ndarray\n",
    "        This is the current value of the c parameter\n",
    "    \n",
    "    Returns:\n",
    "    --------\n",
    "    E : sympy.core.numbers.Float\n",
    "        This number represents the energy value in hartrees\n",
    "\n",
    "    '''\n",
    "    # The scipy minimize function passes the parameters in an array. \n",
    "    # To use it within our mathematical formulae we need to extract it as a separate number, done using the index [0]\n",
    "    c = c[0]\n",
    "\n",
    "    # define the Gaussian wavefunction\n",
    "    psi = sy.exp(-c*r**2)\n",
    "\n",
    "    # define the Hamiltonian\n",
    "    hamiltonian = (-1/2)*(sy.diff(psi, r, r) + (2/r)*sy.diff(psi, r)) - (1/r)*psi\n",
    "\n",
    "    # define the volume of integration\n",
    "    dv = 4*sy.pi*r**2\n",
    "\n",
    "    # define the numerator of the variational principle\n",
    "    num = sy.N(sy.integrate(psi*hamiltonian*dv, (r, 0, sy.oo)))\n",
    "\n",
    "    # define the denominator of the variational principle\n",
    "    den = sy.N(sy.integrate(psi*psi*dv, (r, 0, sy.oo)))\n",
    "\n",
    "    # calculate the energy\n",
    "    E = num/den\n",
    "\n",
    "    # return E to the minimize function\n",
    "    return E\n",
    "\n",
    "# set an initial value for the variational parameter 'c'\n",
    "c = 3\n",
    "\n",
    "# define boundaries on the parameter c\n",
    "bnds = ((0.01,10),)\n",
    "\n",
    "# call the minimize function, tell it the parameter to vary, c, and tell it the tolerance for convergence\n",
    "hydrogen_gs = minimize(min_energy, c, tol=1e-6, bounds=bnds)\n",
    "\n",
    "print(hydrogen_gs)\n",
    "\n",
    "print(\"\\nEnergy = {:.10f} hartrees when c = {:.6f}\".format(hydrogen_gs.fun, hydrogen_gs.x[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Expectation value answer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define gaussian function\n",
    "gaussian = sy.exp(-r**2)\n",
    "\n",
    "# define gaussian function\n",
    "slater = sy.exp(-r)\n",
    "\n",
    "# define the volume of integration\n",
    "dv = 4*sy.pi*r**2\n",
    "    \n",
    "# gaussian expectation value of r\n",
    "expec_num_gaussian = sy.N(sy.integrate(gaussian*r*gaussian*dv, (r, 0, sy.oo)))\n",
    "expec_den_gaussian = sy.N(sy.integrate(gaussian*gaussian*dv, (r, 0, sy.oo)))\n",
    "\n",
    "print(\"<r> using Gaussian = {:.4f}\".format(expec_num_gaussian / expec_den_gaussian))\n",
    "\n",
    "# slater-type expectation value of r\n",
    "expec_num_slater = sy.N(sy.integrate(slater*r*slater*dv, (r, 0, sy.oo)))\n",
    "expec_den_slater = sy.N(sy.integrate(slater*slater*dv, (r, 0, sy.oo)))\n",
    "\n",
    "print(\"<r> using Slater = {:.4f}\".format(expec_num_slater / expec_den_slater))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Further Reading\n",
    "\n",
    "* If you found this example interesting, then I have written similar posts on my T>T blog:\n",
    "    * A [very comprehensive guide to programming Hartree Fock in Python](https://adambaskerville.github.io/posts/HartreeFockGuide/). \n",
    "    * A guide on how to program a gas molecule following [Brownian motion](https://adambaskerville.github.io/posts/BrownianMotion/). \n",
    "    * Two guides on how to program a Langton's ant Turing machine; which produces order (and fractals!) from seemingly random movement patterns. Langton's ant is also one of the challenge problems from Project Euler:\n",
    "\n",
    "    [Langton's ant](https://adambaskerville.github.io/posts/LangtonsAnt/)\n",
    "\n",
    "    [Extension to Langton's ant to produce fractals](https://adambaskerville.github.io/posts/LangtonsAnt2ElectricBoogaloo/)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# References\n",
    "\n",
    "[1] Quantum Mechanics, F. Mandl, 1992, Wiley & Sons, page 187."
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}