Revised MSOFT notebook. Some modifications and corrections to text. M…

…ade it possible to update plots OTF rather than have static images.
osscar-org · May 21, 2023 · 81e425c · 81e425c
1 parent 0ba9cbb
commit 81e425c
Show file tree

Hide file tree

Showing 2 changed files with 93 additions and 23 deletions.
diff --git a/notebook/quantum-mechanics/msoft.ipynb b/notebook/quantum-mechanics/msoft.ipynb
@@ -56,8 +56,7 @@
     "2. What dictates the maximum size of the timestep that we can use for the MSOFT algorithm (consider the main assumptions/approximations that are made when formulating the propagation scheme for MSOFT).\n",
     "    <details>\n",
     "    <summary style=\"color: red\">Solution</summary>\n",
-    "    Recall that the central approximation used in the derivation of the MSOFT propagation scheme is the Trotter product formula: $e^{A+B}  \\underbrace{=}_{P \\to \\infty} (e^{\\frac{A}{P}} e^{\\frac{B}{P}})^{P}$. This approximation becomes more and more accurate the larger we make the value of $P$. In our specific case we have $e^{\\frac{it}{\\hbar} (\\hat{T} + \\hat{V} )}$ and we approximate this via the product \n",
-    "        $(e^{\\frac{idt}{\\hbar}\\hat{T} } e^{\\frac{idt}{\\hbar}\\hat{V} })^{N_{\\text{steps}}}$ where $N_{\\text{steps}}\\cdot dt = t$.\n",
+    "    Recall that the central approximation used in the derivation of the MSOFT propagation scheme is in the truncation of the (symmetric) Trotter product formula: $e^{A+B} = \\lim\\limits_{P \\to \\infty} (e^{\\frac{B}{2P}}e^{\\frac{A}{P}} e^{\\frac{B}{2P}})^{P} \\approx (e^{\\frac{B}{2N}}e^{\\frac{A}{N}} e^{\\frac{B}{2N}})^{N}$ for $N$ sufficiently large. This approximation becomes more and more accurate the larger we make the value of $N$. In our specific case we have $e^{\\frac{it}{\\hbar} (\\hat{T} + \\hat{V} )}$ and we approximate this via the product $(e^{\\frac{it}{2\\hbar N_{\\text{steps}}}\\hat{V} }e^{\\frac{it}{\\hbar N_{\\text{steps}}}\\hat{T} } e^{\\frac{it}{2\\hbar N_{\\text{steps}}}\\hat{V} })^{N_{\\text{steps}}}$ where $N_{\\text{steps}}\\cdot dt = t$. This approximation therefore becomes increasingly more accurate the larger $N_{\\text{steps}}$ (or equivalently the smaller we make $dt$).\n",
     "    </details>\n",
     "\n",
     "3. Why is the use of the MSOFT algorithm unfeasible for larger systems (think about how much data is needed to represent the state of the system on a computer and how many operations are required to carry out the propagation)?\n",
@@ -448,9 +447,25 @@
     "    height=700,\n",
     ")\n",
     "\n",
-    "l = Layout(flex='0 1 auto', height='80px', min_height='80px', width='auto')\n",
-    "pot_eqn=Label(value=r'\\(h(x_i)=\\begin{pmatrix} A(1-e^{-B|\\frac{LX}{2}-x_i|})\\cdot \\text{sgn}(\\frac{LX}{2}-x_i) & Ce^{-D(\\frac{LX}{2}-x_i)^2} \\\\Ce^{-D(\\frac{LX}{2}-x_i)^2} & -A(1-e^{-B|\\frac{LX}{2}-x_i|})\\cdot \\text{sgn}(\\frac{LX}{2}-x_i)\\end{pmatrix}\\)',layout=l)\n",
-    "#pot_eqn=Label(value=r'\\(h(x_i)=\\begin{pmatrix} h(x_i)=A(x_i-B)^2 & Ce^{(-D(x_i-\\frac{LX}{2})^2)} \\\\Ce^{(-D(x_i-\\frac{LX}{2})^2)} & A(x_i-b)^2\\end{pmatrix}\\)',layout=l)\n",
+    "    \n",
+    "\n",
+    "pot_fig = plt.figure()\n",
+    "pot_ax=pot_fig.add_subplot()\n",
+    "pot_ax.set_ylim(-0.05,0.25)\n",
+    "pot_ax.set_xlim(0,50)\n",
+    "\n",
+    "\n",
+    "pot_prev_V1_ad_x_line, = pot_ax.plot([], [], c='purple', label=r'$V_{ad}(x)$')\n",
+    "pot_prev_V2_ad_x_line, = pot_ax.plot([], [], c='cyan')#, label=r'$V_{ad}(x)$')\n",
+    "pot_prev_V1_di_x_line, = pot_ax.plot([], [], c='green', label=r'$V_{di}(x)$')\n",
+    "pot_prev_V2_di_x_line, = pot_ax.plot([], [], c='green')#, label=r'$V_{di}(x)$')\n",
+    "pot_prev_V2_di_coupling_x_line, = pot_ax.plot([], [], c='red', label='$h_{12}(x)$')#, label=r'$V_{di}(x)$')\n",
+    "\n",
+    "pot_leg= pot_ax.legend(prop=dict(size=12))\n",
+    "\n",
+    "l = Layout(flex='0 1 auto', height='80px', min_height='40px', width='auto')\n",
+    "#pot_eqn=Label(value=r'\\(h(x_i)=\\begin{pmatrix} A(1-e^{-B|\\frac{LX}{2}-x_i|})\\cdot \\text{sgn}(\\frac{LX}{2}-x_i) & Ce^{-D(\\frac{LX}{2}-x_i)^2} \\\\Ce^{-D(\\frac{LX}{2}-x_i)^2} & -A(1-e^{-B|\\frac{LX}{2}-x_i|})\\cdot \\text{sgn}(\\frac{LX}{2}-x_i)\\end{pmatrix}\\)',layout=l)\n",
+    "pot_eqn=Label(value=r'\\(h(x_i)=\\begin{pmatrix} &\\color{green} {A(x_i-B)^2} & \\color{red}{Ce^{(-D(x_i-\\frac{LX}{2})^2)}} \\\\&\\color{red}{Ce^{(-D(x_i-\\frac{LX}{2})^2)}} & \\color{green}{A(x_i-b)^2}\\end{pmatrix}\\)',layout=l)\n",
     "\n",
     "\n",
     "pot_img.layout = layout_visible\n",
@@ -468,9 +483,9 @@
     "        s_tully_B.layout = layout_hidden \n",
     "        s_tully_C.layout = layout_hidden \n",
     "        s_tully_D.layout = layout_hidden \n",
-    "        pot_img.layout = layout_visible\n",
+    "        #pot_img.layout = layout_visible\n",
     "        pot_img.value = image1\n",
-    "        pot_eqn.value=r'\\(h(x_i)=\\begin{pmatrix} A(1-e^{-B|\\frac{LX}{2}-x_i|})\\cdot \\text{sgn}(\\frac{LX}{2}-x_i) & Ce^{-D(\\frac{LX}{2}-x_i)^2} \\\\Ce^{-D(\\frac{LX}{2}-x_i)^2} & -A(1-e^{-B|\\frac{LX}{2}-x_i|})\\cdot \\text{sgn}(\\frac{LX}{2}-x_i)\\end{pmatrix}\\)'\n",
+    "        pot_eqn.value = r'\\(h(x_i)=\\begin{pmatrix} &\\color{green} {A(x_i-B)^2} & \\color{red}{Ce^{(-D(x_i-\\frac{LX}{2})^2)}} \\\\&\\color{red}{Ce^{(-D(x_i-\\frac{LX}{2})^2)}} & \\color{green}{A(x_i-b)^2}\\end{pmatrix}\\)'\n",
     "\n",
     "        \n",
     "    elif pot_select.index == 1:\n",
@@ -484,9 +499,10 @@
     "        s_tully_B.layout = layout_visible\n",
     "        s_tully_C.layout = layout_visible\n",
     "        s_tully_D.layout = layout_visible\n",
-    "        pot_img.layout = layout_visible\n",
+    "        # pot_img.layout = layout_visible\n",
     "        pot_img.value = image2         \n",
-    "        pot_eqn.value = r'\\(h(x_i)=\\begin{pmatrix} A(x_i-B)^2 & Ce^{(-D(x_i-\\frac{LX}{2})^2)} \\\\Ce^{(-D(x_i-\\frac{LX}{2})^2)} & A(x_i-b)^2\\end{pmatrix}\\)'\n",
+    "        pot_eqn.value=r'\\(h(x_i)=\\begin{pmatrix} \\color{green}{A(1-e^{-B|\\frac{LX}{2}-x_i|})\\cdot \\text{sgn}(\\frac{LX}{2}-x_i)} & \\color{red}{Ce^{-D(\\frac{LX}{2}-x_i)^2}} \\\\\\color{red}{Ce^{-D(\\frac{LX}{2}-x_i)^2}} & \\color{green}{-A(1-e^{-B|\\frac{LX}{2}-x_i|})\\cdot \\text{sgn}(\\frac{LX}{2}-x_i)}\\end{pmatrix}\\)'\n",
+    "\n",
     "        \n",
     "        \n",
     "def state_change(change):\n",
@@ -500,7 +516,7 @@
     "pot_select.observe(pot_change, names='value', type='change');\n",
     "\n",
     "def on_pot_update(event):\n",
-    "    global h_x\n",
+    "    global h_x, S\n",
     "    if pot_select.index == 0:\n",
     "        h_x = double_harmonic(x, s_dharm_A.value,s_dharm_B.value,s_dharm_b.value,s_dharm_C.value,s_dharm_D.value)\n",
     "    else:\n",
@@ -511,14 +527,39 @@
     "        ax2.set_xlim(0,tot_time)\n",
     "        ax3.set_xlim(0,tot_time)\n",
     "        ax4.set_xlim(0,tot_time)\n",
+    "        \n",
+    "    # global S,h_x,p0,state_idx,pot_leg   \n",
+    "    # # two components (one for each electronic state)\n",
+    "    # psi_x0 = np.zeros((N,2), dtype=np.complex128) \n",
+    "    # # default IC is a Gaussian wavepacket on 1st state and 0 component on second state but \n",
+    "    # # user can choose the reverse of this.\n",
+    "    # psi_x0[:,state_idx] = gauss_x(x, s0, x0, p0) \n",
+    "    S = MSOFT(x = x, dt = dt, psi_0=psi_x0, h=h_x, hbar=hbar, m=m)\n",
+    "    \n",
+    "    # pot_ax.clear()\n",
+    "    # pot_leg.remove()\n",
+    "\n",
+    "    # pot_ax.get_legend().remove()\n",
+    "    # pot_prev_V1_ad_x_line, = pot_ax.plot(S.x, S.E_eigvals[:,0], c='purple', label=r'$V_{ad}(x)$')\n",
+    "    # pot_prev_V2_ad_x_line, = pot_ax.plot(S.x, S.E_eigvals[:,1], c='cyan')#, label=r'$V_{ad}(x)$')\n",
+    "    # pot_prev_V1_di_x_line, = pot_ax.plot(S.x, S.h[:,0,0], c='green', label=r'$V_{di}(x)$')\n",
+    "    # pot_prev_V2_di_x_line, = pot_ax.plot(S.x, S.h[:,1,1], c='green')#, label=r'$V_{di}(x)$')\n",
+    "    # pot_prev_V2_di_coupling_x_line, = pot_ax.plot(S.x,S.h[:1,0], c='red', label='$h_{12}(x)$')#, label=r'$V_{di}(x)$')\n",
+    "    # pot_ax.legend()\n",
+    "    pot_prev_V1_ad_x_line.set_data(S.x, S.E_eigvals[:,0])#, c='purple', label=r'$V_{ad}(x)$')\n",
+    "    pot_prev_V2_ad_x_line.set_data(S.x, S.E_eigvals[:,1])#,c='purple')\n",
+    "    pot_prev_V1_di_x_line.set_data(S.x, S.h[:,0,0])#, c='green', label=r'$V_{di}(x)$')\n",
+    "    pot_prev_V2_di_x_line.set_data(S.x, S.h[:,1,1])#,c ='green')\n",
+    "    #pot_prev_V2_di_coupling_x_line.set_data(S.x,0.1*np.ones(len(S.h)))#,  c='red', label='$h_{12}(x)$')\n",
+    "    S.evolution(1)\n",
     "\n",
     "\n",
     "pot_update = Button(description=\"Update potential\")\n",
     "pot_update.on_click(on_pot_update)\n",
     "\n",
     "# The accordion for the potential\n",
     "label1 = Label(\"(Press update potential and then update parameters to modify the potential.)\")\n",
-    "pot_accordion = Accordion(children=[VBox([pot_select, pot_img, pot_eqn, state_select,\n",
+    "pot_accordion = Accordion(children=[VBox([pot_select, pot_eqn, state_select,\n",
     "                   HBox([s_dharm_A,s_dharm_B,s_dharm_b,s_dharm_C,s_dharm_D]),\n",
     "                   HBox([s_tully_A,s_tully_B,s_tully_C,s_tully_D]),\n",
     "                   HBox([pot_update, label1])])], selected_index = None)\n",
@@ -551,16 +592,32 @@
     "    p0 = w_p0.value;\n",
     "    \n",
     "def on_init_change(b):\n",
-    "    global S,h_x,p0,state_idx    \n",
+    "    global S,h_x,p0,state_idx,pot_leg   \n",
     "    # two components (one for each electronic state)\n",
     "    psi_x0 = np.zeros((N,2), dtype=np.complex128) \n",
     "    # default IC is a Gaussian wavepacket on 1st state and 0 component on second state but \n",
     "    # user can choose the reverse of this.\n",
     "    psi_x0[:,state_idx] = gauss_x(x, s0, x0, p0) \n",
     "    S = MSOFT(x = x, dt = dt, psi_0=psi_x0, h=h_x, hbar=hbar, m=m)\n",
-    "    S.evolution(1)\n",
     "    \n",
-    "    w_mass.observe(on_mass_change, names='value')\n",
+    "    # pot_ax.clear()\n",
+    "    # pot_leg.remove()\n",
+    "\n",
+    "    # pot_ax.get_legend().remove()\n",
+    "    # pot_prev_V1_ad_x_line, = pot_ax.plot(S.x, S.E_eigvals[:,0], c='purple', label=r'$V_{ad}(x)$')\n",
+    "    # pot_prev_V2_ad_x_line, = pot_ax.plot(S.x, S.E_eigvals[:,1], c='cyan')#, label=r'$V_{ad}(x)$')\n",
+    "    # pot_prev_V1_di_x_line, = pot_ax.plot(S.x, S.h[:,0,0], c='green', label=r'$V_{di}(x)$')\n",
+    "    # pot_prev_V2_di_x_line, = pot_ax.plot(S.x, S.h[:,1,1], c='green')#, label=r'$V_{di}(x)$')\n",
+    "    # pot_prev_V2_di_coupling_x_line, = pot_ax.plot(S.x,S.h[:1,0], c='red', label='$h_{12}(x)$')#, label=r'$V_{di}(x)$')\n",
+    "    # pot_ax.legend()\n",
+    "    pot_prev_V1_ad_x_line.set_data(S.x, S.E_eigvals[:,0])#, c='purple', label=r'$V_{ad}(x)$')\n",
+    "    pot_prev_V2_ad_x_line.set_data(S.x, S.E_eigvals[:,1])#,c='purple')\n",
+    "    pot_prev_V1_di_x_line.set_data(S.x, S.h[:,0,0])#, c='green', label=r'$V_{di}(x)$')\n",
+    "    pot_prev_V2_di_x_line.set_data(S.x, S.h[:,1,1])#,c ='green')\n",
+    "    #pot_prev_V2_di_coupling_x_line.set_data(S.x,0.1*np.ones(len(S.h)))#,  c='red', label='$h_{12}(x)$')\n",
+    "    S.evolution(1)\n",
+    "\n",
+    "w_mass.observe(on_mass_change, names='value')\n",
     "w_dt.observe(on_dt_change, names='value')\n",
     "w_p0.observe(on_p0_change, names='value')\n",
     "w_init = Button(description=\"Update parameters\");\n",
@@ -625,6 +682,8 @@
     "V1_di_x_line, = ax1.plot([], [], c='green', label=r'$V_{di}(x)$')\n",
     "V2_di_x_line, = ax1.plot([], [], c='green')#, label=r'$V_{di}(x)$')\n",
     "\n",
+    "\n",
+    "\n",
     "psi1_x_line.set_visible(True)\n",
     "psi2_x_line.set_visible(True)\n",
     "\n",
@@ -660,6 +719,13 @@
     "V1_di_x_line.set_data(S.x, S.h[:,0,0])\n",
     "V2_di_x_line.set_data(S.x, S.h[:,1,1])\n",
     "\n",
+    "pot_prev_V1_ad_x_line.set_data(S.x, S.E_eigvals[:,0])\n",
+    "pot_prev_V2_ad_x_line.set_data(S.x, S.E_eigvals[:,1])\n",
+    "pot_prev_V1_di_x_line.set_data(S.x, S.h[:,0,0])\n",
+    "pot_prev_V2_di_x_line.set_data(S.x, S.h[:,1,1])\n",
+    "pot_prev_V2_di_coupling_x_line.set_data(S.x, S.h[:,1,0])\n",
+    "\n",
+    "\n",
     "psi1_x_line.set_data(S.x, 1. * abs(S.psi_x[:,0])+S.E_eigvals[:,0])\n",
     "psi2_x_line.set_data(S.x, 1. * abs(S.psi_x[:,1])+S.E_eigvals[:,1])\n",
     "\n",
@@ -696,6 +762,11 @@
     "pause = True\n",
     "\n",
     "def init():\n",
+    "    # pot_prev_V1_ad_x_line.set_data([],[])\n",
+    "    # pot_prev_V2_ad_x_line.set_data([],[])\n",
+    "    # pot_prev_V1_di_x_line.set_data([],[])\n",
+    "    # pot_prev_V2_di_x_line.set_data([],[])\n",
+    "    # pot_prev_V2_di_coupling_x_line.set_data([],[])\n",
     "    psi1_x_line.set_data([], [])\n",
     "    psi2_x_line.set_data([], [])\n",
     "    V1_ad_x_line.set_data([],[])\n",
@@ -769,23 +840,22 @@
     "        pause=True\n",
     "        button_pause.description = \"Play\"\n",
     "        anim.event_source.stop()\n",
-    "       \n",
+    "        \n",
     "def onClick(event):\n",
     "    global pause\n",
     "    global init_time\n",
     "    init_time=datetime.now()\n",
-    "    \n",
     "    pause ^= True\n",
     "    if button_pause.description == \"Pause\":\n",
     "        button_pause.description = \"Play\"\n",
     "        anim.event_source.stop()\n",
     "    else:\n",
     "        button_pause.description = \"Pause\"\n",
     "        anim.event_source.start()\n",
-    "\n",
     "        \n",
     "init_time=datetime.now()\n",
     "\n",
+    "\n",
     "anim = animation.FuncAnimation(fig, animate, init_func=init, frames=frames, interval=1, blit=True)\n",
     "button_pause = Button(description=\"Play\");\n",
     "button_pause.on_click(onClick)\n",
@@ -803,9 +873,9 @@
     "\n",
     "# Legend (How to use the interactive visualization)\n",
     "\n",
-    "## Visualization\n",
+    "## Interactive figures\n",
     "\n",
-    "The visualization consists of 4 subplots (clockwise from the top):\n",
+    "The interactive figures consist of 4 subplots (clockwise from the top):\n",
     "1. The main animation pane. This shows the evolution of the two components of the nuclear wavefunction (red and magenta curves) on their respective electronic potential energy surfaces (PESs). The potential energy surfaces are plotted in both the adiabatic (purple) and diabatic bases (green).\n",
     "\n",
     "2. Plot of the kinetic, potential and total energies vs simulation timestep.\n",
@@ -816,7 +886,7 @@
     "\n",
     "Below the plots there is also a display field indicating the current value of the norm of the total wavefunction.\n",
     "\n",
-    "## Interactivity controls\n",
+    "## Controls\n",
     "\n",
     "* The type of potential may be selected by first expanding the \"Select potential and set parameters\" drawer. The two options are: a double harmonic potential and a Tully model potential.\n",
     "* One can choose on which electronic state the Gaussian nuclear wavepacket is initialized by using the \"Choose initial electronic state\" dropdown.\n",
@@ -844,7 +914,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.12"
+   "version": "3.10.6"
   }
  },
  "nbformat": 4,

diff --git a/notebook/quantum-mechanics/theory/theory_msoft.ipynb b/notebook/quantum-mechanics/theory/theory_msoft.ipynb
@@ -51,7 +51,7 @@
     " $\\mathcal{H}_N \\otimes \\mathcal{H} _e$: $|\\Phi(R, \\frac{dt}{2} )\\rangle _{ \\text{Di} } \\to \\text{FFT}[ |\\Phi(R, \\frac{dt}{2} )\\rangle _{ \\text{Di} }] =  |\\tilde{\\Phi}(P, \\frac{dt}{2} )\\rangle _{ \\text{Di} }$<br/>\n",
     " $\\qquad$(5) Apply full-timestep kinetic energy propagator: $|\\tilde{\\Xi}(P, dt )\\rangle _{ \\text{Di} }=   \\exp(- \\frac{i}{\\hbar} \\frac{\\hat{P}^2}{2M} dt ) |\\tilde{\\Phi}(P,dt)\\rangle _{ \\text{Di} }$.<br/>\n",
     "$\\qquad$(6) Apply an inverse FFT to transform state back to direct space: $|\\tilde{\\Xi}(P, dt )\\rangle _{ \\text{Di} } \\to \\text{IFFT}[|\\tilde{\\Xi}(P, dt )\\rangle _{ \\text{Di}}] = |\\Xi(R, dt )\\rangle _{ \\text{Di}}$.<br/>\n",
-    " $\\qquad$(7) Transform to the adiabatic basis: $|\\Xi(R, dt )\\rangle _{ \\text{Di}} \\to |\\Xi(R, dt )\\rangle _{ \\text{Adi}} = \\hat{D_e}^T |\\Xi(R, dt )\\rangle _{ \\text{Di}}|$.<br/>\n",
+    " $\\qquad$(7) Transform to the adiabatic basis: $|\\Xi(R, dt )\\rangle _{ \\text{Di}} \\to |\\Xi(R, dt )\\rangle _{ \\text{Adi}} = \\hat{D_e}^T |\\Xi(R, dt )\\rangle _{ \\text{Di}}\\rangle $.<br/>\n",
     "$\\qquad$(8)  Apply half-step potential propagator:  $|\\Psi(R, dt)\\rangle _{ \\text{Adi} }= \\exp(- \\frac{i}{\\hbar} \\frac{dt}{2} \\hat{E} ) |\\Xi(R,dt)\\rangle _{ \\text{Adi} }$.<br/>\n",
     "$\\qquad$This leaves us with the state of the system evolved by one full timestep $dt$. We then repeat steps (2)-(8) $ \\text{N} _{steps}-1$ times to achieve a total propagation time of $t= \\text{N} _{steps} \\cdot dt$ for the system.\n",
     "</p>\n"
@@ -74,7 +74,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.8"
+   "version": "3.10.6"
   }
  },
  "nbformat": 4,