The Multiple Split Operator Fourier Transform (MSOFT) method¶

Authors: Taylor Baird and Sara Bonella

Go back to index

Source code: https://github.com/osscar-org/quantum-mechanics/blob/master/notebook/quantum-mechanics/msoft.ipynb

This notebook carries out the numerical solution of the 1D time-dependent Schrödinger equation for nuclear evolution on multiple electronic potential energy surfaces via the Multiple Split Operator Fourier Transform (MSOFT) method.

Goals¶

• Understand the different steps in the MSOFT algorithm and how they are translated to code.
• Familiarize yourself with the key assumptions underlying the algorithm and how their validity depends on the values of the parameters used in the simulation.
• Use the MSOFT algorithm to investigate the nonadiabatic dynamics of a nuclear wavepacket propagating in two different 2-level systems (a double harmonic potential and the Tully potential).

Background Theory¶

More in the background theory

Tasks and exercises¶

1. Investigate the dependence on timestep of the stability of the dynamics (try moving the slider for $dt$ and monitor the behaviour of the various control properties of the simulation).

   Solution

   One may observe that as the timestep employed in the simulation is increased, the conservation of total energy of the system degrades until eventually the timestep is so large that the dynamics becomes totally unstable. The reason why this integration scheme does not conserve total energy may be attributed to the non-commutativity of the split-operator propagator with the Hamiltonian. It is worth noting, however, that norm conservation is maintained even as one increases the timestep. This latter fact is due to the unitarity of the propagator in the split-operator scheme.

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).

   Solution

   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$).

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)?

   Solution

   In order to implement the MSOFT algorithm on a computer it is necessary to discretize the nuclear wavefunction. To do so, we must introduce a grid of points that make up the space throughout which the wavefunction extends. Say that for each dimension of the grid we use $N$ grid points. For a system in $d$ dimensions, this means that we shall require $N^d$ points to represent the wavefunction of a single nucleus. Now, if we want to instead consider a system of say $n$ nuclei then we find that a total number of $N^{nd}$ grid points are required. In other words, the amount of data required to represent our system scales exponentially with increasing system size. Moreover, since we must loop over each of these grid points to carry out the time evolution of our system - the number of operations required also scales exponentially. This is what renders the MSOFT algorithm unsuitable for large systems.

4. Investigate the effect of increasing the coupling strength between different electronic states on the evolution of the two state populations (vary the $C$ slider for each potential).

   Solution

   One can note that by increasing the value of $C$ the probability of inducing a transition of the nuclear wavepacket from one electronic potential energy surface to the other is modulated. You can trace this back to the equations of motion given in the theory section.

Legend (How to use the interactive visualization)¶

Interactive figures¶

The interactive figures consist of 4 subplots (clockwise from the top):

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).

2. Plot of the kinetic, potential and total energies vs simulation timestep.

3. Evolution of the populations on the two different electronic PESs.

4. Evolution of the total norm of the nuclear wavepacket.

Below the plots there is also a display field indicating the current value of the norm of the total wavefunction.

Controls¶

• 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.

• One can choose on which electronic state the Gaussian nuclear wavepacket is initialized by using the "Choose initial electronic state" dropdown.

• The parameters of each potential may be varied using the various sliders.

• To update the potential, the "Update potential" button should be pressed, followed by the "Update parameters" button in the "Set simulation parameters" drawer.

• The various simulation parameters (nuclear mass, timestep (dt), and initial momentum of the nuclear wavepacket (p0)) may be set using the sliders in the "Set simulation parameters" drawer.