Numerical Propagation of Schrodinger's Equation

The fortran program (onedschreq.f) and an input data file (onedtimsch1.dat) can be used to numerically propagate Schrodinger's equation. In running this program, you can get a lot of output so make sure you only output the wave function as often as is necessary to see it evolve smoothly in time.

The fortran program time propagates the wave function for an electron in a time independent potential. It uses a number of simple approximations that will give OK results without much fuss. There are a number of ways that the program can be improved by using better approximations for derivatives; however, I feel that it is better to have a simple code that is easily understood instead of a more complicated code for students to use. The code uses atomic units for all parameters (hbar = 1, mass of electron = 1, charge of electron = -1). To convert to normal units, 1 atomic unit of distance is 0.529E-8 cm and 1 atomic unit of time is 2.42E-17 sec and 1 atomic unit of velocity is c/137 (c is the speed of light).

In atomic units the time dependent Scrhodinger Equation in 1 dimension is

i d psi(x,t)/dt = -(1/2)d^2 psi(x,t) /dx^2 + V(x) psi(x,t)

where each of the derivatives is a partial derivative. First, approximate psi(x,t) by only solving on a grid of equally spaced points: x_j = x_0 + j dx. At these points define psi(x_j,t) = psi_j (t) and V(x_j)=V_j. The Schrodinger equation then becomes

i d psi_j(t)/dt=-[psi_{j+1}(t)-2psi_j(t)+psi_{j-1}(t)]/(2 dx^2) + V_j psi_j(t)

Now, the Schrodinger equation is like a first order differential equation in time of a vector with a dimension given by the number of spatial points. Integrals involving psi(x,t) can be approximated using the discretized approximation of psi. For example,

integral |psi(x,t)|^2 dx = dx sum_j |psi_j(t)|^2     and
<x>(t) = integral x |psi(x,t)|^2 dx = dx sum_j x_j |psi_j(t)|^2

There is a difficulty in solving the time dependent Schrodinger equation in that the differential equation is "stiff". Thus if you try the approximation d psi_j(t)/dt = [psi_j(t+dt)-psi_j(t)]/dt with the resulting equation

psi_j(t+dt)=psi_j(t)-i dt {-[psi_{j+1}(t)-2psi_j(t)+psi_{j-1}(t)]/(2 dx^2) + V_j psi_j(t)}

you will find that the norm of the wave function is not constant in time and in fact fairly quickly starts diverging. This is bad. A much better approximation arises when you note that d psi_j(t + dt/2) = [psi_j(t+dt)-psi_j(t)]/dt  is a much better approximation. Also, use psi_j(t + dt/2) = [psi_j(t) + psi_j(t + dt)]/2. These two approximations can be combined to give

psi_j(t+dt) + i dt {-[psi_{j+1}(t+dt)-2psi_j(t+dt)+psi_{j-1}(t+dt)]/(2 dx^2) + V_j psi_j(t+dt)}/2 = psi_j(t)-i dt {-[psi_{j+1}(t)-2psi_j(t)+psi_{j-1}(t)]/(2 dx^2) + V_j psi_j(t)}/2

The right hand side of this equation is something that you can calculate; I call it phi_j

phi_j(t) = psi_j(t)-i dt {-[psi_{j+1}(t)-2psi_j(t)+psi_{j-1}(t)]/(2 dx^2) + V_j psi_j(t)}/2

The wave function at time t+dt can then be obtained from

D_j psi_j(t+dt) + offd psi_{j-1}(t+dt) + offd psi_{j+1}(t+dt) = phi_j(t)

where d_j = 1+(i dt /2)(V_j + dx^-2) and offd = -i dt / (2 dx)^2. This linear equation for psi(t+dt) can be written as

offd psi_2 + D_1 psi_1 = phi_1                           (1)
offd psi_3 + D_2 psi_2 + offd psi_1 = phi_2      (2)
offd psi_4 + D_3 psi_3 + offd psi_2 = phi_3      (3)
etc

These equations can be solved in the following manner. Multiply Eq. (1) by -offd/D_1 and add to Eq. (2) to get

offd psi_2 + D_1 psi_1 = phi_1                           (1)
offd psi_3 + d_2 psi_2  = f_2                              (2')
offd psi_4 + D_3 psi_3 + offd psi_2 = phi_3      (3)
etc

where f_2 = phi_2 - offd phi_1 /D_1 and d_2 = D_2 - offd offd/D_1. Now multiply Eq. (2') by -offd/d_2 and add to Eq. (3) to get

offd psi_2 + D_1 psi_1 = phi_1                           (1)
offd psi_3 + d_2 psi_2  = f_2                              (2')
offd psi_4 + d_3 psi_3  = f_3                              (3')
etc

where f_3 = phi_3 - offd f_2/d_2 and d_3 = D_3 - offd offd/d_2. This can be repeated for all points. The
algorithm can be written as

d_1 = D_1  and f_1 = phi_1
for j greater than or equal to 2 and less then or equal the number of spatial points
f_j = phi_j - offd f_{j-1}/d_{j-1} and d_j = D_j - offd offd/d_{j-1}

This gives the equations

psi_n = f_n/d_n
psi_{n-1} = (f_{n-1} -offd psi_n)/d_{n-1}
psi_{n-2} = (f_{n-2} - offd psi_{n-1})/d_{n-2}
etc

Notice that if you evaluate the equations in the order given then the right hand side in each equation is already calculated by the time it is needed. This sequence of steps gives psi_j(t+dt). An important property of this numerical propagator is that the norm of the wave function does not change with time. An important point to remember with this propagator is that it used psi_0 = 0 and psi_{n+1} = 0. Thus, you are propagating Schrodinger's equation with the potential V(x) with infinite walls at the points x_0 and x_{n+1}=x_0+(n+1)dx.

There are a couple of speed improvements that can be made, but they would make the code a bit more opaque and thus have not been included. The main ones are that the parameters offd/d_j and 1/d_j are recalculated at every time step even though they do not change if the potential is time independent.