Bit-Reversible Version of Milne's Fourth-Order Time-Reversible Integrator for Molecular Dynamics

TL;DR

This paper presents a bit-reversible implementation of Milne's fourth-order integrator, enhancing accuracy and stability for molecular dynamics simulations, especially useful for analyzing Lyapunov instability in dynamical systems.

Contribution

It introduces a bit-reversible version of Milne's integrator, improving accuracy and simplifying velocity and energy definitions for molecular dynamics.

Findings

Enhanced accuracy over Verlet's algorithm

Simplified velocity and energy calculations

Effective for Lyapunov instability analysis

Abstract

We point out that two of Milne's fourth-order integrators are well-suited to bit-reversible simulations. The fourth-order method improves on the accuracy of Levesque and Verlet's algorithm and simplifies the definition of the velocity $v$ and energy $e = (q^2 + v^2)/2$ . ( We use this one-dimensional oscillator problem as an illustration throughout this paper ). Milne's integrator is particularly useful for the analysis of Lyapunov ( exponential ) instability in dynamical systems, including manybody molecular dynamics. We include the details necessary to the implementation of Milne's Algorithms.

Equations

$$q_{t+dt}-2q_t+q_{t-dt}=a_t(dt{)}^2.$$

$$q_{t+(dt/2)}=q_t+v_t(dt/2);v_{t+dt}=v_t+a_{t+(dt/2)}dt;q_{t+dt}=q_{t+(dt/2)}+v_{t+dt}(dt/2),$$

$$q_{t+2dt}-q_{t+dt}-q_{t-dt}+q_{t-2dt}=[5a_{t+dt}+2a_t+5a_{t-dt}](d{t}^2/4).$$

$$q_{t+2dt}-2q_{t+dt}+q_t=[a_{t+2dt}+10a_{t+dt}+a_t](d{t}^2/12).$$

$$\cos (2\omega dt)-\cos (\omega dt)+(d{t}^2/4)[5\cos (\omega dt)+1]=0.$$

$$2C^2+[(5d{t}^2/4)-1]C+(d{t}^2/4)-1=0\mathrm{where}C\equiv \cos (\omega dt).$$

Full text

Bit-Reversible Molecular Dynamics, Lyapunov Instability, Chaotic Dynamical Systems

I Introduction

William Milne’s 1949 work Numerical Calculusb1 was republished by the Princeton University Press in 2015. The book is a particularly valuable source of clear and direct numerical methods. Research workers in statistical mechanics, molecular dynamics, and dynamical systems will find his approach to what is our own research interest, solving and analyzing differential equations for chaotic systems small and large, reliable and useful. Writing about a decade prior to the computer revolution Milne had no particular interest in “reversible computing” and the “bit-reversible” algorithms which make it possible to extend sequences of coordinates forward and backward in time stably and reversibly ad infinitum. Nevertheless his work is directly applicable to such finite-difference applications.

In 1993 Dominique Levesque and Loup Verlet used an integer algorithm to solve problems in Newtonian mechanics with perfect time reversibilityb2 . Loup had popularized Størmer and Newton’s Leapfrog Algorithm a quarter century earlier, in the early days of molecular dynamicsb3 ,

[TABLE]

“Verlet’s algorithm” appears on page 140 of Reference 1. If the righthand side of this finite-difference algorithm is truncated to an integer the resulting acceleration is precisely the same ( to the very last computational “bit” ) in either direction of time. Because this algorithm conserves phase volume when written in a “symplectic” centered-difference form :

[TABLE]

there is no tendency for energy drift. The errors in the velocity and energy in the leapfrog algorithm are unnecessarily large, so that two of Milne’s algorithms ( both of them also on page 140 of Reference 1 ) can provide better accuracy for longer runs :

[TABLE]

The error, $\simeq(17/240)dt^{6}$, is quite tolerable relative to the Størmer error, $\simeq(1/12)dt^{4}$. Milne also gives an even better corrector formula with an error $\simeq(-1/240)dt^{6}$ .

[TABLE]

II Applications

For several years nowb4 ; b5 we have been exploring the differences in Lyapunov spectra forward and backward in time in order to get insight into the Second Law of Thermodynamics. The fractal structures which arise in nonequilibrium deterministic and time-reversible steady-state problems provide explanations to both Loschmidt’s Reversibility paradox and Zermélo’s Recurrence paradoxb6 . Levesque and Verlet’s integer algorithm has proved to be a useful tool in these studies despite its relatively coarse description of particle trajectories.

Integer algorithms are also useful in studies of the effects of finite precision ( single, double, quadruple, … ) on phase-space distributions generated by flows and mapsb7 ; b8 . Mauricio Romero-Bastidab9 suggested the use of the integer-based leapfrog algorithm for generating a reversible reference trajectory of arbitrary length in his studies of the “covariant” Lyapunov exponents. In 2013 we were able to see a qualitative difference between the “important particles” ( those making above-average contributions to the Lyapunov instability ) forward and backward in time in the example inelastic-collision problem of Figure 1b10 . Continuing progress in low-cost computation caused us to revisit these problems in connection with a lecture course delivered at Kharagpur’s Indian Institute of Technology in December 2016b11 . We were very pleased to find that Milne’s work offers an improvement in the precision and accuracy of these Lyapunov studies and believe that others will find this approach useful to their own work. Although these improvements are not at all “new” we do expect that this work will accelerate progress in understanding the time-reversible simulation of irreversible processes.

Figure 1 illustrates the important particles ( those making above-average contributions to the largest Lyapunov exponent ) forward and backward in time for the collision of two 400-particle balls. The 162 important particles forward in time are those blacked in along the interface between the balls while the 120 important particles backward in time are those blacked in in the necking regions where the plastic strain is greatest as the balls are separating. This simulation employed the Levesque-Verlet algorithm for the reference trajectory and a Runge-Kutta fourth-order algorithm for the two satellite trajectories ( one forward and another backward ) as is described in Reference 10.

III Numerical Implementation of Milne’s Algorithms

To illustrate the application of Milne’s Algorithms we consider an integer version of the simpler of his two fourth-order algorithms. We describe a harmonic oscillator with $\ddot{q}=-q$ . The preliminaries, which we give below, provide integer forms for five previous coordinates and the corresponding contributions to the acceleration, all of them multiplied by $10^{15}$. We select an example timestep of $(\pi/50)$ in the Fortran instructions so that an oscillator period corresponds to one hundred timesteps.

We carried out two kinds of tests for the Milne integrator, reversibility, confirming that reversing the four prereversal coordinates exactly reverses the sequence of integers back to the initial value of $10^{15}$. It is easy to show that the algorithm is exactly reversible in this way. Stability can be confirmed by solving for the dependence of the oscillation frequency on the timestep. Numerical work consistent with the linear analysis for the oscillator ( given in more detail in a Postscript ) shows that the dependence of the phase shift is quartic in the timestep for the range $0<dt<1$. A direct simulation of the integer version of the algorithm, using two billion timesteps with $dt=0.2$, showed no tendency toward damping or instability. Similar results can be obtained by solving the floating-point version of the problem, where precise reversibility has to be abandoned ( because roundoff error will spoil it ). These results establish that the Milne algorithm is both reversible and stable for the oscillator. We recommend it to our colleagues for their use.

The implementation of the algorithm is to some extent hardware dependent. On our various Mac computers using the free gnu compiler we had no trouble using 16-byte integers, giving roughly 15 digits for arithmetical operations. The following extract from the setup of the computation generates the inital data ( in this case four points from a cosine curve ) as well as the three integers, proportional to $dt^{2}\times 10^{15}$, and needed for the accelerations. On the following page we summarize the time-stepping loop where the three accelerations are expressed as integers $\{{\tt IAP,IA0,IAM}\}$. We include at the end an indication of the coordinate reversal procedure needed to integrate backward.

IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER16 IQMM,IQM,IQ0,IQP,IQPP ! Contiguous integer coordinates INTEGER16 IAP,IA0,IAM ! Ingredients of the accelerations INTEGER16 I1,I2,I3,I4 ! Storage for coordinate reversal ITMAX = 100 TWOPI = 2.D03.141592653589793D0 DT = TWOPI/ITMAX QMM = DCOS(-2.D0DT) QM = DCOS(-1.D0DT) Q0 = DCOS( 0.D0DT) QP = DCOS(+1.D0DT) CONVERT COORDINATES TO 15-DIGIT INTEGERS IQMM = QMM*(10.D015) IQM = QM *(10.D015) IQ0 = Q0 *(10.D015) IQP = QP *(10.D015) ! Finish of the preliminaries

Here follows a bare-bones evolution loop for the integer coordinates. After ITMAX iterations the coordinate reversal steps make it possible to return precisely to, and beyond, the beginning. In the event that the velocities are to be calculated from Milne’s fifth-order interpolation ( which is one order of overkill ) it is necessary to compute the integer coordinates to include IQPPP and IQMMM, getting IQPPP from the “step” IQPPP = IQPP + IQ0 - IQM - (IAPP + IAP + IA0).

  DO IT = 1,ITMAX
  TIME = IT*DT

COMPUTE INGREDIENTS OF THREE ACCELERATIONS IAP = 0.25D00DTDT*(5.D0IQP) IA0 = 0.25D00DTDT(2.D0IQ0) IAM = 0.25D00DTDT(5.D0*IQM) COORDINATE UPDATES FOR FIVE SUCCESSIVE TIMES IQPP = IQP + IQM - IQMM - (IAP + IA0 + IAM) ! This is the "step". IQMM = IQM IQM = IQ0 IQ0 = IQP IQP = IQPP END DO COORDINATE REVERSAL I1 = IQMM I2 = IQM I3 = IQ0 I4 = IQP IQMM = I4 IQM = I3 IQ0 = I2 IQP = I1 DO IT = 1,IMAX CONTINUE REVERSAL WITH SAME STATEMENTS AS IN THE FORWARD LOOP END DO

There is no difficulty in computing an accurate velocity with Milne’s page 99 formula using six centered coordinates. This fifth-order interpolation gives not only good velocities, but also an accurate energy. Accurate values of these phase variables are a real advantage of the Milne algorithm over that of Levesque and Verlet.

IV Postscript on the Stability of Milne’s Algorithm

The stability analysis for Milne’s algorithm is straightforward. If we substitute the trial solution $q\propto e^{i\omega t}$ into the wholly linear algorithm the result is :

[TABLE]

This simplifies to a quadratic equation in $\cos(\omega dt)$ :

[TABLE]

Figure 2 shows that the dependence of the frequency error on the timestep is quartic, $(1-\omega)\propto dt^{4}$, confirming the stability of the algorithm.

V Acknowledgement

The interest and support of Harald Posch ( Universität Wien ), Clint Sprott ( University of Wisconsin-Madison ), and Karl Travis ( University of Sheffield ) are gratefully acknowledged.