diff --git a/notebook/quantum-mechanics/msoft.ipynb b/notebook/quantum-mechanics/msoft.ipynb index 86d041e..72977ab 100644 --- 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", "
\n", " Solution\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", "
\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,6 +527,31 @@ " 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", @@ -518,7 +559,7 @@ "\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,12 +840,11 @@ " 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", @@ -782,10 +852,10 @@ " 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 index d720b12..b73e3b9 100644 --- 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} }$
\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} }$.
\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}}$.
\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}}|$.
\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 $.
\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} }$.
\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", "

\n" @@ -74,7 +74,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.8" + "version": "3.10.6" } }, "nbformat": 4,