notebook/quantum-mechanics/asymmetricwell.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# **Avoided Crossing in a One-Dimensional Asymmetric Quantum Well**\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/asymmetricwell.ipynb\n",
    "\n",
    "In this notebook we demonstrate the phenomenon of avoided crossing by solving the Shrödinger\n",
    "equation of a one-dimensional asymmetric quantum well.\n",
    "<hr style=\"height:1px;border:none;color:#cccccc;background-color:#cccccc;\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **Goals**\n",
    "  * Familiarize yourself with the phenomenon of avoided crossing.\n",
    "  * Understand the mathematical origin of avoided crossing (i.e. how the behaviour arises from the description of a given system through a dependence on some set of parameters).\n",
    "  * Relate the observation of avoided crossing in this toy model to its manifestation in realistic molecular and material systems."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **Background theory**\n",
    "\n",
    "[**More in the background theory.**](./theory/theory_asymmetricwell.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **Tasks and exercises**\n",
    "\n",
    "1. For which value of $\\mu$ are the two lowest eigenvalues closest to each other? Will the order of the states change by changing the $\\mu$ parameter? And why?\n",
    "    <details>\n",
    "    <summary style=\"color: red\">Solution</summary>\n",
    "    <div style=\"border:blue; border-width:3px; border-style:outset;\">\n",
    "\n",
    "    In the figure, the blue and red lines show the lowest and \n",
    "    second-lowest eigenvalues respectively. The subplot on the bottom shows\n",
    "    the two lowest eigenvalues as a function of the parameter $\\mu$.\n",
    "    One can see that these two eigenvalues are closest to each other\n",
    "    at $\\mu = 0$. For the entire range of plotted $\\mu$ values, the red line is always higher\n",
    "    than the blue line, i.e. the eigenvalues change continuously. However, the nature of the state changes as we change the value of mu! For values $\\mu << 0$, the lowest eigenstate is localized in the right well, while for $\\mu >>0$ is is localized in the left well. For values of $\\mu$ close to zero, the two states mix and we get a symmetric and an antisymmetric solution. The $\\mu x$ term can be considered as a perturbation to the original double well potential, which results in the avoided crossing. Please check the detailed\n",
    "    information in the background theory section.\n",
    "    </div>        \n",
    "    </details>\n",
    "    \n",
    "2. How about other states? Is there any changing of energy order of the states?\n",
    "    <details>\n",
    "    <summary style=\"color: red\">Solution</summary>\n",
    "    <div style=\"border:blue; border-width:3px; border-style:outset;\">\n",
    "    By tuning the $\\mu$ slider we observe that the values of the eigenvalues change continuously with the size of the perturbation. Also in this case, when the value of $\\mu$ would tend to make two eigenvalues identical, an anticrossing effect takes place (e.g. between the second and third state at $\\mu \\sim 0.085$), and we see the appearance of a symmetric and an antisymmetric eigenstate.  \n",
    "    </div>\n",
    "    </details>\n",
    "    \n",
    "3. What type of molecular system could be described by the asymmetric double well model we consider here? \n",
    "\n",
    "    <details>\n",
    "    <summary style=\"color: red\">Solution</summary>\n",
    "    <div style=\"border:blue; border-width:3px; border-style:outset;\">\n",
    "\n",
    "       As demonstrated in the <a href=\"../2quantumwells.ipynb\">notebook</a>, illustrating a double quantum well system, this type of potential gives a reasonable description of diatomic molecules. Adding in asymmetry simply accounts for the two atoms comprising the molecule being different (a heteronuclear molecule). This toy model therefore illustrates the fact that, in diatomic systems, which are parameterized by a single value - namely the separation of the two nuclei - the eigenvalues of the system never cross. One must consider polyatomic molecules consisting of three or more atoms before crossing of electronic energy levels is observed (see <a href=\"https://en.wikipedia.org/wiki/Conical_intersection\">Wikipedia</a>).\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, concatenate\n",
    "from scipy.linalg import eigh_tridiagonal\n",
    "from ipywidgets import FloatSlider, jslink, VBox, HBox, Button, Label, RadioButtons\n",
    "from matplotlib.gridspec import GridSpec\n",
    "import matplotlib.pyplot as plt\n",
    "from math import pi\n",
    "from numpy import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors = ['b', 'r', 'g', 'c', 'm', 'y'] \n",
    "\n",
    "def potential(x, mu):\n",
    "    \"\"\"Potential function for double trough, mu: Parameter\"\"\"\n",
    "    return x**4 - 0.6*x**2 + mu*x # asymmetrical double trough\n",
    "\n",
    "def diagonalization(hquer, L, N, pot=potential, mu = 0.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",
    "         mu: potential function parameter \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, mu)\n",
    "    z = hquer**2 /2.0/dx**2                       \n",
    "    h = (diag(V+2.0*z) + diag(-z*ones(N-1), -1)   \n",
    "                      + diag(-z*ones(N-1), 1) )   #  Hamilton-Matrix\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.2, fak= 2.0):\n",
    "    \"\"\"Plot der niedrigsten Eigenfunktionen 'ef' im Potential 'V(x)'\n",
    "       auf Hoehe der Eigenwerte 'ew' in den Plotbereich 'ax'.\n",
    "       \n",
    "       Der optionale Parameter 'width' (mit Defaultwert 1)\n",
    "       gibt die Linienstaerke beim Plot der Eigenfunktionen\n",
    "       an. 'width' kann auch ein Array von Linienstaerken sein.\n",
    "       'Emax' (mit Default-Wert V_0/10) legt die Energieobergrenze\n",
    "       fuer den Plot fest.\n",
    "       'fak' ist ein Skalierungsfaktor fuer die graphische Darstellung\n",
    "       der Eigenfunktionen.\n",
    "    \"\"\"\n",
    "    fak = fak/100.0; \n",
    "\n",
    "    ax.set_xlim([min(x), max(x)])\n",
    "    ax.set_ylim([min(V)-0.01, ew[5]-0.02])\n",
    "    ax2.set_ylim([min(V)-0.01, ew[5]-0.02])\n",
    "    \n",
    "    #ax[1].yaxis.set_label_position('right')\n",
    "    #ax[1].yaxis.tick_right()\n",
    "    ax.set_xticks([-1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5])\n",
    "    ax.set_xticklabels([-1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5])\n",
    "    ax.set_xlabel(r'$x/a$', fontsize = 10)\n",
    "    ax.set_ylabel(r'$V(x)/V_0\\ \\rm{, Eigenfunctions}$', fontsize = 10)\n",
    "      \n",
    "    indmax = sum(ew<=Emax)                       \n",
    "    \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.plot(x, fak*(ef[:, i])+ew[i], linewidth=width[i]+.1, color=colors[i%len(colors)])\n",
    "        else:\n",
    "            ax.plot(x, fak*abs(ef[:, i])**2+ew[i], linewidth=width[i]+.1, color=colors[i%len(colors)])\n",
    "        ax2.plot(x, x*0.0+ew[i], linewidth=width[i]+2.5, color=colors[i%len(colors)])\n",
    "    \n",
    "    ax.plot(x, V, c='k', linewidth=1.3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu = -0.1                                            # Potential parameter\n",
    "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                                        # The width of the Guassian function\n",
    "zeiten = linspace(0.0, 10.0, 400)                    # Time\n",
    "\n",
    "\n",
    "smu = FloatSlider(value = -0.1, min = -0.1, max = 0.1, step = 0.01, description = r'$\\mu$: ')\n",
    "psi_x = RadioButtons(options=[\"Wavefunction\", \"Probability density\"], value=\"Wavefunction\", description=\"Show:\")\n",
    "\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "fig.canvas.header_visible = False\n",
    "\n",
    "gs = GridSpec(10, 6, figure=fig)\n",
    "\n",
    "ax1 = fig.add_subplot(gs[0:6, 0:5])\n",
    "ax2 = fig.add_subplot(gs[0:6, 5])\n",
    "ax3 = fig.add_subplot(gs[7:, :])\n",
    "\n",
    "fig.suptitle('Avoided Crossing in 1D Quantum Well', fontsize = 13)\n",
    "\n",
    "mu1 = []\n",
    "ew1 = []\n",
    "ew2 = []\n",
    "ew3 = []\n",
    "\n",
    "\n",
    "for i in concatenate((linspace(-0.1, -0.01, 10), linspace(-0.01, 0.01, 20), linspace(0.01, 0.1, 10))):\n",
    "    ew, ef, x, dx, V = diagonalization(hquer, L, N, mu = i)\n",
    "    mu1.append(i)\n",
    "    ew1.append(ew[0])\n",
    "    ew2.append(ew[1])\n",
    "    ew3.append(ew[2])\n",
    "\n",
    "ew, ef, x, dx, V = diagonalization(hquer, L, N, mu = mu)\n",
    "\n",
    "ax3.plot(mu1, ew1, c='b')\n",
    "ax3.plot(mu1, ew2, c='r')\n",
    "ax3.plot(mu1, ew3, c='g')\n",
    "\n",
    "s3 = ax3.plot(mu, ew[2], 'go', label = str(format(ew[0], '.3f')))\n",
    "s2 = ax3.plot(mu, ew[1], 'ro', label = str(format(ew[1], '.3f')))\n",
    "s1 = ax3.plot(mu, ew[0], 'bo', label = str(format(ew[0], '.3f')))\n",
    "\n",
    "ax3.legend()\n",
    "\n",
    "ax3.set_xlim([min(mu1), max(mu1)])\n",
    "ax3.set_ylim([min(ew1), max(ew3)+0.005])\n",
    "\n",
    "ax3.set_xlabel(r'Potential parameter $\\mu$ value', fontsize = 10)\n",
    "#ax3.set_ylabel(r'The values of the two lowest eigenvalues', fontsize = 10)\n",
    "ax3.set_ylabel(r\"E\",fontsize = 10)\n",
    "ax3.set_xticks([-0.1, -0.05, 0.0, 0.05, 0.1])\n",
    "ax3.set_xticklabels([-0.1, -0.05, 0.0, 0.05, 0.1])\n",
    "\n",
    "ax2.yaxis.set_label_position(\"right\")\n",
    "ax2.yaxis.tick_right()\n",
    "ax2.get_xaxis().set_visible(False)\n",
    "ax2.set_ylabel(r'$\\rm{\\ Eigenvalues}$', fontsize = 10)\n",
    "\n",
    "plot_eigenfunctions(ax1, ew, ef, x, V)\n",
    "\n",
    "plt.show()\n",
    "\n",
    "sfak = FloatSlider(value = 2, min = 1.0, max = 5.0, step = 1.0, description = r'Zoom factor: ')\n",
    "\n",
    "def on_mu_change(change):\n",
    "    global ew, ef, x, dx, V\n",
    "    ax1.lines.clear()\n",
    "    ax2.lines.clear()\n",
    "    \n",
    "    try:\n",
    "        ann.remove()\n",
    "        ann1.remove()\n",
    "    except:\n",
    "        pass\n",
    "    \n",
    "    ew, ef, x, dx, V = diagonalization(hquer, L, N, mu = smu.value)\n",
    "    plot_eigenfunctions(ax1, ew, ef, x, V, fak = sfak.value)\n",
    "    \n",
    "    ax3.lines.pop(-1)\n",
    "    ax3.lines.pop(-1)\n",
    "    ax3.lines.pop(-1)\n",
    "\n",
    "    ax3.plot(smu.value, ew[0], 'bo', label = str(format(ew[0], '.3f')))\n",
    "    ax3.plot(smu.value, ew[1], 'ro', label = str(format(ew[1], '.3f')))\n",
    "    ax3.plot(smu.value, ew[2], 'go', label = str(format(ew[2], '.3f')))\n",
    "    ax3.legend()\n",
    "    \n",
    "    for i in ax3.lines:\n",
    "        i.set_alpha(1.0)\n",
    "\n",
    "smu.observe(on_mu_change, names = 'value')\n",
    "psi_x.observe(on_mu_change, names = 'value')\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(ax2.lines)):\n",
    "        ax1.lines[i].set_alpha(0.1)\n",
    "        ax2.lines[i].set_alpha(0.1)\n",
    "        ax1.lines[i].set_linewidth(1.1)\n",
    "        \n",
    "    ax1.lines[ixx].set_alpha(1.0)\n",
    "    ax2.lines[ixx].set_alpha(1.0)\n",
    "    ax1.lines[ixx].set_linewidth(2.0)\n",
    "    \n",
    "    for i in range(3):\n",
    "        if i == ixx:\n",
    "            ax3.lines[i].set_alpha(1.0)\n",
    "        else:\n",
    "            ax3.lines[i].set_alpha(0.1)\n",
    "    \n",
    "    try:\n",
    "        ann.remove()\n",
    "        ann1.remove()\n",
    "    except:\n",
    "        pass\n",
    "    \n",
    "    ann = ax1.annotate(s = 'n = ' + str(ixx+1), xy = (0, ew[ixx]), xytext = (-0.15, ew[ixx]), xycoords = 'data', color='k', size=15)\n",
    "    ann1 = ax2.annotate(s = str(\"{:.3f}\".format(ew[ixx])), xy = (0, ew[ixx]), xytext = (-1.2, ew[ixx]+0.005), xycoords = 'data', color='k', size=11)\n",
    "\n",
    "cid = fig.canvas.mpl_connect('button_press_event', on_press)\n",
    "    \n",
    "def on_xfak_change(change):\n",
    "    ax1.lines.clear()\n",
    "    ax2.lines.clear()\n",
    "    plot_eigenfunctions(ax1, ew, ef, x, V, fak = sfak.value)\n",
    "\n",
    "sfak.observe(on_xfak_change, names = 'value')\n",
    "\n",
    "def on_update_click(b):\n",
    "    for i in ax1.lines:\n",
    "        i.set_alpha(1.0)\n",
    "        i.set_linewidth(1.1)\n",
    "    \n",
    "    for i in ax2.lines:\n",
    "        i.set_alpha(1.0)\n",
    "        \n",
    "    for i in ax3.lines:\n",
    "        i.set_alpha(1.0)\n",
    "        \n",
    "    try:\n",
    "        ann.remove()\n",
    "        ann1.remove()\n",
    "    except:\n",
    "        pass\n",
    "    \n",
    "\n",
    "update = Button(description=\"Show all\")\n",
    "update.on_click(on_update_click)\n",
    "\n",
    "label1 = Label(value=\"(click on a state to select it)\");\n",
    "\n",
    "display(HBox([smu, sfak]), HBox([psi_x, update, label1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* **$\\mu$:** the potential parameter determining the symmetry of the double well potential.\n",
    "* **Zoom factor:** the zoom factor of the eigenfunctions. "
   ]
  },
  {
   "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. In the uppermost subplot, the wide figure on the left shows \n",
    "the well potential alongside the eigenfunctions, $\\psi$, of different states (or their corresponding probability density, the square modulus $|\\psi|^2$). The narrow figure on the right shows \n",
    "the corresponding eigenvalues. \n",
    "\n",
    "The subplot on the bottom shows how the three lowest eigenvalues change with the modification of the potential parameter, $\\mu$.\n",
    "\n",
    "\n",
    "## Controls\n",
    "\n",
    "There is a slider to adjust the $\\mu$ parameter of the potential.\n",
    "The zoom factor slider aids in the inspection of the wavefunctions by\n",
    "multiplying them by a constant (this is purely for visualization purposes). One can highlight the square modulus of a specific eigenfunction, $\\psi$, and its\n",
    "eigenvalue, by clicking on the plot. All the other \n",
    "states shall be hidden from the plot. By clicking the `Show all` button, one can \n",
    "see all states again if they were previously hidden by clicking on one of them to highlight it. Additionally, there is a radio button which enables one to choose whether to plot the wavefunctions $\\psi$, or the 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
}