{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Avoided crossing example\n", "\n", "**Authors:** Dou Du, Taylor James Baird and Giovanni Pizzi \n", "\n", " Go back to index\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", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## **Goals**\n", " * Familiarize oneself 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", " * Appreciate the connection between the presence of regions of avoided crossing in potential energy surfaces and the breakdown of the Born-Oppenheimer approximation for the system being described. \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?\n", "Will the order of the states change by changing the $\\mu$ parameter? And why?\n", "\n", "
\n", " Solution\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. The $\\mu x$ term can be considered as a perturbation\n", " to the original double well potential, which results in the avoided crossing. Please check the detailed\n", " information in the background theory section.\n", "
\n", " \n", "2. How about other states? Is there any changing of energy order of the states?\n", "\n", "
\n", " Solution\n", " By tuning the $\\mu$ slider we observe that the values of the eigenvalues change continuously with the size of the perturbation. However, the order of all the\n", " plotted states remain unchanged (monitor the top right subplot).\n", "
\n", " \n", "3. What is the connection between the phenomenon of avoided crossing and the Born–Oppenheimer approximation?\n", "\n", "
\n", " Solution\n", " The Born–Oppenheimer (BO) approximation is a common approximation employed in quantum chemistry and quantum physics calculations for which one can justify the assumption that the system of interest remains in its electronic ground state. In the\n", " BO approximation, the non-adiabatic coupling vector linking separate electronic states needs to be\n", " negligible:\n", " $$F_{ij}^k(R)=\\frac{<\\psi_i(r;R)|\\nabla_{R_k}H^e|\\psi_j(r;R)>}{V_j-V_i}$$\n", " where $V_j$ and $V_i$ are the eigenvalues of the electronic\n", " Shrödinger equation. If $V_j$ and $V_i$ have close values,\n", " the term is not negligible anymore and BO approximation breaks down. Consequently, we may observe that in the region of an avoided crossing, the validity of the Born-Oppenheimer approximation deteriorates.\n", " A more detailed discussion of the BO approximation can be found at \n", " \n", " Wikipedia.\n", "
\n", " \n", "4. What type of molecular system could be described by the asymmetric double well model we consider here? \n", "\n", "
\n", " Solution\n", " As demonstrated in the notebook 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 Wikipedia).\n", "
\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\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 = [];\n", " ax2.lines = [];\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 = [];\n", " ax2.lines = [];\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": [ "
\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. 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.9.7" }, "voila": { "authors": "Dou Du and Giovanni Pizzi" } }, "nbformat": 4, "nbformat_minor": 4 }