Integrating the time-dependent Schroedinger equation

Code: Linux version and HP-UX version

Choose an integrator

• The classical Runge-Kutta method (ruku.f);
• The Gear predictor-corrector method (gear.f).

In the code: Change the boolean flags ruku and gear. @See

• dynamics.f
• tsa.f

• main.f

There will always some numerical error no matter how small the steplength was chosen to be or how (relatively) accurate the algorithm was devised to be. (We are not hearing a lot of complaints from the molecular simulations community about this problem. Because statistic mechanics is fault tolerant, this problem is not very pressing for molecular dynamics simulations with interests in statistical results instead of individual dynamics. CHARMm's temperature control by use of rescaling or reassigning of velocities is a rough approach, so nothing should be conserved or maintained. But if one uses constraint methods or the Nose-Hoover thermostat, the integration accuracy should be checked.) In electron transfer dynamics, it is important not to violate the conservation law too much, since we are simulating individual processes and the accuracy of the conclusions will be drawn on the basis of that of these individual processes. Because of this consideration, we impose a rescaling scheme that rectifies the population distribution once the violation of conservation law reaches some sort of tolerance. This rescaling procedure will renormalize the orbital populations from time to time, hence (hopefully) diligently prevent the numerical error from propogating.

Quantum mechanics requires that the integration algorithm should also be time reversible. But this is not strictly needed in our simulations.

Saving a restart file

• Gear:save(n,time,dotsa)
• Gear:readrst(n,time0,psi0,dotsa)
• RK4:saverkns(n,time,psi,dotsa)
• RK4:readrkns(n,time0,psi0,dotsa)

In our program, the intermediate Hamiltonian matrix elments are obtained by using the cubic spline interpolation method (see subroutine interpol(nstp,ndim,tmd,h,cof) in util.f, where nstp is the number of discret input points, ndim is the dimension of the matrix to be interpolated h, cof is the array that stores the cubic spline coefficients). An interpolation coefficients table (the array cof) is prepared before integrating the equation of motion, and used in the integration procedure.

Because we have used SHAKE to constrain the bondlengths of the covalent bonds that involve hydrogen atoms in the MD simulations for azurin (therefore actually remove motion of the highest frequency), the choice of the interpolation interval, i.e. the MD time steplength, makes missing the fastest MD motion of the protein less likely.

The spline interpolation method has essentially a boundary problem with respect to where the starting point is to determine the intermediate Hamiltonian matrix elements inside a given interval. We know that the cubic spline interpolation coefficients are decided based on the three nearest neighboring discrete points. Let's say we have four discrete points: 1, 2, 3, 4. Let's first set the starting point at point 1, we have a set of spline coefficents for the intervals [1,2), [2,3) and [3,4). Now, let's select point 2 as the starting point, we will have another set of spline coefficents for the intervals [2,3)' and [3,4)' (prime for distinguishing the two resulting interpolated Hamiltonian time series). We know [3,4)=[3,4)', but [2,3)'!=[2,3), because the basis for determining the spline coefficents has been changed. Here is a postscript picture about how much the deviation is. In order to avoid this problem, the program always reads the Hamiltonian at a step prior to the actual ED starting time just for guaranteeing the agreement of interpolation in the joint area between the present segment and the previous one. (see nstp0 in hamiltonian.f.)

Propogation on the nonorthogonal basis