notebook/quantum-mechanics/2quantumwells.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# **Numerical Solution of the Schrödinger Equation for the Double Square Well Potential**\n",
    "\n",
    "**Authors:** Dou Du, Taylor James Baird and Giovanni Pizzi \n",
    "\n",
    "<i class=\"fa fa-home fa-2x\"></i><a href=\"../index.ipynb\" style=\"font-size: 20px\"> Go back to index</a>\n",
    "\n",
    "**Source code:** https://github.com/osscar-org/quantum-mechanics/blob/master/notebook/quantum-mechanics/2quantumwells.ipynb\n",
    "\n",
    "\n",
    "This notebook displays interactively the eigenfunctions for a double square well potential (DSWP) in one dimension, as obtained from the numerical solution of the associated Schrödinger equation. The double square well\n",
    "potential is a simple but effective model to describe physical systems.\n",
    "\n",
    "<hr style=\"height:1px;border:none;color:#cccccc;background-color:#cccccc;\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **Goals**\n",
    "* Understand one of the numerical methods used to solve the time-independent Schrödinger equation.\n",
    "* Understand the importance of the DSWP model for real material systems.\n",
    "* Investigate quantum tunneling in the double square well potential system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **Background Theory**\n",
    "\n",
    "[**More in the background theory.**](./theory/theory_2quantumwells.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **Tasks and exercises**\n",
    "\n",
    "1. Start from a small gap distance (e.g. 0.1) between the two quantum wells.\n",
    "Fix the depth and width for both wells to their default values. Keep both wells\n",
    "identical. Slowly increase (use 0.1 as step value) the gap distance and record \n",
    "the two lowest eigenvalues at each gap distance. What do you \n",
    "see when the gap distance is increasing? Can you think of a chemical reaction that could be modeled by this process?\n",
    "\n",
    "    <details>\n",
    "    <summary style=\"color: red\">Solution</summary>\n",
    "    <div style=\"border:blue; border-width:3px; border-style:outset;\">\n",
    "    When the gap distance is increased, the energy difference between \n",
    "    the two lowest eigenvalues becomes smaller. When the gap distance becomes significant, \n",
    "    the wavefunctions are localized in both quantum wells, which converge \n",
    "    to single quantum well solutions. This process models the breaking of an H$_2$ bond.\n",
    "    </div>\n",
    "    </details>\n",
    "\n",
    "2. In the accompanying background theory notebook, we employed symmetric double wells to describe the H$_2$ molecule and ammonia maser examples. This means that the width and depth are identical for both wells. However, it is possible to imagine a situation in which the two wells are different. What kinds of chemical systems can be described by two quantum wells with different depths and widths?\n",
    "\n",
    "    <details>\n",
    "    <summary style=\"color: red\">Solution</summary>\n",
    "    <div style=\"border:blue; border-width:3px; border-style:outset;\">\n",
    "    Generally, the double square well potential can be employed to characterize diatomic molecules, which only contain two atoms. There are two categories \n",
    "    of diatomic molecules, namely homonuclear and heteronuclear molecules. \n",
    "    Homonuclear diatomic molecules consist of two identical atoms, \n",
    "    such as H$_2$, O$_2$, and N$_2$. Heteronuclear diatomic molecules have \n",
    "    two different atoms like CO, NO, and HCl. The symmetric DSWP is a reasonable \n",
    "    model for homonuclear diatomic molecules since the electron-nuclear \n",
    "    interactions are the same for both nuclei. On the other hand, one needs \n",
    "    to set two different quantum wells to approximate heteronuclear \n",
    "    diatomic molecules.\n",
    "    </div>\n",
    "    </details>\n",
    "\n",
    "3. What is the disadvantage of this double square well potential model?\n",
    "How would one improve the model?\n",
    "\n",
    "    <details>\n",
    "    <summary style=\"color: red\">Solution</summary>\n",
    "    <div style=\"border:blue; border-width:3px; border-style:outset;\">\n",
    "    Even though the DSWP is successful in the explanation of some quantum \n",
    "    systems, the model is still too simple to be applied in many cases and is far from the real potential. A more realistic potential should not change value so sharply, but should instead move somewhat smoothly from a finite to zero value.\n",
    "    One can replace the square well with another, more realistic potential, \n",
    "    for example, a quartic function as shown in the \n",
    "    <a href=\"./asymmetricwell.ipynb\">\"Avoided Crossing \n",
    "    in 1D Asymmetric Quantum Well\"</a> notebook.\n",
    "    </div>\n",
    "    </details>\n",
    "    \n",
    "4. By making a comparison with the solutions expected from classical considerations, think of how quantum mechanical tunnelling manifests in this system. Furthermore, try to think of the role this phenomenon plays in the operation of an ammonia maser.\n",
    "\n",
    "   <details>\n",
    "    <summary style=\"color: red\">Solution</summary>\n",
    "    <div style=\"border:blue; border-width:3px; border-style:outset;\">  \n",
    "    In the classical picture of this system, one would never expect a particle with energy $E$, smaller than the height of the well, $V$, to be observed in the region between the two wells or to move from its original well to the other. However, we can observe in the figure that in the quantum mechanical treatment of the system, there is in fact a non-zero probability of finding a particle in the region between the two wells and for hopping between wells. This latter behaviour is referred to as tunnelling. This phenomenon is one of the core principles behind an ammonia (NH$_3$) maser which relies on the tunnelling of the nitrogen atom through the potential barrier produced by the plane of hydrogen atoms. See the <a href=\"./theory/theory_2quantumwells.ipynb\">background theory notebook</a> for an illustration of this.\n",
    "    </div>\n",
    "    </details>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"height:1px;border:none;color:#cccccc;background-color:#cccccc;\" />\n",
    "\n",
    "## Interactive visualization\n",
    "(be patient, it might take a few seconds to load)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib widget\n",
    "\n",
    "from numpy import linspace, sqrt, ones, arange, diag, argsort, zeros\n",
    "from scipy.linalg import eigh_tridiagonal\n",
    "import matplotlib.pyplot as plt\n",
    "from ipywidgets import FloatSlider, jslink, VBox, HBox, Button, Label, RadioButtons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors = ['#e41a1c','#377eb8','#4daf4a','#984ea3','#ff7f00','#ffff33','#a65628','#f781bf']\n",
    "ixx = 0\n",
    "\n",
    "def doublewell_potential(x, width1, depth1, width2, depth2, dist):\n",
    "    xa = zeros(len(x))\n",
    "    xb = zeros(len(x))\n",
    "    \n",
    "    for i in range(len(x)):\n",
    "        if x[i] > -dist/2.0 - width1 and x[i] < -dist/2.0:\n",
    "            xa[i] = depth1\n",
    "\n",
    "    for i in range(len(x)):\n",
    "        if x[i] > dist/2.0 and x[i] < dist/2.0 + width2:\n",
    "            xb[i] = depth2\n",
    "    \n",
    "    return xa + xb\n",
    "            \n",
    "    \n",
    "def diagonalization(hquer, L, N, pot = doublewell_potential, width1 = 0.5, depth1 = -0.2,\n",
    "                     width2 = 0.5, depth2 = -0.2, dist = 1.0):\n",
    "    \"\"\"Calculated sorted eigenvalues and eigenfunctions. \n",
    "\n",
    "       Input:\n",
    "         hquer: Planck constant\n",
    "         L: set viewed interval [-L,L] \n",
    "         N: number of grid points i.e. size of the matrix \n",
    "         pot: potential function of the form pot\n",
    "         x0: center of the quantum well\n",
    "         width: the width of the quantum well\n",
    "         depth: the depth of the quantum well\n",
    "       Ouput:\n",
    "         ew: sorted eigenvalues (array of length N)\n",
    "         ef: sorted eigenfunctions, ef[:,i] (size N*N)\n",
    "         x:  grid points (arry of length N)\n",
    "         dx: grid space\n",
    "         V:  Potential at positions x (array of length N)\n",
    "    \"\"\"\n",
    "    x = linspace(-L, L, N+2)[1:N+1]               # grid points \n",
    "    dx = x[1] - x[0]                              # grid spacing\n",
    "    V = pot(x, width1, depth1, width2, depth2, dist)\n",
    "    z = hquer**2 /2.0/dx**2                       # second diagonals\n",
    "\n",
    "    ew, ef = eigh_tridiagonal(V+2.0*z, -z*ones(N-1))    \n",
    "    ew = ew.real                                  # real part of the eigenvalues\n",
    "    ind = argsort(ew)                             # Indizes f. sort. Array\n",
    "    ew = ew[ind]                                  # Sort the ew by ind\n",
    "    ef = ef[:, ind]                               # Sort the columns\n",
    "\n",
    "    ef = ef/sqrt(dx)                              # Correct standardization \n",
    "    return ew, ef, x, dx, V\n",
    "\n",
    "\n",
    "def plot_eigenfunctions(ax, ew, ef, x, V, width=1, Emax=0.05, fak=2.0, single = 0):\n",
    "    \"\"\"Plot of the lowest eigenfunctions 'ef' in the potential 'V (x)'\n",
    "       at the level of the eigenvalues 'ew' in the plot area 'ax'.\n",
    "    \"\"\"\n",
    "    \n",
    "    if psi_x.value == \"Wavefunction\":\n",
    "        y_text = r'$V(x)$ and $\\psi(x)$'\n",
    "        fig.suptitle('Numerical Solution ($\\psi$) of One Dimension Schrödinger Equation', fontsize = 13)\n",
    "    else:\n",
    "        y_text = r'$V(x)$ and $|\\psi(x)|^2$'\n",
    "        fig.suptitle('Numerical Solution ($\\psi^2$) of One Dimension Schrödinger Equation', fontsize = 13)\n",
    "    \n",
    "    fak = fak/100.0;\n",
    "    \n",
    "    ax[0].axhspan(0.0, Emax, facecolor='lightgrey')\n",
    "    \n",
    "    ax[0].set_xlim([min(x), max(x)])\n",
    "    ax[0].set_ylim([min(V)-0.05, Emax])\n",
    "    \n",
    "    ax[0].set_xlabel(r'$x/a$', fontsize = 10)\n",
    "    ax[0].set_ylabel(y_text, fontsize = 10)\n",
    "    \n",
    "    ax[1].set_xlim([min(x), max(x)])\n",
    "    ax[1].set_ylim([min(V)-0.05, Emax])\n",
    "    \n",
    "    ax[1].yaxis.set_label_position(\"right\")\n",
    "    ax[1].yaxis.tick_right()\n",
    "    \n",
    "    ax[1].get_xaxis().set_visible(False)\n",
    "    ax[1].set_ylabel(r'$\\rm{\\ Eigenvalues}$', fontsize = 10)\n",
    "    \n",
    "    indmax = sum(ew<=0.0)                         \n",
    "    if not hasattr(width, \"__iter__\"):           \n",
    "        width = width*ones(indmax)               \n",
    "    for i in arange(indmax):\n",
    "        if psi_x.value == \"Wavefunction\":\n",
    "            ax[0].plot(x, fak*(ef[:, i])+ew[i], linewidth=width[i]+.1, color=colors[i%len(colors)])\n",
    "        else:\n",
    "            ax[0].plot(x, fak*abs(ef[:, i])**2+ew[i], linewidth=width[i]+.1, color=colors[i%len(colors)])\n",
    "\n",
    "        ax[1].plot(x, x*0.0+ew[i], linewidth=width[i]+2.5, color=colors[i%len(colors)])\n",
    "        \n",
    "    ax[0].plot(x, V, c='k', linewidth=1.6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L = 1.5                                              # x range [-L,L]\n",
    "N = 200                                              # Number of grid points \n",
    "hquer = 0.06                                         # Planck constant\n",
    "sigma_x = 0.1                                        # Width of the Guassian function\n",
    "zeiten = linspace(0.0, 10.0, 400)                    # Time\n",
    "\n",
    "style = {'description_width': 'initial'}\n",
    "\n",
    "swidth1 = FloatSlider(value = 0.5, min = 0.1, max = 1.0, description = 'Width (left): ', style = style)\n",
    "sdepth1 = FloatSlider(value = -0.2, min = -1.0, max = -0.1, description = 'Depth (left): ', style = style)\n",
    "\n",
    "swidth2 = FloatSlider(value = 0.5, min = 0.1, max = 1.0, description = 'Width (right): ', style = style)\n",
    "sdepth2 = FloatSlider(value = -0.2, min = -1.0, max = -0.1, description = 'Depth (right): ', style = style)\n",
    "\n",
    "sdist = FloatSlider(value = 0.5, min = 0.0, max = L, description = r'Gap distance: ', style = style)\n",
    "sfak = FloatSlider(value = 2, min = 1.0, max = 5.0, step = 0.5, description = r'Zoom factor: ', style = style)\n",
    "\n",
    "psi_x = RadioButtons(options=[\"Wavefunction\", \"Probability density\"], value=\"Wavefunction\", description=\"Show:\")\n",
    "\n",
    "ew, ef, x, dx, V = diagonalization(hquer, L, N)\n",
    "    \n",
    "fig, ax = plt.subplots(1, 2, figsize=(7,5), gridspec_kw={'width_ratios': [10, 1]})\n",
    "fig.canvas.header_visible = False\n",
    "fig.canvas.layout.width = \"750px\"\n",
    "\n",
    "fig.canvas.header_visible = False\n",
    "plot_eigenfunctions(ax, ew, ef, x, V)\n",
    "\n",
    "plt.show()\n",
    "\n",
    "def on_width_change1(change):\n",
    "    global ew, ef, x, dx, V\n",
    "    \n",
    "    try:\n",
    "        ann.remove()\n",
    "        ann1.remove()\n",
    "    except:\n",
    "        pass\n",
    "    \n",
    "    ax[0].lines.clear()\n",
    "    ax[1].lines.clear()\n",
    "    ew, ef, x, dx, V = diagonalization(hquer, L, N, \n",
    "                                        width1 = swidth1.value, depth1 = sdepth1.value,\n",
    "                                        width2 = swidth2.value, depth2 = sdepth2.value,\n",
    "                                        dist = sdist.value)\n",
    "    plot_eigenfunctions(ax, ew, ef, x, V)\n",
    "\n",
    "def on_depth_change1(change):\n",
    "    global ew, ef, x, dx, V\n",
    "       \n",
    "    try:\n",
    "        ann.remove()\n",
    "        ann1.remove()\n",
    "    except:\n",
    "        pass\n",
    "    \n",
    "    ax[0].clear()\n",
    "    ax[1].clear()\n",
    "    ew, ef, x, dx, V = diagonalization(hquer, L, N, \n",
    "                                        width1 = swidth1.value, depth1 = sdepth1.value,\n",
    "                                        width2 = swidth2.value, depth2 = sdepth2.value,\n",
    "                                        dist = sdist.value)\n",
    "    plot_eigenfunctions(ax, ew, ef, x, V)\n",
    "    \n",
    "\n",
    "def on_width_change2(change):\n",
    "    global ew, ef, x, dx, V\n",
    "    \n",
    "    try:\n",
    "        ann.remove()\n",
    "        ann1.remove()\n",
    "    except:\n",
    "        pass\n",
    "    \n",
    "    ax[0].lines.clear()\n",
    "    ax[1].lines.clear()\n",
    "    ew, ef, x, dx, V = diagonalization(hquer, L, N, \n",
    "                                        width1 = swidth1.value, depth1 = sdepth1.value,\n",
    "                                        width2 = swidth2.value, depth2 = sdepth2.value,\n",
    "                                        dist = sdist.value)\n",
    "    plot_eigenfunctions(ax, ew, ef, x, V)\n",
    "    \n",
    "def on_depth_change2(change):\n",
    "    global ew, ef, x, dx, V\n",
    "    \n",
    "    try:\n",
    "        ann.remove()\n",
    "        ann1.remove()\n",
    "    except:\n",
    "        pass\n",
    "    \n",
    "    ax[0].lines.clear()\n",
    "    ax[1].lines.clear()\n",
    "    ew, ef, x, dx, V = diagonalization(hquer, L, N, \n",
    "                                        width1 = swidth1.value, depth1 = sdepth1.value,\n",
    "                                        width2 = swidth2.value, depth2 = sdepth2.value,\n",
    "                                        dist = sdist.value)\n",
    "    plot_eigenfunctions(ax, ew, ef, x, V)\n",
    "    \n",
    "def on_dist_change(change):\n",
    "    global ew, ef, x, dx, V\n",
    "    \n",
    "    try:\n",
    "        ann.remove()\n",
    "        ann1.remove()\n",
    "    except:\n",
    "        pass\n",
    "    \n",
    "    ax[0].lines.clear()\n",
    "    ax[1].lines.clear()\n",
    "    ew, ef, x, dx, V = diagonalization(hquer, L, N, \n",
    "                                        width1 = swidth1.value, depth1 = sdepth1.value,\n",
    "                                        width2 = swidth2.value, depth2 = sdepth2.value,\n",
    "                                        dist = sdist.value)\n",
    "    plot_eigenfunctions(ax, ew, ef, x, V)\n",
    "    \n",
    "def on_xfak_change(change):\n",
    "    \n",
    "    try:\n",
    "        ann.remove()\n",
    "        ann1.remove()\n",
    "    except:\n",
    "        pass\n",
    "    \n",
    "    ax[0].lines.clear()\n",
    "    ax[1].lines.clear()\n",
    "\n",
    "    plot_eigenfunctions(ax, ew, ef, x, V, fak = sfak.value, single = ixx)\n",
    "\n",
    "def on_press(event):\n",
    "    global ann, ann1, ixx\n",
    "    \n",
    "    ixx = min(enumerate(ew), key = lambda x: abs(x[1]-event.ydata))[0]\n",
    "    \n",
    "    for i in range(len(ax[1].lines)):\n",
    "        ax[0].lines[i].set_alpha(0.1)\n",
    "        ax[1].lines[i].set_alpha(0.1)\n",
    "        ax[0].lines[i].set_linewidth(1.1)\n",
    "        \n",
    "    ax[0].lines[ixx].set_alpha(1.0)\n",
    "    ax[1].lines[ixx].set_alpha(1.0)\n",
    "    ax[0].lines[ixx].set_linewidth(2.0)\n",
    "    \n",
    "    try:\n",
    "        ann.remove()\n",
    "        ann1.remove()\n",
    "    except:\n",
    "        pass\n",
    "    \n",
    "    ann = ax[0].annotate(s = 'n = ' + str(ixx+1), xy = (0, ew[ixx]), xytext = (-0.15, ew[ixx]), xycoords = 'data', color='k', size=15)\n",
    "    ann1 = ax[1].annotate(s = str(\"{:.3f}\".format(ew[ixx])), xy = (0, ew[ixx]), xytext = (-1.2, ew[ixx]+0.005), xycoords = 'data', color='k', size=9)\n",
    "\n",
    "cid = fig.canvas.mpl_connect('button_press_event', on_press)\n",
    "\n",
    "def on_update_click(b):\n",
    "    for i in ax[0].lines:\n",
    "        i.set_alpha(1.0)\n",
    "    for i in ax[1].lines:\n",
    "        i.set_alpha(1.0)\n",
    "\n",
    "    try:\n",
    "        ann.remove()\n",
    "        ann1.remove()\n",
    "    except:\n",
    "        pass\n",
    "    \n",
    "update = Button(description=\"Show all\")\n",
    "update.on_click(on_update_click)\n",
    "\n",
    "swidth1.observe(on_width_change1, names = 'value')\n",
    "sdepth1.observe(on_depth_change1, names = 'value')\n",
    "swidth2.observe(on_width_change2, names = 'value')\n",
    "sdepth2.observe(on_depth_change2, names = 'value')\n",
    "sdist.observe(on_dist_change, names = 'value')\n",
    "sfak.observe(on_xfak_change, names = 'value')\n",
    "psi_x.observe(on_width_change1, names = 'value')\n",
    "\n",
    "label1 = Label(value=\"(click on a state to select it)\");\n",
    "\n",
    "display(HBox([swidth1, sdepth1]), HBox([swidth2, sdepth2]), HBox([sdist, sfak]), HBox([psi_x, update, label1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* **Width:** the width of the quantum well.\n",
    "* **Depth:** the depth of the quantum well.\n",
    "* **Zoom factor:** the zoom factor of the eigenfunctions. \n",
    "* **Gap distance:** the distance between two quantum wells (edge to edge)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"height:1px;border:none;color:#cccccc;background-color:#cccccc;\" />\n",
    "\n",
    "## **Legend**\n",
    "\n",
    "(How to use the interactive visualization)\n",
    "\n",
    "### Interactive figures \n",
    "\n",
    "There are two subplots shown above. The wide figure on the left shows \n",
    "the well potential and the square moduli of the eigenfunctions, $|\\psi|^2$, of different \n",
    "states (or just their wavefunctions if the \"Wavefunction\" radio button is selected). The narrow figure on the right shows the eigenvalues.\n",
    "One can highlight the square modulus of a particular eigenfunction, $|\\psi|^2$, along with its corresponding eigenvalue \n",
    "by clicking on the plot. All the other states will \n",
    "be hidden from the plot. \n",
    "\n",
    "### Controls\n",
    "\n",
    "There are sliders to adjust the depth and width of the quantum well potentials. \n",
    "One can also tune the distance between the two quantum wells. \n",
    "By clicking the \"Show all\" button, one can see all states again if they were previously hidden by clicking on one of them to highlight it.\n",
    "The radio button is used to choose whether to show the wavefunctions $\\psi$ or the \n",
    "probability densities $|\\psi|^2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.0"
  },
  "voila": {
   "authors": "Dou Du and Giovanni Pizzi"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}