ipynb/ParticleIn1DBox.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#The-One-Dimensional-Particle-in-a-Box\" data-toc-modified-id=\"The-One-Dimensional-Particle-in-a-Box-3\"><span class=\"toc-item-num\">3&nbsp;&nbsp;</span>The One-Dimensional Particle in a Box</a></span><ul class=\"toc-item\"><li><span><a href=\"#🥅-Learning-Objectives\" data-toc-modified-id=\"🥅-Learning-Objectives-3.1\"><span class=\"toc-item-num\">3.1&nbsp;&nbsp;</span>🥅 Learning Objectives</a></span></li><li><span><a href=\"#Cyanine-Dyes\" data-toc-modified-id=\"Cyanine-Dyes-3.2\"><span class=\"toc-item-num\">3.2&nbsp;&nbsp;</span>Cyanine Dyes</a></span></li><li><span><a href=\"#Postulate:-The-squared-magnitude-of-the-wavefunction-is-proportional-to-probability\" data-toc-modified-id=\"Postulate:-The-squared-magnitude-of-the-wavefunction-is-proportional-to-probability-3.3\"><span class=\"toc-item-num\">3.3&nbsp;&nbsp;</span>Postulate: The squared magnitude of the wavefunction is proportional to probability</a></span></li><li><span><a href=\"#The-Wavefunctions-of-the-Particle-in-a-Box-(boundary-conditions)\" data-toc-modified-id=\"The-Wavefunctions-of-the-Particle-in-a-Box-(boundary-conditions)-3.4\"><span class=\"toc-item-num\">3.4&nbsp;&nbsp;</span>The Wavefunctions of the Particle in a Box (boundary conditions)</a></span></li><li><span><a href=\"#Postulate:-The-wavefunction-of-a-system-is-determined-by-solving-the-Schrödinger-equation\" data-toc-modified-id=\"Postulate:-The-wavefunction-of-a-system-is-determined-by-solving-the-Schrödinger-equation-3.5\"><span class=\"toc-item-num\">3.5&nbsp;&nbsp;</span>Postulate: The wavefunction of a system is determined by solving the Schrödinger equation</a></span></li><li><span><a href=\"#The-Wavefunctions-of-the-Particle-in-a-Box-(solution)\" data-toc-modified-id=\"The-Wavefunctions-of-the-Particle-in-a-Box-(solution)-3.6\"><span class=\"toc-item-num\">3.6&nbsp;&nbsp;</span>The Wavefunctions of the Particle in a Box (solution)</a></span></li><li><span><a href=\"#Normalization-of-Wavefunctions\" data-toc-modified-id=\"Normalization-of-Wavefunctions-3.7\"><span class=\"toc-item-num\">3.7&nbsp;&nbsp;</span>Normalization of Wavefunctions</a></span></li><li><span><a href=\"#Normalization-Check\" data-toc-modified-id=\"Normalization-Check-3.8\"><span class=\"toc-item-num\">3.8&nbsp;&nbsp;</span>Normalization Check</a></span></li><li><span><a href=\"#The-Energies-of-the-Particle-in-a-Box\" data-toc-modified-id=\"The-Energies-of-the-Particle-in-a-Box-3.9\"><span class=\"toc-item-num\">3.9&nbsp;&nbsp;</span>The Energies of the Particle in a Box</a></span></li><li><span><a href=\"#Zero-Point-Energy\" data-toc-modified-id=\"Zero-Point-Energy-3.10\"><span class=\"toc-item-num\">3.10&nbsp;&nbsp;</span>Zero-Point Energy</a></span></li><li><span><a href=\"#Atomic-Units\" data-toc-modified-id=\"Atomic-Units-3.11\"><span class=\"toc-item-num\">3.11&nbsp;&nbsp;</span>Atomic Units</a></span></li><li><span><a href=\"#Postulate:-The-wavefunction-contains-all-the-physically-meaningful-information-about-a-system.\" data-toc-modified-id=\"Postulate:-The-wavefunction-contains-all-the-physically-meaningful-information-about-a-system.-3.12\"><span class=\"toc-item-num\">3.12&nbsp;&nbsp;</span>Postulate: The wavefunction contains all the physically meaningful information about a system.</a></span></li><li><span><a href=\"#Postulate:-Observable-Quantities-Correspond-to-Linear,-Hermitian-Operators.\" data-toc-modified-id=\"Postulate:-Observable-Quantities-Correspond-to-Linear,-Hermitian-Operators.-3.13\"><span class=\"toc-item-num\">3.13&nbsp;&nbsp;</span>Postulate: Observable Quantities Correspond to Linear, Hermitian Operators.</a></span></li><li><span><a href=\"#Heisenberg-Uncertainty-Principle\" data-toc-modified-id=\"Heisenberg-Uncertainty-Principle-3.14\"><span class=\"toc-item-num\">3.14&nbsp;&nbsp;</span>Heisenberg Uncertainty Principle</a></span></li><li><span><a href=\"#Double-Checking-the-Energy-of-a-Particle-in-a-Box\" data-toc-modified-id=\"Double-Checking-the-Energy-of-a-Particle-in-a-Box-3.15\"><span class=\"toc-item-num\">3.15&nbsp;&nbsp;</span>Double-Checking the Energy of a Particle-in-a-Box</a></span></li><li><span><a href=\"#Visualizing-the-Particle-in-a-Box-Wavefunctions,-Probabilities,-etc.\" data-toc-modified-id=\"Visualizing-the-Particle-in-a-Box-Wavefunctions,-Probabilities,-etc.-3.16\"><span class=\"toc-item-num\">3.16&nbsp;&nbsp;</span>Visualizing the Particle-in-a-Box Wavefunctions, Probabilities, etc.</a></span></li><li><span><a href=\"#🪞-Self-Reflection\" data-toc-modified-id=\"🪞-Self-Reflection-3.17\"><span class=\"toc-item-num\">3.17&nbsp;&nbsp;</span>🪞 Self-Reflection</a></span></li><li><span><a href=\"#🤔-Thought-Provoking-Questions\" data-toc-modified-id=\"🤔-Thought-Provoking-Questions-3.18\"><span class=\"toc-item-num\">3.18&nbsp;&nbsp;</span>🤔 Thought-Provoking Questions</a></span></li><li><span><a href=\"#🔁-Recapitulation\" data-toc-modified-id=\"🔁-Recapitulation-3.19\"><span class=\"toc-item-num\">3.19&nbsp;&nbsp;</span>🔁 Recapitulation</a></span></li><li><span><a href=\"#🔮-Next-Up...\" data-toc-modified-id=\"🔮-Next-Up...-3.20\"><span class=\"toc-item-num\">3.20&nbsp;&nbsp;</span>🔮 Next Up...</a></span></li><li><span><a href=\"#📚-References\" data-toc-modified-id=\"📚-References-3.21\"><span class=\"toc-item-num\">3.21&nbsp;&nbsp;</span>📚 References</a></span></li></ul></li></ul></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The One-Dimensional Particle in a Box\n",
    "![Cyanine Dyes](https://github.com/PaulWAyers/IntroQChem/blob/main/linkedFiles/cyanine.jfif?raw=true \"The absorption spectrum of cyanine dyes can be qualitatively described by a particle-in-a-box model. Image from Jeff Canaria, University of Rhode Island\")\n",
    "\n",
    "## &#x1f945; Learning Objectives\n",
    "- Determine the energies and eigenfunctions of the particle-in-a-box.\n",
    "- Learn how to normalize a wavefunction.\n",
    "- Learn how to compute expectation values for quantum-mechanical operators.\n",
    "- Learn the postulates of quantum mechanics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cyanine Dyes\n",
    "![Cyanine Dyes](https://github.com/PaulWAyers/IntroQChem/blob/main/linkedFiles/CyanineDyes.svg?raw=true? \"Resonance structures of two different types of cyanine dye\")\n",
    "Cyanine dye molecules are often modelled as one-dimension particles in a box. To understand why, start by thinking classically.  You learn in organic chemistry that electrons can more “freely” along alternating double bonds.  If this is true, then you can imagine that the electrons can more from one Nitrogen to the other, almost without resistance.  On the other hand, there are sp<sup>3</sup>-hybridized functional groups attached to the Nitrogen atom, so once the electron gets to Nitrogen atom, it has to turn around and go back whence it came.   A very, very, very simple model would be to imagine that the electron is totally free between the Nitrogen atoms, and totally forbidden from going much beyond the Nitrogen atoms.  This suggests modeling these systems a potential energy function like:\n",
    "$$\n",
    "V(x) = \n",
    "\\begin{cases}\n",
    "    +\\infty & x\\leq 0\\\\\n",
    "    0       & 0\\lt x \\lt a\\\\\n",
    "    +\\infty & a \\leq x\n",
    "\\end{cases}\n",
    "$$\n",
    "where $a$ is the length of the box. A reasonable approximate formula for $a$ is\n",
    "$$\n",
    "a = \\left(5.67 + 2.49 (k + 1)\\right) \\cdot 10^{-10} \\text{ m}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Postulate: The squared magnitude of the wavefunction is proportional to probability\n",
    "What is the interpretation of the wavefunction? The Born postulate indicates that the squared magnitude of the wavefunction is proportional to the probability of observing the system at that location. E.g., if $\\psi(x)$ is the wavefunction for an electron as a function of $x$, then\n",
    "$$ p(x) = |\\psi(x)|^2 $$\n",
    "is the probability of observing an electron at the point $x$. This is called the Born Postulate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Wavefunctions of the Particle in a Box (boundary conditions)\n",
    "![Particle in a Box Potential](https://upload.wikimedia.org/wikipedia/commons/thumb/1/13/Infinite_potential_well-en.svg/308px-Infinite_potential_well-en.svg.png \"Particle in a box potential, licensed CC-SA3 by Papa.\")\n",
    "\n",
    "The nice thing about this “particle in a box” model is that it is easy to solve the time-independent Schrödinger equation in this case.  Because there is no chance that the particle could ever “escape” an infinite box like this (such an electron would have infinite potential energy!), $|\\psi(x)|^2$ must equal zero outside the box. Therefore the wavefunction can only be nonzero inside the box. In addition, the wavefunction should be zero at the edges of the box, because otherwise the wavefunction will not be continuous.  So we should have a wavefunction like\n",
    "$$\n",
    "\\psi(x) = \n",
    "\\begin{cases}\n",
    "    0& x\\le 0\\\\\n",
    "    \\text{?????}    & 0 < x < a\\\\\n",
    "    0 & a \\le x\n",
    "\\end{cases}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Postulate: The wavefunction of a system is determined by solving the Schr&ouml;dinger equation\n",
    "How do we find the wavefunction for the particle-in-a-box or, for that matter, any other system? The wavefunction can be determined by solving the time-independent (when the potential is time-independent) or time-dependent (when the potential is time-dependent) Schr&ouml;dinger equation. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Wavefunctions of the Particle in a Box (solution)\n",
    "![Wavefunction gif](https://upload.wikimedia.org/wikipedia/commons/8/8f/InfiniteSquareWellAnimation.gif \"The wavefunctions of the 1-dimensional particle in a box, both classically and quantum-mechanically\")\n",
    "To find the wavefunctions for a system, one solves the Schr&ouml;dinger equation. For a particle of mass $m$ in a one-dimensional box, the (time-independent) Schr&ouml;dinger equation is:\n",
    "$$ \n",
    "\\left(-\\frac{\\hbar^2}{2m} \\frac{d^2}{dx^2} + V(x) \\right)\\psi_n(x) = E_n \\psi_n(x)\n",
    "$$\n",
    "where \n",
    "$$\n",
    "V(x) = \n",
    "\\begin{cases}\n",
    "    +\\infty & x\\leq 0\\\\\n",
    "    0       & 0\\lt x \\lt a\\\\\n",
    "    +\\infty & a \\leq x\n",
    "\\end{cases}\n",
    "$$\n",
    "We already deduced that $\\psi(x) = 0$ except when the electron is inside the box ($0 < x < a$), so we only need to consider the Schr&ouml;dinger equation inside the box: \n",
    "$$ \n",
    "\\left(-\\frac{\\hbar^2}{2m} \\frac{d^2}{dx^2} \\right)\\psi_n(x) = E_n \\psi_n(x)\n",
    "$$\n",
    "There are systematic ways to solve this equation, but let's solve it by inspection. That is, we need to know:\n",
    "> Question: What function(s), when differentiated twice, are proportional to themselves?\\\n",
    "\n",
    "This suggests that the eigenfunctions of the 1-dimensional particle-in-a-box must be some linear combination of sine and cosine functions, \n",
    "$$\n",
    "\\psi_n(x) = A \\sin(cx) + B  \\cos(dx)\n",
    "$$\n",
    "We know that the wavefunction must be zero at the edges of the box, $\\psi(0) = 0$ and $\\psi(a) = 0$. These are called the *boundary conditions* for the problem. Examining the first boundary condition,\n",
    "$$\n",
    "0 = \\psi(0) = A \\sin(0) + B \\cos(0) = 0 + B\n",
    "$$\n",
    "indicates that $B=0$. The second boundary condition \n",
    "$$\n",
    "0 = \\psi(a) = A \\sin(ca)\n",
    "$$\n",
    "requires us to recall that $\\sin(x) = 0$ whenever $x$ is an integer multiple of $\\pi$. So $c=n\\pi$ where $n=1,2,3,\\ldots$. The wavefunction for the particle in a box is thus,\n",
    "$$\n",
    "\\psi_n(x) = A_n \\sin\\left(\\tfrac{n \\pi x}{a}\\right) \\qquad \\qquad n=1,2,3,\\ldots\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Normalization of Wavefunctions\n",
    "As seen in the previous section, if a wavefunction solves the Schr&ouml;dinger equation, any constant multiple of the wavefunction also solves the Schr&ouml;dinger equation,\n",
    "$$\n",
    "\\hat{H} \\psi(x) = E \\psi(x) \\quad \\longleftrightarrow \\quad \\hat{H} \\left(A\\psi(x)\\right) = E \\left(A\\psi(x)\\right)\n",
    "$$\n",
    "Owing to the Born postulate, the complex square of the wavefunction can be interpreted as probability. Since the probability of a particle being at *some* point in space is one, we can define the normalization constant, $A$, for the wavefunction through the requirement that:\n",
    "$$\n",
    "\\int_{-\\infty}^{\\infty} \\left|\\psi(x)\\right|^2 dx = 1.\n",
    "$$\n",
    "In the case of a particle in a box, this is:\n",
    "\\begin{align}\n",
    "1 &= \\int_{-\\infty}^{\\infty} \\left|\\psi_n(x)\\right|^2 dx  \\\\\n",
    "&= \\int_0^a \\psi_n(x) \\psi_n^*(x) dx \\\\\n",
    "&= \\int_0^a A_n \\sin\\left(\\tfrac{n \\pi x}{a}\\right) \\left(A_n \\sin\\left(\\tfrac{n \\pi x}{a}\\right) \\right)^* dx \\\\\n",
    "&= \\left|A_n\\right|^2\\int_0^a \\sin^2\\left(\\tfrac{n \\pi x}{a}\\right) dx\n",
    "\\end{align}\n",
    "To evaluate this integral, it is useful to remember some [trigonometric identities](https://en.wikipedia.org/wiki/List_of_trigonometric_identities). (You can learn more about how I remember trigonometric identities [here](https://github.com/PaulWAyers/IntroQChem/blob/main/linkedFiles/TrigIdentities.md).) The specific identity we need here is $\\sin^2 x = \\tfrac{1}{2}(1-\\cos 2x)$:\n",
    "\n",
    "\\begin{align}\n",
    "1 &= \\left|A_n\\right|^2\\int_0^a \\sin^2\\left(\\tfrac{n \\pi x}{a}\\right) \\,dx \\\\\n",
    " &= \\left|A_n\\right|^2\\int_0^a \\tfrac{1}{2}\\left(1-\\cos \\left(\\tfrac{2n \\pi x}{a}\\right)\\right) \\,dx \\\\\n",
    " &=\\tfrac{\\left|A_n\\right|^2}{2} \\left( \\int_0^a 1 \\,dx - \\int_0^a \\cos \\left(\\tfrac{2n \\pi x}{a}\\right)\\,dx  \\right) \\\\\n",
    " &=\\tfrac{\\left|A_n\\right|^2}{2} \\left( \\left[ x \\right]_0^a - \\left[\\frac{-a}{2 n \\pi}\\sin \\left(\\tfrac{2n \\pi x}{a}\\right) \\right]_0^a  \\right) \\\\\n",
    " &=\\tfrac{\\left|A_n\\right|^2}{2} \\left( a - 0  \\right) \n",
    "\\end{align}\n",
    "So \n",
    "$$\n",
    " \\left|A_n\\right|^2 = \\tfrac{2}{a} \n",
    "$$\n",
    "Note that this does not completely determine $A_n$. For example, any of the following normalization constants are allowed,\n",
    "$$\n",
    "A_n = \\sqrt{\\tfrac{2}{a}}\n",
    "= - \\sqrt{\\tfrac{2}{a}}\n",
    "= i \\sqrt{\\tfrac{2}{a}}\n",
    "= -i \\sqrt{\\tfrac{2}{a}}\n",
    "$$\n",
    "In general, any [square root of unity](https://en.wikipedia.org/wiki/Root_of_unity) can be used, \n",
    "$$\n",
    "A_n = \\left( k \\pm i \\sqrt{1-k^2} \\right) \\sqrt{\\tfrac{2}{a}}\n",
    "$$\n",
    "where $k$ is any real number. The arbitrariness of the *phase* of the wavefunction is an important feature. Because the wavefunction can be imaginary (e.g., if you choose $A_n = i \\sqrt{\\tfrac{2}{a}}$), it is obvious that the wavefunction is not an observable property of a system. **The wavefunction is only a mathematical tool for quantum mechanics; it is not a physical object.**\n",
    "\n",
    "Summarizing, the (normalized) wavefunction for a particle with mass $m$ confined to a one-dimensional box with length $a$ can be written as:\n",
    "$$\n",
    "\\psi_n(x) = \\sqrt{\\tfrac{2}{a}} \\sin\\left(\\tfrac{n \\pi x}{a}\\right) \\qquad \\qquad n=1,2,3,\\ldots\n",
    "$$\n",
    "\n",
    "Note that in this case, the normalization condition is the same for all $n$; that is an unusual property of the particle-in-a-box wavefunction. \n",
    "\n",
    "While this normalization convention is used 99% of the time, there are some cases where it is more convenient to make a different choice for the amplitude of the wavefunctions. I say this to remind you that normalization the wavefunction is something we do for convenience; it is not required by physics!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Normalization Check\n",
    "One advantage of using Jupyter is that we can easily check our (symbolic) mathematics. Let's confirm that the wavefunction is normalized by evaluating,\n",
    "$$\n",
    "\\int_0^a \\left| \\psi_n(x) \\right|^2 \\, dx\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute this code block to import required objects.\n",
    "\n",
    "# Note: The numpy library from autograd is imported, which behaves the same as\n",
    "#       importing numpy directly. However, to make automatic differentiation work,\n",
    "#       you should NOT import numpy directly by `import numpy as np`.\n",
    "\n",
    "import autograd.numpy as np\n",
    "from autograd import elementwise_grad as egrad\n",
    "\n",
    "# import numpy as np\n",
    "from scipy.integrate import trapezoid, quad\n",
    "from scipy import constants\n",
    "\n",
    "import ipywidgets as widgets\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# set the size of the plot\n",
    "# plt.rcParams['figure.figsize'] = [10, 5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3dd3453b61924d02a4739312452be4e0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(VBox(children=(FloatSlider(value=1.0, description='a (length):', max=10.0, min=0.01, step=0.01)…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#The next few lines just set up the sliders for setting parameters.\n",
    "#Principle quantum number slider:\n",
    "n = widgets.IntSlider(\n",
    "    value=1,\n",
    "    min=1,\n",
    "    max=10,\n",
    "    step=1,\n",
    "    description='n (quant. #):',\n",
    "    disabled=False,\n",
    "    continuous_update=True,\n",
    "    orientation='horizontal',\n",
    "    readout=True,\n",
    "    readout_format='d')\n",
    "\n",
    "#Box length slider:\n",
    "a = widgets.FloatSlider(\n",
    "    value=1,\n",
    "    min=.01,\n",
    "    max=10.0,\n",
    "    step=0.01,\n",
    "    description='a (length):',\n",
    "    disabled=False,\n",
    "    continuous_update=True,\n",
    "    orientation='horizontal',\n",
    "    readout=True,\n",
    "    readout_format='.2f',\n",
    ")\n",
    "\n",
    "# Define a function for the wavefunction\n",
    "def compute_wavefunction(x, n, a):\n",
    "    \"\"\"Compute 1-dimensional particle-in-a-box wave-function value(s).\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    x: float or np.ndarray\n",
    "        Position of the particle.\n",
    "    n: int\n",
    "        Quantum number value.\n",
    "    a: float \n",
    "        Length of the box.\n",
    "    \"\"\"\n",
    "    # check argument n\n",
    "    if not (isinstance(n, int) and n > 0):\n",
    "        raise ValueError(\"Argument n should be a positive integer.\")\n",
    "    # check argument a\n",
    "    if a <= 0.0:\n",
    "        raise ValueError(\"Argument a should be positive.\")\n",
    "    # check argument x\n",
    "    if not (isinstance(x, float) or hasattr(x, \"__iter__\")):\n",
    "        raise ValueError(\"Argument x should be a float or an array!\")\n",
    "        \n",
    "    # compute wave-function\n",
    "    value = np.sqrt(2 / a) * np.sin(n * np.pi * x / a)\n",
    "\n",
    "    # set wave-function values out of the box equal to zero\n",
    "    if hasattr(x, \"__iter__\"):\n",
    "        value[x > a] = 0.0\n",
    "        value[x < 0] = 0.0\n",
    "    else:\n",
    "        if x < 0.0 or x > a:\n",
    "            value = 0.0\n",
    "    return value\n",
    "\n",
    "# Define a function for the wavefunction squared\n",
    "def compute_probability(x, n, a):\n",
    "    \"\"\"Compute 1-dimensional particle-in-a-box probablity value(s).\n",
    "    \n",
    "    See `compute_wavefunction` parameters.\n",
    "    \"\"\"\n",
    "    return compute_wavefunction(x, n, a)**2\n",
    "\n",
    "#This next bit of code just prints out the normalization error\n",
    "def check_normalization(a, n):\n",
    "    #check the computed values of the moments against the analytic formula\n",
    "    normalization,error = quad(compute_probability, 0, a, args=(n, a))\n",
    "    print(\"Normalization of wavefunction = \", normalization)\n",
    "\n",
    "\n",
    "out = widgets.interactive_output(check_normalization, {'a': a, 'n': n})\n",
    "\n",
    "widgets.HBox([widgets.VBox([a, n]),out])  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Energies of the Particle in a Box\n",
    "How do we compute the energy of a particle in a box? All we need to do is substitute the eigenfunctions of the Hamiltonian, $\\psi_n(x)$ back into the Schr&ouml;dinger equation to determine the eigenenergies, $E_n$. That is, from\n",
    "$$\n",
    "\\hat{H} \\psi_n(x) = E_n \\psi_n(x) \n",
    "$$\n",
    "we deduce\n",
    "\\begin{align}\n",
    "-\\frac{\\hbar^2}{2m} \\frac{d^2}{dx^2} \\left( A_n \\sin \\left( \\frac{n \\pi x}{a}\\right) \\right)\n",
    " &= E_n  \\left( A_n \\sin \\left( \\frac{n \\pi x}{a}\\right) \\right) \\\\\n",
    " -A_n \\frac{\\hbar^2}{2m} \\frac{d}{dx} \\left( \\frac{n \\pi}{a} \\cos \\left( \\frac{n \\pi x}{a}\\right) \\right)\n",
    " &= E_n  \\left( A_n \\sin \\left( \\frac{n \\pi x}{a}\\right) \\right) \\\\\n",
    "  A_n \\frac{\\hbar^2}{2m} \\left( \\frac{n \\pi}{a} \\right)^2 \\sin \\left( \\frac{n \\pi x}{a}\\right) \n",
    " &= E_n  \\left( A_n \\sin \\left( \\frac{n \\pi x}{a}\\right) \\right) \\\\\n",
    "   \\frac{\\hbar^2 n^2 \\pi^2}{2ma^2}  \n",
    " &= E_n   \n",
    "\\end{align}\n",
    "Using the definition of $\\hbar$, we can rearrange this to:\n",
    "\\begin{align}\n",
    "E_n &= \\frac{\\hbar^2 n^2 \\pi^2}{2ma^2} \\qquad \\qquad n=1,2,3,\\ldots\\\\\n",
    "&= \\frac{h^2 n^2}{8ma^2}\n",
    "\\end{align}\n",
    "Notice that only certain energies are allowed. This is a fundamental principle of quantum mechanics, and it is related to the \"waviness\" of particles. Certain \"frequencies\" are resonant, and other \"frequencies\" cannot be observed. *The **only** energies that can be observed for a particle-in-a-box are the ones given by the above formula.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Zero-Point Energy\n",
    "Na&iuml;vely, you might expect that the lowest-energy state of a particle in a box has zero energy. (The potential in the box is zero, after all, so shouldn't the lowest-energy state be the state with zero kinetic energy? And if the kinetic energy were zero and the potential energy were zero, then the total energy would be zero.) \n",
    "\n",
    "But this doesn't happen. It turns out that you can never \"stop\" a quantum particle; it always has a zero-point motion, typically a resonant oscillation about the lowest-potential-energy location(s). Indeed, the more you try to confine a particle to stop it, the bigger its kinetic energy becomes. This is clear in the particle-in-a-box, which has only kinetic energy. There the (kinetic) energy increases rapidly, as $a^{-2}$, as the box becomes smaller:\n",
    "$$ \n",
    "T_n = E_n = \\frac{h^2n^2}{8ma^2}\n",
    "$$\n",
    "The residual energy in the electronic ground state is called the **zero-point energy**, \n",
    "$$\n",
    "E_{\\text{zero-point energy}} = \\frac{h^2}{8ma^2}\n",
    "$$\n",
    "The existence of the zero-point energy, and the fact that zero-point kinetic energy is always positive, is a general feature of quantum mechanics. \n",
    "> **Zero-Point Energy Principle:** Let $V(x)$ be a nonnegative potential. The ground-state energy is always greater than zero.\n",
    "\n",
    "More generally, for any potential that is bound from below,\n",
    "$$\n",
    "V_{\\text{min}}= \\min_x V(x)\n",
    "$$\n",
    "the ground-state energy of the system satisfies $E_{\\text{zero-point energy}} > V_{\\text{min}}$. \n",
    ">Nuance: There is a tiny mathematical footnote here; there are some $V(x)$ for which there are *no* bound states. In such cases, e.g., $V(x) = \\text{constant}$, it is possible for $E = V_{\\text{min}}$.)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Atomic Units\n",
    "Because Planck's constant and the electron mass are tiny numbers, it is often useful to use [atomic units](https://en.wikipedia.org/wiki/Hartree_atomic_units) when performing calculations. We'll learn more about atomic units later but, for now, we only need to know that, in atomic units, $\\hbar$, the mass of the electron, $m_e$, the charge of the electron, $e$, and the average (mean) distance of an electron from the nucleus in the Hydrogen atom, $a_0$, are all defined to be equal to 1.0 in atomic units. \n",
    "\\begin{align}\n",
    "\\hbar &= \\frac{h}{2 \\pi} = 1.0545718 \\times 10^{-34} \\text{ J s}= 1 \\text{ a.u.} \\\\\n",
    "m_e &= 9.10938356 \\times 10^{-31} \\text{ kg} = 1 \\text{ a.u.} \\\\\n",
    "e &= 1.602176634 \\times 10^-19 \\text{ C} = 1  \\text{ a.u.} \\\\\n",
    "a_0 &= \t5.291772109 \\times 10^{-11} \\text{ m} = 1 \\text{ Bohr} = 1 \\text{ a.u.}\n",
    "\\end{align}\n",
    "The unit of energy in atomic units is called the Hartree, \n",
    "$$\n",
    "E_h =4.359744722 \\times 10^{-18} \\text{ J} = 1 \\text{ Hartree} = 1 \\text{ a.u.}\n",
    "$$\n",
    "and the ground-state (zero-point) energy of the Hydrogen atom is $-\\tfrac{1}{2} E_h$. \n",
    "\n",
    "We can now define functions for the eigenenergies of the 1-dimensional particle in a box:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "code_folding": []
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8bbf39a36b7f4ba6a4c33c63f15ad021",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(VBox(children=(FloatSlider(value=1.0, description='a (length):', max=10.0, min=0.01, step=0.01)…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#The next few lines just set up the sliders for setting parameters.\n",
    "#Principle quantum number slider:\n",
    "n = widgets.IntSlider(\n",
    "    value=1,\n",
    "    min=1,\n",
    "    max=10,\n",
    "    step=1,\n",
    "    description='n (quant. #):',\n",
    "    disabled=False,\n",
    "    continuous_update=True,\n",
    "    orientation='horizontal',\n",
    "    readout=True,\n",
    "    readout_format='d')\n",
    "\n",
    "#Box length slider:\n",
    "a = widgets.FloatSlider(\n",
    "    value=1,\n",
    "    min=.01,\n",
    "    max=10.0,\n",
    "    step=0.01,\n",
    "    description='a (length):',\n",
    "    disabled=False,\n",
    "    continuous_update=True,\n",
    "    orientation='horizontal',\n",
    "    readout=True,\n",
    "    readout_format='.2f',\n",
    ")\n",
    "\n",
    "# Define a function for the energy of a particle in a box\n",
    "# with length a and quantum number n [in atomic units!]\n",
    "# The length is input in Bohr (atomic units)\n",
    "def compute_energy(n, a):\n",
    "    \"Compute 1-dimensional particle-in-a-box energy.\"\n",
    "    return n**2 * np.pi**2 / (2 * a**2)\n",
    "\n",
    "# Define a function for the energy of an electron in a box\n",
    "# with length a and quantum number n [in SI units!]. \n",
    "# The length is input in meters.\n",
    "def compute_energy_si(n, a):\n",
    "    \"Compute 1-dimensional particle-in-a-box energy.\"\n",
    "    return n**2 * constants.h**2 / (8 * constants.m_e* a**2)\n",
    "\n",
    "#Define variable for atomic unit of length in meters\n",
    "a0 = constants.value('atomic unit of length')\n",
    "\n",
    "\n",
    "#This next bit of code just prints out the energy in atomic and SI units\n",
    "def print_energy(a, n):\n",
    "    print(f'The energy of an electron in a box of length {a:.2f} a.u. with '\n",
    "       f'quantum number {n} is {compute_energy(n, a):.2f} a.u..')\n",
    "    print(f'The energy of an electron in a box of length {a*a0:.2e} m with '\n",
    "       f'quantum number {n} is {compute_energy_si(n, a*a0):.2e} Joules.')\n",
    "\n",
    "out = widgets.interactive_output(print_energy, {'a': a, 'n': n})\n",
    "\n",
    "widgets.VBox([widgets.VBox([a, n]),out])  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Postulate: The wavefunction contains all the physically meaningful information about a system.\n",
    "While the wavefunction is not itself observable, all observable properties of a system can be determined from the wavefunction. However, just because the wavefunction encapsulates all the *observable* properties of a system does not mean that it contains *all information* about a system. In quantum mechanics, some things are not observable. Consider that for the ground ($n=1$) state of the particle in a box, the root-mean-square average momentum,\n",
    "$$\n",
    "\\bar{p}_{rms} = \\sqrt{2m \\cdot T} = \\sqrt{(2m)\\frac{h^2n^2}{8ma^2}} = \\frac{hn}{2a}\n",
    "$$\n",
    "increases as you squeeze the box. That is, the more you try to constrain the particle in space, the faster it moves. You can't \"stop\" the particle no matter how hard you squeeze it, so it's impossible to exactly know where the particle is located. You can only determine its *average* position."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Postulate: Observable Quantities Correspond to Linear, Hermitian Operators.\n",
    "The *correspondence* principle says that for every classical observable there is a linear, Hermitian, operator that allows computation of the quantum-mechanical observable. An operator, $\\hat{C}$ is linear if for any complex numbers $a$ and $b$, and any wavefunctions $\\psi_1(x)$ and $\\psi_2(x)$, \n",
    "$$\n",
    "\\hat{C} \\left(a \\psi_1(x,t) + b \\psi_2(x,t) \\right) = a \\hat{C} \\psi_1(x,t) + b \\hat{C} \\psi_2(x,t)\n",
    "$$\n",
    "Similarly, an operator is Hermitian if it satisfies the relation,\n",
    "$$\n",
    "\\int \\psi_1^*(x,t) \\hat{C} \\psi_2(x,t) \\, dx = \\int \\left( \\hat{C}\\psi_1(x,t) \\right)^*  \\psi_2(x,t) \\, dx \n",
    "$$\n",
    "or, equivalently,\n",
    "$$\n",
    "\\int \\psi_1^*(x,t) \\left( \\hat{C} \\psi_2(x,t) \\right) \\, dx = \\int \\psi_2(x,t)\\left( \\hat{C} \\psi_1(x,t) \\right)^* \\, dx \n",
    "$$\n",
    "That is for a linear operator, the linear operator applied to a sum of wavefunctions is equal to the sum of the linear operators directly applied to the wavefunctions separately, and the linear operator applied to a constant times a wavefunction is the constant times the linear operator applied directly to the wavefunction. A Hermitian operator can apply forward (towards $\\psi_2(x,t)$) or backwards (towards $\\psi_1(x,t)$). This is very useful, because sometimes it is much easier to apply an operator in one direction.\n",
    "\n",
    "We've already been introduced to the quantum-mechanical operators for the momentum, \n",
    "$$\n",
    "\\hat{p} = -i \\hbar \\tfrac{d}{dx}\n",
    "$$\n",
    "and the kinetic energy,\n",
    "$$\n",
    "\\hat{T} = -\\tfrac{\\hbar^2}{2m} \\tfrac{d^2}{dx^2}\n",
    "$$\n",
    "These operators are linear because the derivative of a sum is the sum of the derivatives, and the derivative of a constant times a function is that constant times the derivative of the function. These operators are also Hermitian. For example, to show that the momentum operator is Hermitian:\n",
    "\n",
    "\\begin{align}\n",
    "\\int_{-\\infty}^{\\infty} \\psi_1^*(x,t) \\hat{p} \\psi_2(x,t) dx &= \\int_{-\\infty}^{\\infty} \\psi_1^*(x,t) \\left( -i \\hbar \\tfrac{d}{dx} \\right) \\psi_2(x,t) dx \\\\\n",
    "&= -i \\hbar \\int_{-\\infty}^{\\infty} \\tfrac{d}{dx} \\left(\\psi_1^*(x,t)\\psi_2(x,t)  \\right) - \\left( \n",
    " \\psi_2(x,t) \\tfrac{d}{dx} \\psi_1^*(x,t)\\right)  dx \n",
    "\\end{align}\n",
    "Here we used the product rule for derivatives, $f(x)\\tfrac{dg}{dx} = \\tfrac{d f(x) g(x)}{dx} - g(x) \\tfrac{df}{dx}$. Using the fundamental theorem of calculus and the fact the probability of observing a particle at $\\pm \\infty$ is zero, and therefore the wavefunctions at infinity are also zero, one knows that\n",
    "$$\n",
    "\\int_{-\\infty}^{\\infty} \\tfrac{d}{dx} \\left(\\psi_1^*(x,t)\\psi_2(x,t)  \\right) = \\left[  \\psi_1^*(x,t)\\psi_2(x,t)\\right]_{-\\infty}^{\\infty} = 0\n",
    "$$\n",
    "Therefore the above equation can be simplified to\n",
    "\\begin{align}\n",
    "\\int_{-\\infty}^{\\infty} \\psi_1^*(x,t) \\hat{p} \\psi_2(x,t) dx &=i \\hbar \\int_{-\\infty}^{\\infty}   \\psi_2(x,t) \\tfrac{d}{dx} \\psi_1^*(x,t)  dx \\\\\n",
    "&= \\int_{-\\infty}^{\\infty}   \\psi_2(x,t) i \\hbar  \\tfrac{d}{dx} \\psi_1^*(x,t) dx \\\\\n",
    "&= \\int_{-\\infty}^{\\infty}  \\psi_2(x,t) \\left( -i \\hbar  \\tfrac{d}{dx} \\psi_1^*(x,t)\\right)  dx \\\\\n",
    "&= \\int_{-\\infty}^{\\infty}  \\psi_2(x,t) \\left( \\hat{p} \\psi_1(x,t)\\right)^*  dx \n",
    "\\end{align}\n",
    "\n",
    "The expectation value of the momentum of a particle-in-a-box is always zero. This is intuitive, since electrons (on average) are neither moving to the right nor to the left inside the box: if they were, then the box itself would need to be moving. Indeed, for any real wavefunction, the average momentum is always zero. This follows directly from the previous derivation with $\\psi_1^*(x,t) = \\psi_2(x,t)$. Thus: \n",
    "\\begin{align}\n",
    "\\int_{-\\infty}^{\\infty} \\psi_2(x,t) \\hat{p} \\psi_2(x,t) dx \n",
    "&=i \\hbar \\int_{-\\infty}^{\\infty}   \\psi_2(x,t) \\tfrac{d}{dx} \\psi_2(x,t)  dx \\\\\n",
    "&=-i \\hbar \\int_{-\\infty}^{\\infty}  \\psi_2(x,t) \\left( \\tfrac{d}{dx} \\psi_2(x,t)\\right)  dx \\\\\n",
    "&= 0\n",
    "\\end{align}\n",
    "The last line follows because the only number that is equal to its negative is zero. (That is, $x=-x$ if and only if $x=0$.) It is a subtle feature that the eigenfunctions of a real-valued Hamiltonian operator can always be chosen to be real-valued themselves, so their average momentum is clearly zero. We often denote quantum-mechanical expectation values with the shorthand,\n",
    "$$\n",
    "\\langle \\hat{p} \\rangle =0\n",
    "$$\n",
    "\n",
    "The momentum-squared of the particle-in-a-box is easily computed from the kinetic energy, \n",
    "$$\n",
    "\\langle \\hat{p}^2 \\rangle = \\int_0^a \\psi_n(x) \\hat{p}^2 \\psi_n(x) dx = \\int_0^a \\psi_n(x) \\left(2m\\hat{T}\\right) \\psi_n(x) dx = 2m E_n = \\frac{h^2n^2}{4a^2}\n",
    "$$\n",
    "\n",
    "Intuitively, since the box is symmetric about $x=\\tfrac{a}{2}$, the particle has equal probability of being in the first-half and the second-half of the box. So the average position is expected to be \n",
    "$$ \\langle x \\rangle =\\tfrac{a}{2}$$\n",
    "\n",
    "We can confirm this by explicit integration,\n",
    "\n",
    "\\begin{align}\n",
    "\\langle x \\rangle &= \\int_0^a \\psi_n^*(x)\\, x \\,\\psi_n(x) dx \\\\\n",
    "&= \\int_0^a \\left(\\sqrt{\\tfrac{2}{a}} \\sin\\left(\\tfrac{n \\pi x}{a} \\right)\\right) x \\left(\\sqrt{\\tfrac{2}{a}}\\sin\\left(\\tfrac{n \\pi x}{a} \\right)\\right) dx \\\\\n",
    "&= \\tfrac{2}{a} \\int_0^a  x \\sin^2\\left(\\tfrac{n \\pi x}{a} \\right) dx \\\\\n",
    "&= \\tfrac{2}{a} \\left[ \\tfrac{x^2}{4} - x \\tfrac{\\sin \\tfrac{2n \\pi x}{a}}{\\tfrac{4 n \\pi}{a}}  \n",
    "- \\tfrac{\\cos \\tfrac{2n \\pi x}{a}}{\\tfrac{8 n^2 \\pi^2}{a^2}}  \n",
    " \\right]_0^a \\\\\n",
    " &= \\tfrac{2}{a} \\left[ \\tfrac{a^2}{4} - 0 - 0 \\right] \\\\\n",
    " &= \\tfrac{a}{2}\n",
    "\\end{align} \n",
    "\n",
    "Similarly, we expect that the expectation value of $\\langle x^2 \\rangle$ will be proportional to $a^2$. We can confirm this by explicit integration, \n",
    "\\begin{align}\n",
    "\\langle x^2 \\rangle &= \\int_0^a \\psi_n^*(x)\\, x^2 \\,\\psi_n(x) dx \\\\\n",
    "&= \\int_0^a \\left(\\sqrt{\\tfrac{2}{a}} \\sin\\left(\\tfrac{n \\pi x}{a} \\right)\\right) x^2 \\left(\\sqrt{\\tfrac{2}{a}}\\sin\\left(\\tfrac{n \\pi x}{a} \\right)\\right) dx \\\\\n",
    "&= \\tfrac{2}{a} \\int_0^a  x^2 \\sin^2\\left(\\tfrac{n \\pi x}{a} \\right) dx \\\\\n",
    "&= \\tfrac{2}{a} \\left[ \\tfrac{x^3}{6} \n",
    " - x^2 \\tfrac{\\sin \\tfrac{2n \\pi x}{a}}{\\tfrac{4 n \\pi}{a}}\n",
    " - x \\tfrac{\\cos \\tfrac{2n \\pi x}{a}}{\\tfrac{4 n^2 \\pi^2}{a^2}}\n",
    " - \\tfrac{\\sin \\tfrac{2n \\pi x}{a}}{\\tfrac{8 n^3 \\pi^3}{a^3}}\n",
    " \\right]_0^a \\\\\n",
    " &= \\tfrac{2}{a} \\left[ \\tfrac{a^3}{6} - 0 - \\tfrac{a}{{\\tfrac{4 n^2 \\pi^2}{a^2}}} - 0 \\right] \\\\\n",
    "  &= \\tfrac{2}{a} \\left[ \\tfrac{a^3}{6} - \\tfrac{a^3}{4 n^2 \\pi^2} \\right] \\\\\n",
    "  &= a^2\\left[ \\tfrac{1}{3} - \\tfrac{1}{2 n^2 \\pi^2} \\right] \n",
    "\\end{align} \n",
    "\n",
    "We can verify these formulas by explicit integration. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3e34c291de69453ca0c4f2fff7206d3d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(VBox(children=(FloatSlider(value=1.0, description='a (length):', max=10.0, min=0.01, step=0.01)…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#The next few lines just set up the sliders for setting parameters.\n",
    "#Principle quantum number slider:\n",
    "n = widgets.IntSlider(\n",
    "    value=1,\n",
    "    min=1,\n",
    "    max=10,\n",
    "    step=1,\n",
    "    description='n (quant. #):',\n",
    "    disabled=False,\n",
    "    continuous_update=True,\n",
    "    orientation='horizontal',\n",
    "    readout=True,\n",
    "    readout_format='d')\n",
    "\n",
    "#Box length slider:\n",
    "a = widgets.FloatSlider(\n",
    "    value=1,\n",
    "    min=.01,\n",
    "    max=10.0,\n",
    "    step=0.01,\n",
    "    description='a (length):',\n",
    "    disabled=False,\n",
    "    continuous_update=True,\n",
    "    orientation='horizontal',\n",
    "    readout=True,\n",
    "    readout_format='.2f',\n",
    ")\n",
    "\n",
    "#Compute <x^power>, the expectation value of x^power\n",
    "def compute_moment(x, n, a, power):\n",
    "    \"\"\"Compute the x^power moment of the 1-dimensional particle-in-a-box.\n",
    "    \n",
    "    See `compute_wavefunction` parameters.\n",
    "    \"\"\"\n",
    "    return compute_probability(x, n, a)*x**power\n",
    "\n",
    "#This next bit of code just prints out the values.\n",
    "def check_moments(a, n):\n",
    "    #check the computed values of the moments against the analytic formula\n",
    "    avg_r,error = quad(compute_moment, 0, a, args=(n, a, 1))\n",
    "    avg_r2,error = quad(compute_moment, 0, a, args=(n, a, 2))\n",
    "    print(f\"<r> computed = {avg_r:.5f}\")\n",
    "    print(f\"<r> analytic = {a/2:.5f}\")\n",
    "    print(f\"<r^2> computed = {avg_r2:.5f}\")\n",
    "    print(f\"<r^2> analytic = {a**2*(1/3 - 1./(2*n**2*np.pi**2)):.5f}\")\n",
    "\n",
    "out = widgets.interactive_output(check_moments, {'a': a, 'n': n})\n",
    "\n",
    "widgets.HBox([widgets.VBox([a, n]),out])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Heisenberg Uncertainty Principle\n",
    "The previous example gives a first example of the more general [Heisenberg Uncertainty Principle](https://en.wikipedia.org/wiki/Uncertainty_principle). One specific manifestation of the Heisenberg Uncertaity Principle is that the variance of the position, $\\sigma_x^2 = \\langle x^2 \\rangle - \\langle x \\rangle^2$, times the variance of the momentum. $\\sigma_p^2 = \\langle \\hat{p}^2 \\rangle - \\langle \\hat{p} \\rangle^2$ is greater than $\\tfrac{\\hbar^2}{4}$. We can verify this formula for the particle in a box. \n",
    "\\begin{align}\n",
    "\\tfrac{\\hbar^2}{4} &\\le \\sigma_x^2 \\sigma_p^2 \\\\\n",
    "&= \\left( a^2\\left[ \\tfrac{1}{3} - \\tfrac{1}{2 n^2 \\pi^2} \\right] - \\tfrac{a^2}{4} \\right)\n",
    " \\left( \\frac{h^2n^2}{4a^2} - 0 \\right) \\\\\n",
    "&= \\left( a^2\\left[ \\tfrac{1}{3} - \\tfrac{1}{2 n^2 \\pi^2} \\right] - \\tfrac{a^2}{4} \\right)\n",
    " \\left( \\frac{\\hbar^2 \\pi^2 n^2}{a^2} \\right) \\\\\n",
    "&= \\hbar^2 \\pi^2 n^2 \\left(\\tfrac{1}{12} - \\tfrac{1}{2 n^2 \\pi^2} \\right) \n",
    "\\end{align}\n",
    "The right-hand-side gets larger and larger as $n$ increases, so the largest value occurs where:\n",
    "\\begin{align}\n",
    "\\tfrac{\\hbar^2}{4} &\\le \\hbar^2 \\pi^2 \\left(\\tfrac{1}{12} - \\tfrac{1}{2 \\pi^2} \\right) \\\\\n",
    "&= 0.32247 \\hbar^2\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Double-Checking the Energy of a Particle-in-a-Box\n",
    "To check the energy of the particle in a box, we can compute the kinetic energy density, then integrate it over all space. That is, we define:\n",
    "$$\n",
    "\\tau_n(x) = \\psi_n^*(x) \\left(-\\tfrac{\\hbar^2}{2m}\\tfrac{d^2}{dx^2}\\right) \\psi_n(x)\n",
    "$$\n",
    "and then the kinetic energy (which is the energy for the particle in a box) is\n",
    "$$\n",
    "T_n = \\int_0^a \\tau_n(x) dx\n",
    "$$\n",
    "\n",
    "> *Note:* In fact there are several different equivalent definitions of the kinetic energy density, but this is not very important in an introductory quantum chemistry course. All of the kinetic energy densities give the same total kinetic energy. However, because the kinetic energy density, $\\tau(x)$, represents the kinetic energy of a particle at the point $x$, and it is impossible to know the momentum (or the momentum-squared, ergo the kinetic energy) exactly at a given point in space according to the Heisenberg uncertainty principle), there can be no unique definition for $\\tau(x)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e69800832d224d97a05b816119b1e150",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(VBox(children=(FloatSlider(value=1.0, description='a (length):', max=10.0, min=0.01, step=0.01)…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#The next few lines just set up the sliders for setting parameters.\n",
    "#Principle quantum number slider:\n",
    "n = widgets.IntSlider(\n",
    "    value=1,\n",
    "    min=1,\n",
    "    max=10,\n",
    "    step=1,\n",
    "    description='n (quant. #):',\n",
    "    disabled=False,\n",
    "    continuous_update=True,\n",
    "    orientation='horizontal',\n",
    "    readout=True,\n",
    "    readout_format='d')\n",
    "\n",
    "#Box length slider:\n",
    "a = widgets.FloatSlider(\n",
    "    value=1,\n",
    "    min=.01,\n",
    "    max=10.0,\n",
    "    step=0.01,\n",
    "    description='a (length):',\n",
    "    disabled=False,\n",
    "    continuous_update=True,\n",
    "    orientation='horizontal',\n",
    "    readout=True,\n",
    "    readout_format='.2f',\n",
    ")\n",
    "\n",
    "#The next few lines just set up the sliders for setting parameters.\n",
    "#Principle quantum number slider:\n",
    "\n",
    "def compute_wavefunction_derivative(x, n, a, order=1):\n",
    "    \"\"\"Compute 1-dimensional particle-in-a-box kinetic energy density.\n",
    "    \"\"\"\n",
    "    if not (isinstance(order, int) and n > 0):\n",
    "        raise ValueError(\"Argument order is expected to be a positive integer!\")\n",
    "\n",
    "    def wavefunction(x):\n",
    "        v = np.sqrt(2 / a) * np.sin(n * np.pi * x / a)\n",
    "        return v\n",
    "\n",
    "    # compute derivative\n",
    "    deriv = egrad(wavefunction)\n",
    "    for _ in range(order - 1):\n",
    "        deriv = egrad(deriv)\n",
    "\n",
    "    # return zero for x values out of the box\n",
    "    deriv = deriv(x)\n",
    "#     deriv[x < 0] = 0.0\n",
    "#     deriv[x > a] = 0.0\n",
    "\n",
    "    if hasattr(x, \"__iter__\"):\n",
    "        deriv[x > a] = 0.0\n",
    "        deriv[x < 0] = 0.0\n",
    "    else:\n",
    "        if x < 0.0 or x > a:\n",
    "            deriv = 0.0\n",
    "    \n",
    "    return deriv\n",
    "\n",
    "\n",
    "def compute_kinetic_energy_density(x, n, a):\n",
    "    \"\"\"Compute 1-dimensional particle-in-a-box kinetic energy density.\n",
    "    \n",
    "    See `compute_wavefunction` parameters.\n",
    "    \"\"\"\n",
    "    # evaluate wave-function and its 2nd derivative w.r.t. x\n",
    "    wf = compute_wavefunction(x, n, a)\n",
    "    d2 = compute_wavefunction_derivative(x, n, a, order=2)\n",
    "    return -0.5 * wf * d2\n",
    "\n",
    "#This next bit of code just prints out the values.\n",
    "def check_energy(a, n):\n",
    "    #check the computed values of the moments against the analytic formula\n",
    "    ke,error = quad(compute_kinetic_energy_density, 0, a, args=(n, a))\n",
    "    energy = compute_energy(n, a)\n",
    "    print(f\"The energy computed by integrating the k.e. density is {ke:.5f}\")\n",
    "    print(f\"The energy computed directly is {energy:.5f}\")\n",
    "\n",
    "out = widgets.interactive_output(check_energy, {'a': a, 'n': n})\n",
    "\n",
    "widgets.VBox([widgets.VBox([a, n]),out])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing the Particle-in-a-Box Wavefunctions, Probabilities, etc.\n",
    "In the next code block, the wavefunction, probability density, derivative, second-derivative, and kinetic energy density for the particle-in-a-box are shown. Notice that the kinetic-energy-density is proportional to the probability density, and that the first and second derivatives are not zero at the edge of the box, but the wavefunction and probability density are. It's useful to change the parameters in the below figures to build your intuition for the particle-in-a-box."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0ddf05e29ca74f189e6c46d6a2377de0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(VBox(children=(FloatSlider(value=1.0, description='a (length):', max=10.0, min=0.01, step=0.01)…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#The next few lines just set up the sliders for setting parameters.\n",
    "#Principle quantum number slider:\n",
    "n = widgets.IntSlider(\n",
    "    value=1,\n",
    "    min=1,\n",
    "    max=10,\n",
    "    step=1,\n",
    "    description='n (quant. #):',\n",
    "    disabled=False,\n",
    "    continuous_update=True,\n",
    "    orientation='horizontal',\n",
    "    readout=True,\n",
    "    readout_format='d')\n",
    "\n",
    "#Box length slider:\n",
    "a = widgets.FloatSlider(\n",
    "    value=1,\n",
    "    min=.01,\n",
    "    max=10.0,\n",
    "    step=0.01,\n",
    "    description='a (length):',\n",
    "    disabled=False,\n",
    "    continuous_update=True,\n",
    "    orientation='horizontal',\n",
    "    readout=True,\n",
    "    readout_format='.2f',\n",
    ")\n",
    "\n",
    "#This next bit of code makes the plots and prints out the energy\n",
    "def make_plots(a, n):\n",
    "    #check the computed values of the moments against the analytic formula\n",
    "    energy = compute_energy(n, a)\n",
    "    print(f\"The energy computed directly is {energy:.5f}\")\n",
    "    # sample x coordinates\n",
    "    x = np.arange(-0.6, a + 0.6, 0.01)\n",
    "\n",
    "    # evaluate wave-function & probability\n",
    "    wf = compute_wavefunction(x, n, a)\n",
    "    pr = compute_probability(x, n, a)\n",
    "\n",
    "    # evaluate 1st & 2nd derivative of wavefunction w.r.t. x\n",
    "    d1 = compute_wavefunction_derivative(x, n, a, order=1)\n",
    "    d2 = compute_wavefunction_derivative(x, n, a, order=2)\n",
    "\n",
    "    # evaluate kinetic energy density\n",
    "    kin = compute_kinetic_energy_density(x, n, a)\n",
    "    #print(\"Integrate KED = \", trapezoid(kin, x))\n",
    "\n",
    "\n",
    "    # set the size of the plot\n",
    "    plt.rcParams['figure.figsize'] = [15, 10]\n",
    "    plt.rcParams['font.family'] = 'DejaVu Sans'\n",
    "    plt.rcParams['font.serif'] = ['Times New Roman']\n",
    "    plt.rcParams['mathtext.fontset'] = 'stix'\n",
    "    plt.rcParams['xtick.labelsize'] = 15\n",
    "    plt.rcParams['ytick.labelsize'] = 15\n",
    "\n",
    "    # define four subplots\n",
    "    fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2)\n",
    "    fig.suptitle(f'a={a}  n={n}  E={compute_energy(n, a):.4f} a.u.', fontsize=35, fontweight='bold')\n",
    "\n",
    "    # plot 1\n",
    "    ax1.plot(x, wf, linestyle='--', label=r'$\\psi(x)$', linewidth=3)\n",
    "    ax1.plot(x, pr, linestyle='-', label=r'$\\|\\psi(x)\\|^2$', linewidth=3)\n",
    "    ax1.legend(frameon=False, fontsize=20)\n",
    "    ax1.set_xlabel('x coordinates', fontsize=20)\n",
    "\n",
    "    # plot 2\n",
    "    ax2.plot(x, d1, linewidth=3, c='k')\n",
    "    ax2.set_xlabel('x coordinates', fontsize=20)\n",
    "    ax2.set_ylabel(r'$\\frac{d\\psi(x)}{dx}$', fontsize=30, rotation=0, labelpad=25)\n",
    "\n",
    "    # plot 3\n",
    "    ax3.plot(x, d2, linewidth=3, c='g', )\n",
    "    ax3.set_xlabel('x coordinates', fontsize=20)\n",
    "    ax3.set_ylabel(r'$\\frac{d^2\\psi(x)}{dx^2}$', fontsize=25, rotation=0, labelpad=25)\n",
    "\n",
    "    # plot 4\n",
    "    ax4.plot(x, kin, linewidth=3, c='r')\n",
    "    ax4.set_xlabel('x coordinates', fontsize=20)\n",
    "    ax4.set_ylabel('Kinetic Energy Density', fontsize=16)\n",
    "\n",
    "    # adjust spacing between plots\n",
    "    plt.subplots_adjust(left=0.125,\n",
    "                    bottom=0.1, \n",
    "                    right=0.9, \n",
    "                    top=0.9, \n",
    "                    wspace=0.35, \n",
    "                    hspace=0.35)\n",
    "    #Show Plot\n",
    "    plt.show()\n",
    "    \n",
    "    \n",
    "out = widgets.interactive_output(make_plots, {'a': a, 'n': n})\n",
    "\n",
    "widgets.VBox([widgets.VBox([a, n]),out])\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## &#x1fa9e; Self-Reflection\n",
    "- Can you think of other physical or chemical systems where the particle-in-a-box Hamiltonian would be appropriate?\n",
    "- Can you think of another property density, besides the kinetic energy density, that cannot be uniquely defined in quantum mechanics?\n",
    "\n",
    "## &#x1f914; Thought-Provoking Questions\n",
    "- How would the wavefunction and ground-state energy change if you made a small change in the particle-in a box Hamiltonian, so that the right-hand-side of the box was a little higher than the left-hand-side? \n",
    "- Any system with a zero-point energy of zero is classical. Why? \n",
    "- How would you compute the probability of observing an electron at the center of a box if the box contained 2 electrons? If it contained 4 electrons? If it contained 8 electrons? The probability of observing an electron at the center of a box with 3 electrons is sometimes *lower* than the probability of observing an electron at the center of a box with 2 electrons. Why? \n",
    "- Demonstrate that the kinetic energy operator is linear and Hermitian.\n",
    "\n",
    "## &#x1f501; Recapitulation\n",
    "- Write the Hamiltonian, time-independent Schr&ouml;dinger equation, eigenfunctions, and eigenvalues for the one-dimensional particle in a box.\n",
    "- Play around with the eigenfunctions and eigenenergies of the particle-in-a-box to build intuition for them.\n",
    "- How would you compute the uncertainty in $x^4$?\n",
    "- Practice your calculus by explicitly computing, using integration by parts, $\\langle x \\rangle$ and $\\langle x^2 \\rangle$.\n",
    "\n",
    "## &#x1f52e; Next Up...\n",
    "- Postulates of Quantum Mechanics\n",
    "- Multielectron particle-in-a-box.\n",
    "- Multidimensional particle-in-a-box.\n",
    "- Harmonic Oscillator\n",
    "\n",
    "## &#x1f4da; References\n",
    "My favorite sources for this material are:\n",
    "- [Randy's book](https://github.com/PaulWAyers/IntroQChem/blob/main/documents/DumontBook.pdf?raw=true) has an excellent treatment of the particle-in-a-box model, including several extensions to the material covered here. (Chapter 3)\n",
    "- Also see my (pdf) class [notes]\n",
    "(https://github.com/PaulWAyers/IntroQChem/blob/main/documents/PinBox.pdf?raw=true).\n",
    "- Also see my notes on the [mathematical structure of quantum mechanics](https://github.com/PaulWAyers/IntroQChem/blob/main/documents/LinAlgAnalogy.pdf?raw=true).\n",
    "- D. A. MacQuarrie, Quantum Chemistry (University Science Books, Mill Valley California, 1983)\n",
    "- [An excellent explanation of the link to the spectrum of cyanine dyes](https://pubs.acs.org/doi/10.1021/ed084p1840)\n",
    "- Chemistry Libre Text: [one dimensional](https://chem.libretexts.org/Bookshelves/Physical_and_Theoretical_Chemistry_Textbook_Maps/Map%3A_Physical_Chemistry_for_the_Biosciences_(Chang)/11%3A_Quantum_Mechanics_and_Atomic_Structure/11.08%3A_Particle_in_a_One-Dimensional_Box)\n",
    "and [multi-dimensional](https://chem.libretexts.org/Bookshelves/Physical_and_Theoretical_Chemistry_Textbook_Maps/Supplemental_Modules_(Physical_and_Theoretical_Chemistry)/Quantum_Mechanics/05.5%3A_Particle_in_Boxes/Particle_in_a_3-Dimensional_box)\n",
    "- [McQuarrie and Simon summary](https://chem.libretexts.org/Bookshelves/Physical_and_Theoretical_Chemistry_Textbook_Maps/Map%3A_Physical_Chemistry_(McQuarrie_and_Simon)/03%3A_The_Schrodinger_Equation_and_a_Particle_in_a_Box)\n",
    "- [Mathematical Features of Quantum Mechanics]() (my notes).\n",
    "\n",
    "There are also some excellent wikipedia articles:\n",
    "- [Particle in a Box](https://en.wikipedia.org/wiki/Particle_in_a_box)\n",
    "- [Particle on a Ring](https://en.wikipedia.org/wiki/Particle_in_a_ring)\n",
    "- [Postulates of Quantum Mechanics](https://en.wikipedia.org/wiki/Mathematical_formulation_of_quantum_mechanics#Postulates_of_quantum_mechanics)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  },
  "toc": {
   "base_numbering": "3",
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "343px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}