notebook/band-theory/pseudopotential.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# **Norm-Conserving Pseudopotentials**\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/band-theory/pseudopotential.ipynb\n",
    "\n",
    "The pseudopotential method is a technique employed to simplify the description of a system of interacting electrons and nuclei. It is used to construct an effective potential\n",
    "that includes both the effects of a nucleus and of its core electrons,\n",
    "allowing one to consider explicitly only the valence electrons.\n",
    "This notebook illustrates a method of constructing norm-conserving pseudopotentials\n",
    "and to display them interactively, together with the resulting pseudowavefunctions.\n",
    "\n",
    "<hr style=\"height:1px;border:none;color:#cccccc;background-color:#cccccc;\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "## **Goals**\n",
    "\n",
    "* Understand why pseudopotentials are needed.\n",
    "* Learn how to construct pseudopotentials using Kerker's method.\n",
    "* Examine the results for various values of the principal quantum number <em>n</em> and of the angular quantum number <em>l</em>.\n",
    "* Examine the effect of changing the cutoff radius."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **Background theory** \n",
    "\n",
    "[More on the background theory.](./theory/theory_pseudopotential.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## **Tasks and exercises**\n",
    "\n",
    "1. Investigate the role of the cutoff radius by varying the $R_c$ slider.\n",
    "    \n",
    "    <details>\n",
    "    <summary style=\"color: red\">Solution</summary>\n",
    "        Move the slider for $R_c$ and press the button \"Compute pseudopotential\" to obtain \n",
    "        the results. Check if there are values for which no solutions can be found. \n",
    "        Inspect how different the pseudopotential is from the Coulomb potential.\n",
    "    </details>\n",
    "\n",
    "2. Investigate how the pseudopotential changes for different values of the quantum numbers\n",
    "    \n",
    "    <details>\n",
    "    <summary style=\"color: red\">Solution</summary>\n",
    "        Try to construct the pseudopotential for various values of <em>n</em> and <em>l</em>.\n",
    "        Check what happens when constructing a pseudopotential for a nodeless wavefunction \n",
    "        (e.g. $n=1$ and $l=0$, or $n=2$ and $l=1$).\n",
    "    </details>\n",
    "\n",
    "3. Why do we need pseudopotentials?\n",
    "    \n",
    "    <details>\n",
    "    <summary style=\"color: red\">Solution</summary>\n",
    "        Wavefunction oscillates rapidly in the core region. In a plane-wave approach, \n",
    "        this would require a huge basis set (i.e., a huge number of plane waves) to \n",
    "        be described accurately. What is most relevant, however, is that while the \n",
    "        largest part of the contribution to the total energy of the system comes from \n",
    "        the core electrons, these electrons are essentially frozen and do not participate \n",
    "        in chemistry and the creation of bonds; the electronic structure is instead \n",
    "        determined by the valence electrons. Avoiding to treat explicitly core \n",
    "        electrons avoids that small relative errors on the core electrons completely \n",
    "        spoil the calculation of the energy of the full system and, in particular, \n",
    "        of the (small, but crucial) energy differences between different atomic \n",
    "        configurations or crystalline phases.\n",
    "    </details>\n",
    "\n",
    "4. What is the meaning of the norm-conservation condition?\n",
    "    \n",
    "    <details>\n",
    "    <summary style=\"color: red\">Solution</summary>\n",
    "        The condition ensures that the total charge inside the cutoff radius $R_c$ is \n",
    "        correct. However, there are more profound consequences that are implied by \n",
    "        this condition: it turns out that imposing norm conservation implies also \n",
    "        that the first energy derivative of the logarithmic derivatives of the \n",
    "        all-electron wavefunction and the pseduowafevuntion agree at $R_c$. This \n",
    "        is a very important condition for the transferability of the pseudopotential. \n",
    "        It can be shown that this means that the true atom with its electrons and the\n",
    "        pseudopotential generate the same phase shift when a plane wave is scattered \n",
    "        into a spherical wave. A detailed discussion can be found in Section 11.4 of \n",
    "        the book <a href=\"https://www.cambridge.org/ch/academic/subjects/physics/condensed-matter-physics-nanoscience-and-mesoscopic-physics/electronic-structure-basic-theory-and-practical-methods?format=PB&isbn=9780521534406\" target=\"_blank\">\"Electronic Structure: \n",
    "        Basic Theory and Practical Methods\"</a> by Richard M. Martin, Cambridge \n",
    "        University Press (2004).\n",
    "    </details>\n",
    "\n",
    "5. Are the pseudopotentials constructed with this method local?\n",
    "    \n",
    "    <details>\n",
    "    <summary style=\"color: red\">Solution</summary>\n",
    "        As discussed earlier, we obtain a different pseudopotential for different \n",
    "        values of the quantum numbers <em>n</em> and <em>l</em>. Therefore, the \n",
    "        pseudopotential is not local (its action is not not just the product of \n",
    "        the same a single function $V^{PS}(r)$ times the wavefunction).\n",
    "    </details>\n",
    "\n",
    "<hr style=\"height:1px;border:none;color:#cccccc;background-color:#cccccc;\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interactive visualization\n",
    "(be patient, it might take a few seconds to load)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy.physics.hydrogen import R_nl, E_nl\n",
    "from sympy.abc import r\n",
    "from sympy.functions import exp\n",
    "from sympy import lambdify, diff, log\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sympy.solvers.solveset import linsolve\n",
    "from ipywidgets import FloatSlider, Button, IntSlider, Layout, HBox, VBox, Label, Tab, Layout, Text\n",
    "from scipy.optimize import newton\n",
    "import matplotlib.gridspec as gridspec\n",
    "\n",
    "%matplotlib widget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sn: the slider to control the quantum number n\n",
    "# sl: the slider to control the quantum number l\n",
    "sn = IntSlider(value=3, min=1, max=5, description=\"$n$\")\n",
    "sl = IntSlider(value=0, min=0, max=sn.value-1, description=\"$l$\")\n",
    "\n",
    "labeln = Label(value=\"(principal quantum number)\")\n",
    "labell = Label(value=\"(angular quantum number)\")\n",
    "\n",
    "#init the quantum number n and l\n",
    "n = sn.value\n",
    "l = sl.value\n",
    "# We consider a H atom, with charge +1\n",
    "Z = 1\n",
    "\n",
    "#ho: is the Hydrogen radial wavefunction\n",
    "#Ea: is the eigenvalue\n",
    "#rf: make the analytical formula numerically\n",
    "ho = R_nl(n, l, r, Z=Z)\n",
    "Ea = E_nl(n, Z=Z)\n",
    "rf = lambdify(r, ho, \"numpy\")\n",
    "\n",
    "\n",
    "def update_functions_and_plot():\n",
    "    \"\"\"Update the functions and then the plot.\n",
    "    \"\"\"\n",
    "    global ho, Ea, rf\n",
    "    ho = R_nl(n, l, r, Z=Z)\n",
    "    # Check that the value of the radial wavefunction at large radius is positive\n",
    "    # (i.e., that the wavefunction goes to zero from positive values as r->infinity)\n",
    "    # If negative: flip it! (This is because due to the formula of the polynomial\n",
    "    # that we want to use for r < R_c, $\\Psi(R_c)R_c$ must be positive at R_c, and thus after R_c as well).\n",
    "    # Because of the form of R_nl, the sign needs to be swapped if n-l is even.\n",
    "    if (n - l) % 2 == 0:\n",
    "        ho = -ho\n",
    "    Ea = E_nl(n, Z=Z)\n",
    "    rf = lambdify(r, ho, \"numpy\")\n",
    "    update_plot()\n",
    "\n",
    "def nvalue_change(c):\n",
    "    \"\"\"Observe the change of the sn and update plot.\n",
    "    \"\"\"\n",
    "    global n\n",
    "    n = c[\"new\"]\n",
    "    sl.value = 0\n",
    "    sl.max = n - 1\n",
    "    update_functions_and_plot()\n",
    "    \n",
    "sn.observe(nvalue_change, names=\"value\")\n",
    "\n",
    "def lvalue_change(c):\n",
    "    \"\"\"Observe the change of the sl and upate plot.\n",
    "    \"\"\"\n",
    "    global l\n",
    "    l = c[\"new\"]\n",
    "    update_functions_and_plot()\n",
    "\n",
    "sl.observe(lvalue_change, names=\"value\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## NOTE: the widgets in this cell are currently not shown\n",
    "\n",
    "text_l0 = Text(description = r\"Lt. $\\psi^{(0)}(r_c)$:\")\n",
    "text_l1 = Text(description = r\"Lt. $\\psi^{(1)}(r_c)$:\")\n",
    "text_l2 = Text(description = r\"Lt. $\\psi^{(2)}(r_c)$:\")\n",
    "\n",
    "text_r0 = Text(description = r\"Rt. $\\psi^{(0)}(r_c)$:\")\n",
    "text_r1 = Text(description = r\"Rt. $\\psi^{(1)}(r_c)$:\")\n",
    "text_r2 = Text(description = r\"Rt. $\\psi^{(2)}(r_c)$:\")\n",
    "\n",
    "cof_0 = Text(description = r\"$\\lambda_0$:\")\n",
    "cof_2 = Text(description = r\"$\\lambda_2$:\")\n",
    "cof_3 = Text(description = r\"$\\lambda_3$:\")\n",
    "cof_4 = Text(description = r\"$\\lambda_4$:\")\n",
    "\n",
    "def clear_texts():\n",
    "    text_l0.value = \"\"\n",
    "    text_l1.value = \"\"\n",
    "    text_l2.value = \"\"\n",
    "    text_r0.value = \"\"\n",
    "    text_r1.value = \"\"\n",
    "    text_r2.value = \"\"\n",
    "    cof_0.value = \"\"\n",
    "    cof_2.value = \"\"\n",
    "    cof_3.value = \"\"\n",
    "    cof_4.value = \"\"\n",
    "\n",
    "output1 = VBox([HBox([text_l0, text_r0]), HBox([text_l1, text_r1]), HBox([text_l2, text_r2])]);\n",
    "output2 = VBox([cof_0, cof_2, cof_3, cof_4])\n",
    "\n",
    "tab = Tab(layout=Layout(width='700px'))\n",
    "tab.children = [output1, output2]\n",
    "tab.set_title(0, r\"Continuity at Rc\")\n",
    "tab.set_title(1, r\"Polynomial coeff.\")\n",
    "\n",
    "## Uncomment this if you want to also see some text boxes with the values of the wavefunction and its derivatives\n",
    "## at the cutoff radius\n",
    "#display(tab)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_rc = FloatSlider(value = 21.0, min = 1.0, max = 50, description = \"$R_c$\", \n",
    "                   layout={'width':'300px'})\n",
    "\n",
    "compute = Button(description=\"Compute pseudopotential\", style={'description_width': 'initial'}, layout={'width':'200px'})\n",
    "\n",
    "\n",
    "display(HBox([sn, labeln]), HBox([sl, labell]))\n",
    "display(HBox([s_rc, compute]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "img = plt.figure(tight_layout=True, figsize=(9,9))\n",
    "img.canvas.header_visible = False\n",
    "\n",
    "gs = gridspec.GridSpec(2, 2)\n",
    "\n",
    "ax1 = img.add_subplot(gs[0, 0])\n",
    "ax2 = img.add_subplot(gs[0, 1])\n",
    "ax3 = img.add_subplot(gs[1, 0]) # subplot for log derivatives\n",
    "ax4 = img.add_subplot(gs[1, 1])\n",
    "\n",
    "x1 = np.arange(0, 50.0, 0.01)\n",
    "y1 = rf(x1)*x1\n",
    "\n",
    "line_rho, = ax1.plot(x1, y1, 'r-', label=\"$R_{\"+str(sn.value)+str(sl.value)+\"}(r)r$\")\n",
    "#ax1.fill_between(x1, y1, 0, where=x1<s_rc.value, facecolor='yellow', alpha=0.5)\n",
    "ax1.set_xlim([0, 50.0])\n",
    "\n",
    "ax1.hlines(0, 0, 50, 'k','--')\n",
    "\n",
    "\n",
    "# subplot for R^2r^2\n",
    "\n",
    "x2 = np.arange(0, 50.0, 0.01)\n",
    "y2 = (rf(x2)*x2)**2\n",
    "\n",
    "line_rho2, = ax2.plot(x2, y2, 'r-', label=\"$R^2_{\"+str(sn.value)+str(sl.value)+\"}(r)r^2$\")\n",
    "\n",
    "ax2.fill_between(x2, y2, 0, where=x2<s_rc.value, facecolor='yellow', alpha=0.5)\n",
    "ax2.set_xlim([0, 50.0])\n",
    "\n",
    "ax2.hlines(0, 0, 50, 'k','--')\n",
    "\n",
    "line_rc1 = ax1.axvline(s_rc.value)\n",
    "line_pswf, = ax1.plot([],[],'b-', linewidth=1.5, label=\"$R^{PS}_{\"+str(sn.value)+str(sl.value)+\"}(r)r$\")\n",
    "ann_rc = ax1.annotate(\"$R_c$\", xy=(s_rc.value + 1.0, rf(s_rc.value)*s_rc.value), fontsize=20)\n",
    "point, = ax1.plot(s_rc.value, rf(s_rc.value)*s_rc.value, 'ko')\n",
    "\n",
    "line_rc2 = ax2.axvline(s_rc.value)\n",
    "line_pswf2, = ax2.plot([],[],'b-', linewidth=1.5, label=\"$R^{PS}_{\"+str(sn.value)+str(sl.value)+\"}(r)r$\")\n",
    "ann_rc2 = ax2.annotate(\"$R_c$\", xy=(s_rc.value + 1.0, (rf(s_rc.value)*s_rc.value)**2), fontsize=20)\n",
    "point2, = ax2.plot(s_rc.value, (rf(s_rc.value)*s_rc.value)**2, 'ko')\n",
    "\n",
    "ann_norm1 = ax2.annotate(\"Yellow area:\", xy=(50, 30), xycoords='axes points', fontsize=9, \n",
    "                         bbox=dict(boxstyle='round', facecolor='yellow', alpha=0.5))\n",
    "ann_norm2 = ax2.annotate(\"Green area:\", xy=(50, 10), xycoords='axes points', fontsize=9,\n",
    "                         bbox=dict(boxstyle='round', facecolor='green', alpha=0.5))\n",
    "\n",
    "# ann_logl = ax3.annotate(\"Logarithmic deriv. $\\psi(r_c)^{PS}$:\", xy=(0,0), xycoords='axes points', fontsize=9, \n",
    "#                          bbox=dict(boxstyle='round', facecolor='blue', alpha=0.5))\n",
    "# ann_logr = ax3.annotate(\"Logarithmic deriv. $\\psi(r_c)$:\", xy=(0,0), xycoords='axes points', fontsize=9,\n",
    "#                          bbox=dict(boxstyle='round', facecolor='red', alpha=0.5))\n",
    "\n",
    "# ann_norm1 = ax2.annotate(\"\", xy=(250, 30), xycoords='axes points', fontsize=9, \n",
    "#                          bbox=dict(boxstyle='round', facecolor='yellow', alpha=0.5))\n",
    "# ann_norm2 = ax2.annotate(\"\", xy=(250, 10), xycoords='axes points', fontsize=9,\n",
    "#                          bbox=dict(boxstyle='round', facecolor='green', alpha=0.5))\n",
    "\n",
    "# ann_logl = ax3.annotate(\"\", xy=(0,0), xycoords='axes points', fontsize=9, \n",
    "#                          bbox=dict(boxstyle='round', facecolor='blue', alpha=0.5))\n",
    "# ann_logr = ax3.annotate(\"\", xy=(0,0), xycoords='axes points', fontsize=9,\n",
    "#                          bbox=dict(boxstyle='round', facecolor='red', alpha=0.5))\n",
    "\n",
    "ax1.set_xlabel(\"r\", fontsize = 15)\n",
    "ax1.set_ylabel(\"$R_{nl}(r)r$\", fontsize = 15)\n",
    "\n",
    "ax1.legend(loc=1, fontsize=15)\n",
    "\n",
    "ax2.set_xlabel(\"r\", fontsize = 15)\n",
    "ax2.set_ylabel(\"$R_{nl}(r)^2r^2$\", fontsize = 15)\n",
    "\n",
    "ax2.legend(loc=1, fontsize=15)\n",
    "#ax1.set_ylim([-0.3, 0.5])\n",
    "\n",
    "x3 = np.arange(0, 50.0, 0.01)\n",
    "line_rc3 = ax3.axvline(s_rc.value)\n",
    "line_logdr_rho, = ax3.plot([],[],'r-', linewidth=1.5, label=r\"$\\frac{d}{dr} \\ln(R_{\"+str(sn.value)+str(sl.value)+r\"}(r)r)$\")\n",
    "line_logdr_ps, = ax3.plot([],[],'b-', linewidth=1.5, label=r\"$\\frac{d}{dr} \\ln (R^{PS}_{\"+str(sn.value)+str(sl.value)+r\"}(r)r)$\")\n",
    "# line_logdr_rho, = ax3.plot([],[],'r-', linewidth=1.5, label=r\"$\\frac{d}{dr} \\ln \\psi$\"+\"${}_{\"+str(sn.value)+\"}$\"+\"${}_{\"+str(sl.value)+\"}$\")\n",
    "# line_logdr_ps, = ax3.plot([],[],'b-', linewidth=1.5, label=r\"$\\psi^{PS}$\")\n",
    "\n",
    "\n",
    "ax3.set_xlim([0, 50.0])\n",
    "\n",
    "ax3.hlines(0, 0, 50, 'k','--')\n",
    "point3, = ax3.plot( [],[], 'ko')\n",
    "ann_rc3 = ax3.annotate(\"$R_c$\", xy=(20,100) , fontsize=20)\n",
    "# ann_rc3.set_visible(False)\n",
    "\n",
    "ax3.set_xlabel(\"r\", fontsize = 15)\n",
    "ax3.set_ylabel(\"Log deriv.\", fontsize = 15)\n",
    "\n",
    "ax3.legend(loc=1, fontsize=15)\n",
    "\n",
    "x4 = np.linspace(0.001, 50, 500)\n",
    "y4 = -Z/x4\n",
    "\n",
    "ax4.plot(x4, y4, 'r--', linewidth=1.0, label=\"$-Z/r$\")\n",
    "\n",
    "ax4.set_xlim([0, 50.0])\n",
    "#ax4.set_ylim([-1.0, 0.10])\n",
    "line_rc4 = ax4.axvline(s_rc.value)\n",
    "line_psv, = ax4.plot([],[], 'b-', linewidth=1.5, label=\"$V^{PS}(r)$\")\n",
    "ax4.hlines(0, 0, 50, 'k','--')\n",
    "\n",
    "ax4.set_xlabel(\"r\", fontsize = 15)\n",
    "ax4.set_ylabel(\"$V$\", fontsize = 15)\n",
    "\n",
    "ax4.legend(loc=1, fontsize=15)\n",
    "\n",
    "def update_plot():\n",
    "    \"\"\" Update the plot when quantum number n and l changing.\n",
    "    \"\"\"\n",
    "    x1 = np.arange(0, 50.0, 0.01)\n",
    "    y1 = rf(x1)*x1\n",
    "    line_rho.set_data([x1, y1])\n",
    "    line_rc1.set_data([s_rc.value], [-1, 1])\n",
    "    ann_rc.set_position((s_rc.value + 1.0, rf(s_rc.value)*s_rc.value))\n",
    "\n",
    "    x2 = np.arange(0, 50.0, 0.01)\n",
    "    y2 = (rf(x2)*x2)**2\n",
    "    line_rho2.set_data([x2, y2])\n",
    "    line_rc2.set_data([s_rc.value], [-1, 1])\n",
    "    ann_rc2.set_position((s_rc.value + 1.0, (rf(s_rc.value)*s_rc.value)**2))\n",
    "  \n",
    "\n",
    "    x3 = np.arange(0, 50.0, 0.01)\n",
    "    \n",
    "    # logl_r = lambdify(r, diff(ps,r)/ps, \"numpy\")\n",
    "    # logr_r = lambdify(r, diff(ho,r)/ho, \"numpy\")\n",
    "\n",
    "    # logl = diff(ps, r).subs(r, rc).evalf()/ps.subs(r, rc).evalf()\n",
    "    # logr = diff(ho, r).subs(r, rc).evalf()/ho.subs(r, rc).evalf()\n",
    "    \n",
    "    # line_logdr_rho.set_data(x3,logr_r(x3))\n",
    "    # line_logdr_ps.set_data(x3,logl_r(x3))\n",
    "    \n",
    "    line_logdr_rho.set_data([],[])\n",
    "    line_logdr_ps.set_data([],[])\n",
    "    \n",
    "    # line_logdr_rho.set_data([x3, y3])\n",
    "    line_rc3.set_data([s_rc.value], [-1, 1])\n",
    "    line_rc4.set_data([s_rc.value], [-1, 1])\n",
    "\n",
    "    ann_rc2.set_position((s_rc.value + 1.0, (rf(s_rc.value)*s_rc.value)**2))\n",
    "  \n",
    "    for line in ax1.collections:\n",
    "        line.remove()\n",
    "        \n",
    "    for line in ax2.collections:\n",
    "        line.remove()\n",
    "\n",
    "    for line in ax3.collections:\n",
    "        line.remove()\n",
    "\n",
    "    ax2.fill_between(x2, y2, 0, where=x2<s_rc.value, facecolor='yellow', alpha=0.5)\n",
    "    ax1.hlines(0, 0, 50, 'k','--')\n",
    "    ax2.hlines(0, 0, 50, 'k','--')\n",
    "\n",
    "    line_pswf.set_data([],[])\n",
    "    line_pswf2.set_data([],[])\n",
    "\n",
    "    line_psv.set_data([],[])\n",
    "    point.set_data([s_rc.value], [rf(s_rc.value)*s_rc.value])\n",
    "    point2.set_data([s_rc.value], [(rf(s_rc.value)*s_rc.value)**2])\n",
    "\n",
    "    ax1.set_ylim([y1.min()-0.04, y1.max()+0.04])\n",
    "    ax2.set_ylim([y2.min()-0.04, y2.max()+0.04])\n",
    "\n",
    "    line_rho.set_label(\"$R_{\"+str(sn.value)+str(sl.value)+\"}(r)r$\")\n",
    "    line_pswf.set_label(\"$R^{PS}_{\"+str(sn.value)+str(sl.value)+\"}(r)r$\")\n",
    "    \n",
    "    line_rho2.set_label(\"$R_{\"+str(sn.value)+str(sl.value)+\"}(r)r$\")\n",
    "    line_pswf2.set_label(\"$R^{PS}_{\"+str(sn.value)+str(sl.value)+\"}(r)r$\")\n",
    "    ax1.legend(loc=1, fontsize=15)\n",
    "    ax2.legend(loc=1, fontsize=15)\n",
    "\n",
    "    ax4.legend(loc=1, fontsize=15)\n",
    "    \n",
    "    #clear_texts()\n",
    "\n",
    "def on_rc_change(b):\n",
    "    \"\"\" Update the plot when the slider of the Rc changing.\n",
    "    \"\"\"\n",
    "    x1 = np.arange(0, 50.0, 0.01)\n",
    "    y1 = rf(x1)*x1\n",
    "    x2 = np.arange(0, 50.0, 0.01)\n",
    "    y2 = (rf(x2)*x2)**2\n",
    "    \n",
    "    line_rc1.set_data([s_rc.value], [-1, 1])\n",
    "    line_rc2.set_data([s_rc.value], [-1, 1])\n",
    "    line_rc3.set_data([s_rc.value], [-1, 1])\n",
    "    line_rc4.set_data([s_rc.value], [-1, 1])\n",
    "\n",
    "    ann_rc.set_position((s_rc.value + 1.0, rf(s_rc.value)*s_rc.value))\n",
    "    \n",
    "    for line in ax1.collections:\n",
    "        line.remove()\n",
    "        \n",
    "    for line in ax2.collections:\n",
    "        line.remove()\n",
    "\n",
    "    for line in ax3.collections:\n",
    "        line.remove()\n",
    "\n",
    "\n",
    "    ax2.fill_between(x2, y2, 0, where=x2<s_rc.value, facecolor='yellow', alpha=0.5)\n",
    "    ax2.hlines(0, 0, 50, 'k','--')\n",
    "    line_pswf.set_data([],[])\n",
    "    line_pswf2.set_data([],[])\n",
    "    line_psv.set_data([],[])\n",
    "    point.set_data([s_rc.value], rf(s_rc.value)*s_rc.value)\n",
    "    point2.set_data([s_rc.value], (rf(s_rc.value)*s_rc.value)**2)\n",
    "    ax1.hlines(0, 0, 50, 'k','--')\n",
    "\n",
    "    clear_texts()\n",
    "     \n",
    "s_rc.observe(on_rc_change, names='value')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_right_derivative(rc):\n",
    "    k0 = rf(rc)\n",
    "    k1 = diff(ho, r).subs(r, rc).evalf()\n",
    "    k2 = diff(ho, r, 2).subs(r, rc).evalf()\n",
    "    return np.array([float(k0), float(k1), float(k2)])\n",
    "\n",
    "def compute_left_derivative(ps, rc):\n",
    "    psf = lambdify(r, ps, \"numpy\")\n",
    "    k0 = psf(rc)\n",
    "    k1 = diff(ps, r).subs(r, rc).evalf()\n",
    "    k2 = diff(ps, r, 2).subs(r, rc).evalf()\n",
    "    return np.array([float(k0), float(k1), float(k2)])\n",
    "    \n",
    "def solver_kernel(devs, rc, l, b):\n",
    "    A = np.zeros([3,4])\n",
    "    A[0, :] = np.array([1, rc**2, rc**3, rc**4]);\n",
    "    A[1, :] = np.array([0, 2*rc, 3*rc**2, 4*rc**3]);\n",
    "    A[2, :] = np.array([0, 2, 6*rc, 12*rc**2])\n",
    "    \n",
    "    B = np.zeros(3)\n",
    "    B[0] = log(devs[0]/rc**l);\n",
    "    B[1] = devs[1]/devs[0] - l/rc;\n",
    "    B[2] = devs[2]/devs[0] - devs[1]**2/devs[0]**2 + l/rc**2;\n",
    "    \n",
    "    B-=b*A[:, 1];\n",
    "    A = np.delete(A, (1), axis=1);\n",
    "    \n",
    "    coff = np.linalg.solve(A, B);\n",
    "    coff = np.insert(coff, 1, b)\n",
    "    \n",
    "    return coff\n",
    "\n",
    "def diff_norms(b, rc, l):\n",
    "    devs = compute_right_derivative(rc)\n",
    "    coff = solver_kernel(devs, rc, l, b)\n",
    "    \n",
    "    ps  = r**l*exp(coff[0] + coff[1]*r**2 + coff[2]*r**3 + coff[3]*r**4)\n",
    "    psf = lambdify(r, ps, \"numpy\")\n",
    "    psr = lambdify(r, ps*ps*r*r, \"numpy\")\n",
    "    hor = lambdify(r, ho*ho*r*r, \"numpy\")\n",
    "        \n",
    "    x1 = np.linspace(0, rc, 800);\n",
    "        \n",
    "    norm1 = np.sum(hor(x1))*(x1[1]-x1[0])\n",
    "    norm2 = np.sum(psr(x1))*(x1[1]-x1[0])\n",
    "    \n",
    "    return float(norm1 - norm2)\n",
    "\n",
    "def compute_norms(b, rc, l):\n",
    "    devs = compute_right_derivative(rc)\n",
    "    coff = solver_kernel(devs, rc, l, b)\n",
    "    \n",
    "    ps  = r**l*exp(coff[0] + coff[1]*r**2 + coff[2]*r**3 + coff[3]*r**4)\n",
    "    psf = lambdify(r, ps, \"numpy\")\n",
    "    psr = lambdify(r, ps*ps*r*r, \"numpy\")\n",
    "    hor = lambdify(r, ho*ho*r*r, \"numpy\")\n",
    "        \n",
    "    x1 = np.linspace(0, rc, 800);\n",
    "        \n",
    "    norm1 = np.sum(hor(x1))*(x1[1]-x1[0])\n",
    "    norm2 = np.sum(psr(x1))*(x1[1]-x1[0])\n",
    "    \n",
    "    return norm1, norm2\n",
    "\n",
    "def compute_potential(l):\n",
    "    psf = Ea - l*(l+1)/(2*r*r) + 1/(2*ho*r)*diff(ho*r, r, r)\n",
    "    return lambdify(r, psf, \"numpy\")\n",
    "\n",
    "\n",
    "def plot_ps_wavefunction(b, rc, l):\n",
    "    global ann_rc3\n",
    "    devs = compute_right_derivative(rc)\n",
    "    coff = solver_kernel(devs, rc, l, b)\n",
    "    \n",
    "    ps  = r**l*exp(coff[0] + coff[1]*r**2 + coff[2]*r**3 + coff[3]*r**4)\n",
    "    ps2  = ps**2\n",
    "\n",
    "    psf = lambdify(r, ps*r, \"numpy\")\n",
    "    psf2 = lambdify(r, ps2*r**2, \"numpy\")\n",
    "\n",
    "    devl = compute_left_derivative(ps, rc)\n",
    "    \n",
    "    x1 = np.linspace(0, rc, 800);\n",
    "    x2 = np.linspace(0, rc, 800);\n",
    "    # x3 = np.linspace(0, rc, 800);\n",
    "    x3 = np.arange(0, 50.0, 0.01)\n",
    "\n",
    "    line_pswf.set_data([x1], [psf(x1)]);\n",
    "    line_pswf2.set_data([x2], [psf2(x2)]);\n",
    "\n",
    "    ax2.fill_between(x2, psf2(x2), 0, where=x2<s_rc.value, facecolor='green', alpha=0.5)\n",
    "    \n",
    "    logl_r = lambdify(r, diff(ps,r)/ps, \"numpy\")\n",
    "    logr_r = lambdify(r, diff(ho,r)/ho, \"numpy\")\n",
    "\n",
    "    logl = diff(ps, r).subs(r, rc).evalf()/ps.subs(r, rc).evalf()\n",
    "    logr = diff(ho, r).subs(r, rc).evalf()/ho.subs(r, rc).evalf()\n",
    "    \n",
    "        \n",
    "    line_logdr_rho.set_data([x3],[logr_r(x3)])\n",
    "    line_logdr_ps.set_data([x3],[logl_r(x3)])\n",
    "    \n",
    "    point3.set_data([s_rc.value], [logr_r(s_rc.value)])\n",
    "    ann_rc3.remove()\n",
    "    ann_rc3 = ax3.annotate(\"$R_c$\", xy=(s_rc.value + 1.0,logr_r(s_rc.value)) , fontsize=20)\n",
    "    # ann_rc3.set_position((s_rc.value + 1.0,logr_r(s_rc.value)))\n",
    "    # ann_rc3.set_visible(True)\n",
    "    #ax3.set_ylim([list(line_logdr_rho.get_data())[1].min()-0.04, list(line_logdr_rho.get_data())[1].max()+0.04])\n",
    "    ax3.set_ylim(logr_r(s_rc.value)-0.1,logr_r(s_rc.value)+0.1)\n",
    "\n",
    "    # ann_logl.set_text(\"Logarithmic deriv. $\\psi(r_c)^{PS}$:\" + str(\"{:.10f}\".format(logl)))\n",
    "    # ann_logr.set_text(\"Logarithmic deriv. $\\psi(r_c)$:\" + str(\"{:.10f}\".format(logr)))\n",
    "    # ann_logl.set_text(\"\")\n",
    "    # ann_logr.set_text(\"\")\n",
    "    \n",
    "    \n",
    "    #text_l0.value = str(ps.subs(r, rc).evalf())\n",
    "    #text_l1.value = str(diff(ps, r).subs(r, rc).evalf())\n",
    "    #text_l2.value = str(diff(ps, r, r).subs(r, rc).evalf())\n",
    "    \n",
    "    #text_r0.value = str(ho.subs(r, rc).evalf())\n",
    "    #text_r1.value = str(diff(ho, r).subs(r, rc).evalf())\n",
    "    #text_r2.value = str(diff(ho, r, r).subs(r, rc).evalf())\n",
    "    \n",
    "    #cof_0.value = str(coff[0])\n",
    "    #cof_2.value = str(coff[1])\n",
    "    #cof_3.value = str(coff[2])\n",
    "    #cof_4.value = str(coff[3])\n",
    "\n",
    "    \n",
    "def plot_log_derivatives(b,rc,l):\n",
    "    return\n",
    "    \n",
    "def plot_ps_potential(b, rc, l):\n",
    "    devs = compute_right_derivative(rc)\n",
    "    coff = solver_kernel(devs, rc, l, b)\n",
    "    \n",
    "    pf  = r**(l+1)*exp(coff[0] + coff[1]*r**2 + coff[2]*r**3 + coff[3]*r**4)\n",
    "    psf = Ea - l*(l+1)/(2*r*r) + 1/(2*pf)*diff(pf, r, r)\n",
    "    psfnl = lambdify(r, psf, \"numpy\")\n",
    "    psfnr = compute_potential(l)\n",
    "        \n",
    "    x1 = np.linspace(0.001, rc, 800);\n",
    "    x2 = np.linspace(rc, 50, 800);\n",
    "    line_psv.set_data([np.concatenate((x1,x2))], [np.concatenate((psfnl(x1),psfnr(x2)))]);\n",
    "    ax4.set_ylim([psfnl(x1).min(axis=0)-0.05, max(psfnl(x1).max(axis=0)+0.1, 0.05)])\n",
    "\n",
    "def compute_pseudopotential(c):\n",
    "    global compute\n",
    "    compute.disabled = True\n",
    "    old_description = compute.description\n",
    "    compute.description = \"Computing...\"\n",
    "    \n",
    "    try:\n",
    "        try:\n",
    "            b = newton(lambda x: diff_norms(x, s_rc.value, l), x0 = 0.00001, tol = 1e-10, maxiter=100)\n",
    "\n",
    "            if abs(diff_norms(b, s_rc.value, l)) > 0.001:\n",
    "                ann_norm1.set_text(\"No numerical solution found!\");\n",
    "                ann_norm2.set_text(\"Please change $R_c$!\");\n",
    "                # ann_logl.set_text(\"\");\n",
    "                # ang_logr.set_text(\"\");\n",
    "                return None            \n",
    "        except Exception:\n",
    "            ann_norm1.set_text(\"No numerical solution found!\");\n",
    "            ann_norm2.set_text(\"Please change $R_c$!\");\n",
    "            # ann_logl.set_text(\"\");\n",
    "            # ann_logr.set_text(\"\");\n",
    "            return None\n",
    "        update_plot()\n",
    "        plot_ps_wavefunction(b, s_rc.value, l)\n",
    "        plot_ps_potential(b, s_rc.value, l)\n",
    "        #plot_log_derivatives(b, s_rc.value, l)\n",
    "        norm1, norm2 = compute_norms(b, s_rc.value, l)\n",
    "        # ann_norm1.set_text(\"\")\n",
    "        # ann_norm2.set_text(\"\")\n",
    "        ann_norm1.set_text(\"Yellow area: \" + str(\"{:.10f}\".format(norm1)))\n",
    "        ann_norm2.set_text(\"Green area: \" + str(\"{:.10f}\".format(norm2)))\n",
    "    finally:\n",
    "        compute.disabled = False\n",
    "        compute.description = old_description\n",
    "        \n",
    "compute.on_click(compute_pseudopotential)\n",
    "compute_pseudopotential(\"init\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"height:1px;border:none;color:#cccccc;background-color:#cccccc;\" />\n",
    "\n",
    "## **Legend**\n",
    "(How to use the interactive visualization)\n",
    "\n",
    "### Controls\n",
    "\n",
    "Here, we consider the wavefunctions from the solution of the \n",
    "Schrödinger equation of one hydrogen atom. One can choose the \n",
    "state with quantum numbers <em>n</em> and <em>l</em> via the sliders.\n",
    "The position of the cutoff distance $R_c$ can be tuned with the slider. \n",
    "Click the \"Compute pseudopotentials\" button to calculate the pseudowavefunction\n",
    "and pseudopotential. \n",
    "\n",
    "### Interactive figure\n",
    "The top left panel shows the calculated radial \n",
    "pseudowavefunction times $r$: $r R^{PS}_{nl}(r)$ (blue line) and the radial wavefunction \n",
    "times r: $rR_{nl}(r)$ (red line). The top right subplot shows their respective square moduli. To fulfill the norm-conserving\n",
    "condition, the green and yellow areas illustrated in this subplot should be equal. The square of the area\n",
    "should be 1 when $r$ goes to $+\\infty$. The bottom left panel shows plots of the logarithmic derivatives of the all-electron (in red) and pseudowavefunctions (in blue) respectively. The bottom right panel shows the calculated\n",
    "pseudopotential (in blue)."
   ]
  },
  {
   "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.7"
  },
  "voila": {
   "authors": "Dou Du and Giovanni Pizzi"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}