Numerical Techniques

This page contains brief descriptions of various numerical and computational topics related to my project.

Time Propagation

Analytic solutions of the Schrodinger equation only exist for simple, usually highly symmetric, situations. For realistic problems, it is only feasible to solve the Schrodinger equation numerically. Instead of obtaining an expression for \(\Psi\left(\vec{r}, t\right)\) where \(\vec{r}, t\) are continuous variables, most numerical methods iteratively approximate \(\Psi\) at discrete \(t\) (and \(\vec{r}\)). given the initial state of a system.

In atomic units, the time-dependent Schrodinger equation is:

\[\frac{\partial}{\partial t} \Psi \left( \vec{r}, t \right) = -i \hat{H} \Psi \left(\vec{r}, t \right)\]

In general, the process of integrating the time-dependent Schrodinger equation can be numerically approximated over a small time step by:

\[\Psi \left(\vec{r}, t + \Delta t\right) = e^{-i \hat{H}\left(t\right)\Delta t} \Psi \left( \vec{r}, t \right)\]

Since directly calculating a (matrix) exponential of the Hamiltonian at each time step is normally computationally unfeasible, various time propagation methods indirectly calculate or approximate the exponential in a variety of ways, including grid-based discretization and expansion in terms of various basis functions. The method discussed here uses the Crank-Nicholson scheme, which is a finite difference method on a numerical grid.

Cayley Approximation and Crank-Nicholson Method

It can be shown that:

\[\Psi\left(\vec{r}, t + \Delta t\right) = \frac {1 - i \hat{H}\left(t + \frac{\Delta t}{2}\right)\frac{\Delta t}{2}} {1 + i \hat{H}\left(t + \frac{\Delta t}{2}\right)\frac{\Delta t}{2}} \Psi\left(\vec{r},t\right) + \mathcal{O}\left(\Delta t^3 \right)\]

which provides an accurate scheme for numerically approximating the time propagation of an arbitrary wavefunction given sufficiently small time steps. This approximation is known as the Cayley Approximation.

If the Cayley approximation is used in conjunction with a discrete spatial numerical grid (which requires the use of numerical differentiation), the resulting scheme is known as the Crank-Nicholson approximation. For each propagated time step, using numerical derivatives for the Hamiltonian discretized over the numerical grid results in a tridiagonal (for a 3-point second derivative) system of linear equations.

In our case, we used a internuclear numerical grid 80 a.u. in size with spatial grid steps 0.01 a.u. and time steps on the order of 10 fs. Numerical convergence was manually verified across different grid sizes and time steps.

Numerical Optimizations

As described above, the Crank-Nicholson scheme on a numerical grid with \(n\) points involves solving a general system of linear equations of the form \(A \Psi = B\), where \(A\) and \(B\) are \(n \times n\) matrices. A numerically optimized representation of the Cayley Approximation exists which has algorithmic complexity on the order of solving a single system of linear equations (instead of solving \(n\) different systems) for each time step. This form is:

\[\Psi\left(t + \Delta t\right) = \Phi\left(t\right) - \Psi\left(t\right)\]

where an auxillary wave function \(\Phi\) is the solution to the following system of equations:

\[\Psi\left(t\right) = \frac{1}{2} \left[1 + i \hat{H}\frac{\Delta t}{2} \right] \Phi\left(t\right)\]