A far more efficient algorithm (by a factor of four) than the Runge-Kutta
algorithm is due to Gear [1970]. This is a multistep method, requiring
higher order derivatives of 
to be known. As a consequence, a second
order Runge-Kutta algorithm must be used to evaluate the higher order
derivatives at time t=0. This was tried for Gaussian systems, and was found
to give the correct trajectory 99.5% of the time. The cause of the erroneous
trajectories was never resolved. It was thought that the error involved in
neglecting was statistically insignificant, however, the value calculated for
colour conductivity disagreed with the accepted value by about 10%, or 2
standard errors. When the theory based on a single equilibrium trajectory was
developed, it no longer mattered what the initial condition of the system are,
so a self-starting method was no longer required.
A listing of the Gear integrator which is 4
order in h follows. The
SPM_quot1" in the name refers to the fact that a first order differential
equation is being solved. Algorithms are listed in Gear's book for second and
higher order equation. The SPM_quotZ" refers to the equations of motion being
integrated as the current and thermo-statted system listed as SPM_quotNORTON".
It was found that in order to efficiently vectorize the code on the Fujitsu
VP100, each subroutine had to be expanded inline at the compilation stage
in order to avoid any recursive data dependencies (procedure integration).
Unfortunately, procedure integration could not be performed on functional
parameters, so the name of the equations of motion function had to be hard
coded into the routine. For similar reasons, the phase space vector
SPM_quotX" is now passed as a common block, so the only parameters left are
SPM_quotH" and SPM_quotNTSTEP". To alter this routine to do anything else only
requires the strings SPM_quotNORTON", SPM_quotGAMMAG" and SPM_quotDIMPS+2" to be
globally changed with a text editor.
C Solve the equations of motion by the 4th order Gear predictor-
C corrector method. Gear, C.W. (1971) Numerical Initial Value Problems
C in Differential Equations. (Prentice Hall: Englewood Cliffs, N.J.)
C This subroutine treats the equations of motion as a 1st order d.e.
SUBROUTINE GEAR1Z( H, NTSTEP )
INCLUDE (PARAM)
C XT = x(t) = final state after time t = H*NTSTEP.
C H = step size
C NTSTEP = Number of time steps
C n
C H n n
C Xn = -- (d /dt ) XT
C n!
C
C XnPRED = the predicted values of the above using a truncated Taylor
C series
INTEGER NTSTEP, I, T
DOUBLE PRECISION H, CORR, NORTON
DOUBLE PRECISION XT(DIMPS+2), XTPRED(DIMPS+2), X3PRED(DIMPS+2)
DOUBLE PRECISION X2(DIMPS+2), X2PRED(DIMPS+2), X3(DIMPS+2)
DOUBLE PRECISION X4(DIMPS+2)
DOUBLE PRECISION X1(DIMPS+2), X1PRED(DIMPS+2)
C Local storage is shared with SETGEG, GEARN, SETGEN
COMMON /GSTORE/ XTPRED,X1PRED,X2PRED,X3PRED
C Derivatives filled of SETGEG
COMMON /GSTORZ/ X1,X2,X3,X4
COMMON /GAMMAG/ XT
C Gear coefficients
DOUBLE PRECISION F0,F1,F2,F3,F4
PARAMETER( F0=251./720, F1=1, F2=11./12, F3=1./3, F4=1./24)
C These statements are required for Procedure Integration on the FACOM
DOUBLE PRECISION F(D*NP),V
COMMON /FORCE/F,V
C Calculation loop
DO 2 T=1,NTSTEP
C Predictor loop
*VOCL LOOP,NOVREC(XTPRED)
DO 3 I=1,DIMPS+2
XTPRED(I) = X1(I) + XT(I)+X2(I)+X3(I)+X4(I)
X1PRED(I) = X1(I) +
& 2 * X2(I) + 3 * X3(I) + 4 * X4(I)
X2PRED(I) = X2(I) + 3*X3(I) + 6*X4(I)
X3PRED(I) = X3(I) + 4*X4(I)
3 CONTINUE
C Evaluation step
CALL FORCES(XTPRED)
C Corrector loop
*VOCL LOOP,NOVREC(XT)
DO 4 I=1,DIMPS+2
CORR = NORTON(I,XTPRED) * H - X1PRED(I)
XT(I) = XTPRED(I) + F0 * CORR
X1(I) = X1PRED(I) + F1 * CORR
X2(I) = X2PRED(I) + F2 * CORR
X3(I) = X3PRED(I) + F3 * CORR
X4(I) = X4(I) + F4 * CORR
4 CONTINUE
2 CONTINUE
CALL CUBE( XT)
RETURN
END