← Go back to the interactive notebook

Background Theory: Numerical Solution of 1D Time Dependent Schrödinger Equation for nuclear evolution on multiple electronic potential energy surfaces by the Multiple Split Operator Fourier Transform (MSOFT) method

Original notebook: theory_msoft.ipynb

The MSOFT algorithm is one approach to exactly (to within numerical accuracy) simulate the nonadiabatic dynamical evolution of quantum systems. Nonadiabatic effects must be accounted for when a given quantum mechanical system possesses more than one electronic state and the coupling between the different states is sufficiently large so as to induce transitions between them. For simplicity, in this notebook we shall consider the application of MSOFT to a two level system. It is, however, entirely possible (at least in theory) to use the method in the presence of arbitrarily many electronic states.

Propagation Scheme

Working in the adiabatic basis {∣R⟩∣ϕα(R)⟩}\{ |R \rangle | \phi_{\alpha} (R)\rangle \} (within which the electronic Hamiltonian h^el(R^)\hat{h}_{\text{el}} ( \hat{R}) is diagonal), the dynamical equations governing the evolution of the nuclear wavepacket components, c1(t)c_1(t) and c2(t)c_2(t), on each of the two electronic states are

iℏ∂∂tc1(R,t)={−ℏ22M∂2∂R2+E1(R)}c1(R,t)+ℏMD12(R)∂∂Rc2(R,t)i \hbar \frac{\partial }{\partial t} c_1(R,t) = \left\{ - \frac{\hbar^2}{2M} \frac{\partial^2 }{\partial R^2 } + E_1(R) \right\} c_1(R,t) + \frac{\hbar}{M} D _{12}(R) \frac{\partial }{\partial R} c_2(R,t)

and

iℏ∂∂tc2(R,t)={−ℏ22M∂2∂R2+E2(R)}c2(R,t)+ℏMD21(R)∂∂Rc1(R,t)i \hbar \frac{\partial }{\partial t} c_2(R,t) = \left\{ - \frac{\hbar^2}{2M} \frac{\partial^2 }{\partial R^2 } + E_2(R) \right\} c_2(R,t) + \frac{\hbar}{M} D _{21}(R) \frac{\partial }{\partial R} c_1(R,t)

where D12(R)=D21∗(R)=⟨ϕ1(R)∣ddRϕ2(R)⟩D _{12}(R) = D _{21}^*(R) = \langle \phi _1(R) | \frac{d}{dR} \phi_2(R) \rangle.

We may also consider introducing the so-called diabatic basis {∣R⟩∣λ⟩}\{| R \rangle | \lambda \rangle\} in which we can write the state of the system as ∣Ψ(t)⟩=∑λ∫dRχλ(R,t)∣R⟩∣λ⟩|\Psi (t) \rangle = \sum \limits_{\lambda}^{} \int dR \chi _{\lambda}(R,t) |R\rangle | \lambda \rangle and for which the kinetic energy operator is now diagonal (but the electronic Hamiltonian is non-diagonal). In this basis the above evolution equations take the form

iℏ∂∂tχ1(R,t)={−ℏ22M∂2∂R2+h11(R)}χ1(R,t)+h12(R)χ2(R,t)i \hbar \frac{\partial }{\partial t} \chi_1(R,t) = \left\{ - \frac{\hbar^2}{2M} \frac{\partial^2 }{\partial R^2 } + h_{11}(R) \right\} \chi_1(R,t) + h _{12}(R) \chi_2(R,t) iℏ∂∂tχ2(R,t)={−ℏ22M∂2∂R2+h22(R)}χ2(R,t)+h21(R)χ1(R,t)i \hbar \frac{\partial }{\partial t} \chi_2(R,t) = \left\{ - \frac{\hbar^2}{2M} \frac{\partial^2 }{\partial R^2 } + h_{22}(R) \right\} \chi_2(R,t) + h _{21}(R) \chi_1(R,t)

Steps of the MSOFT algorithm

Starting with the initial state ∣Ψ(R,0)⟩|\Psi(R,0)\rangle of the system and the electronic Hamiltonian h^el\hat{h}_{\text{el}} expressed in the diabatic basis the MSOFT algorithm proceeds as follows:

Perform the following loop for a specifies number of timesteps Nsteps\text{N} _{ \text{steps} }:

  1. Transform to the adiabatic basis: ∣Ψ(R,0)⟩Di→∣Ψ(R,0)⟩Adi=De^∣Ψ(R,0)⟩Di|\Psi(R,0)\rangle _{ \text{Di} } \to |\Psi(R,0)\rangle _{ \text{Adi} } = \hat{D_e} |\Psi(R,0)\rangle _{ \text{Di} }
  2. Apply half-step potential propagator: ∣Φ(R,dt2)⟩Adi=exp⁡(−iℏdt2E^)∣Ψ(R,0)⟩Adi|\Phi(R, \frac{dt}{2} )\rangle _{ \text{Adi} }= \exp(- \frac{i}{\hbar} \frac{dt}{2} \hat{E} ) |\Psi(R,0)\rangle _{ \text{Adi} }.
  3. Transform back to the diabatic basis ∣Φ(R,dt2)⟩Adi→∣Φ(R,dt2)⟩Di=De^T∣Φ(R,dt2)⟩Adi|\Phi(R, \frac{dt}{2} )\rangle _{ \text{Adi} } \to |\Phi(R, \frac{dt}{2} )\rangle _{ \text{Di} } = \hat{D_e}^T |\Phi(R, \frac{dt}{2} )\rangle _{ \text{Adi} }
  4. Apply fast Fourier transform (FFT) to transform state to reciprocal space and render kinetic energy operator diagonal in nuclear part of total tensor product space HN⊗He\mathcal{H}_N \otimes \mathcal{H} _e: ∣Φ(R,dt2)⟩Di→FFT[∣Φ(R,dt2)⟩Di]=∣Φ~(P,dt2)⟩Di|\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} }
  5. Apply full-timestep kinetic energy propagator: ∣Ξ~(P,dt)⟩Di=exp⁡(−iℏP^22Mdt)∣Φ~(P,dt)⟩Di|\tilde{\Xi}(P, dt )\rangle _{ \text{Di} }= \exp(- \frac{i}{\hbar} \frac{\hat{P}^2}{2M} dt ) |\tilde{\Phi}(P,dt)\rangle _{ \text{Di} }.
  6. Apply an inverse FFT to transform state back to direct space: ∣Ξ~(P,dt)⟩Di→IFFT[∣Ξ~(P,dt)⟩Di]=∣Ξ(R,dt)⟩Di|\tilde{\Xi}(P, dt )\rangle _{ \text{Di} } \to \text{IFFT}[|\tilde{\Xi}(P, dt )\rangle _{ \text{Di}}] = |\Xi(R, dt )\rangle _{ \text{Di}}.
  7. Transform to the adiabatic basis: ∣Ξ(R,dt)⟩Di→∣Ξ(R,dt)⟩Adi=De^T∣Ξ(R,dt)⟩Di∣|\Xi(R, dt )\rangle _{ \text{Di}} \to |\Xi(R, dt )\rangle _{ \text{Adi}} = \hat{D_e}^T |\Xi(R, dt )\rangle _{ \text{Di}}|.
  8. Apply half-step potential propagator: ∣Ψ(R,dt)⟩Adi=exp⁡(−iℏdt2E^)∣Ξ(R,dt)⟩Adi|\Psi(R, dt)\rangle _{ \text{Adi} }= \exp(- \frac{i}{\hbar} \frac{dt}{2} \hat{E} ) |\Xi(R,dt)\rangle _{ \text{Adi} }.

This leaves us with the state of the system evolved by one full timestep dtdt. We then repeat steps (2)-(8) Nsteps−1\text{N} _{steps}-1 times to achieve a total propagation time of t=Nsteps⋅dtt= \text{N} _{steps} \cdot dt for the system.