Verlet method symplecticness

Fig. 2. The time evolution of the total energy of four water molecules (potential-energy details are given in [48]) as propagated by the symplectic Verlet method (solid) and the nonsymplectic fourth-order Runge-Kutta method (dashed pattern) for Newtonian dynamics at two timestep values.

Example 2.5 (Verlet is Symplectic) For the Symplectic Euler method /, and its adjoint method consider the composition... 

This is an exact transcription of the formulas in Newmark s original paper, only substituting M p for the velocity v wherever it appears. In practice the choice a = 1/2 is used to avoid spurious damping (it can be demonstrated for a simple model problem) this certainly would appear to be desirable in the setting of molecular dynamics. For rj = 0 we then arrive at the Verlet method. For other values of r the scheme is clearly implicit, which likely is the reason it is rarely used in molecular simulation, although it is popular in structural mechanics. The implicit Newmark methods are not symplectic, but a related family of symplectic methods can be constructed by replacing interpolated forces by forces evaluated at interpolated positions [395]. 

More complicated symmetric second-order schemes can be devised by proceeding almost arbitrarily and maintaining a symmetric composition, but it is not found that alternative approaches improve on the two Verlet schemes, at least for molecular applications. The Verlet methods are seen as the gold-standard for molecular dynamics computations both require only one evaluation of VU q) per iteration (where the velocity Verlet scheme can reuse VU Q) for the next iteration), and offer a second-order symplectic evolution. 

Fig. 4.1 Stability regions for Euler s method (l ) and Symplectic Euler/Verlet (right). When a harmonic oscillator is treated using these methods, the origin is unstable for Euler s method, regardless of stepsize—this means that there is no choice of scaling h which will allow us to ensure that 11 + ftA, < 1. On the other hand, the Verlet method has an interval of stability on the imaginary axis, and it is always possible to find a value of h which guarantees that hQ < 2...

The symmetry suggests that this will be second order. In order to understand the symplectic property associated to this method, we need to define the updates for both positions and momenta. A natural choice is to consider the phase space formulation of the Stormer-Verlet method for q = M p,p = F q), then replace F by F —... 

Recall that the Stormer-Verlet method could be constructed by composing steps using Symplectic Euler and its adjoint method. Using more complicated methods it is possible to build higher order schemes. It seems natural that a similar procedure should be possible in the constrained setting. But what, precisely, is the adjoint method in the case of (4.20)-(4.24) ... 

But the methods have not really changed. The Verlet algorithm to solve Newton s equations, introduced by Verlet in 1967 [7], and it s variants are still the most popular algorithms today, possibly because they are time-reversible and symplectic, but surely because they are simple. The force field description was then, and still is, a combination of Lennard-Jones and Coulombic terms, with (mostly) harmonic bonds and periodic dihedrals. Modern extensions have added many more parameters but only modestly more reliability. The now almost universal use of constraints for bonds (and sometimes bond angles) was already introduced in 1977 [8]. That polarisability would be necessary was realized then [9], but it is still not routinely implemented today. Long-range interactions are still troublesome, but the methods that now become popular date back to Ewald in 1921 [10] and Hockney and Eastwood in 1981 [11]. 

This article is organized as follows Sect. 2 explains why it seems important to use symplectic integrators. Sect. 3 describes the Verlet-I/r-RESPA impulse MTS method, Sect. 4 presents the 5 femtosecond time step barrier. Sect. 5 introduce a possible solution termed the mollified impulse method (MOLLY), and Sect. 6 gives the results of preliminary numerical tests with MOLLY. 

Verlet, L. Computer Experiments" on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules. Physical Review 159 (1967) 98-103 Janezic, D., Merzel, F. Split Integration Symplectic Method for Molecular Dynamics Integration. J. Chem. Inf. Comput. Sci. 37 (1997) 1048-1054 McLachlan, R. I. On the Numerical Integration of Ordinary Differential Equations by Symplectic Composition Methods. SIAM J. Sci. Comput. 16 (1995) 151-168... 

An example of a symplectic/time-reversible method is the Verlet (leap-frog) scheme. This method is applicable to separataP Hamiltonian systems of the... 

The symmetry T p) = T[—p) implies that reversing the order of these three steps and changing the sign of r and p results in exactly the same method. In other words, Verlet is time-reversible. (In practice, the equations are usually reduced to equations for the positions at time-steps and the momenta at halfsteps, only, but for consideration of time-reversibility or symplecticness, the method should be formulated as a mapping of phase space.)... 

Fig. 3.2 We compare the symplectic Euler, Verlet and Yoshida schemes in application to a Lennard-Jones oscillator. The plot shows the absolute deviation in the computed Hamiltonian (left) as a function of time, for each scheme at a fixed timestep h = 0.005. Moreover we simulate the system using different stepsizes and compute for each stepsize the mtiximum deviation in the Hamiltonian (right), comparing the results with guide lines associated to various powers of the step size. The scheme used as the base method for the Yoshida composition methods is denoted in the parenthesis (either position or velocity Verlet)...

Among explicit symplectic Partitioned Runge-Kutta methods this is the maximum stability threshold [74]. In a similar way one can analyze the stability of the Verlet and other methods and one thus obtains conditions on the stepsize that must hold for the equilibrium points to be stable in the linearization. Analyzing the stability of both continuous and discrete iteration is much more compUcated for... 

The projected symplectic constrained method (4.20)-(4.24) is only first order accurate. We forego providing a detailed proof of this fact, but note that it could be demonstrated using standard methods [164]. Note that (4.20)-(4.24) reduces to the symplectic Euler method in the absence of constraints, and the projection of the momenta would not alter this fact. There are several constraint-preserving, second-order alternatives which generalize the Stormer-Verlet scheme. One of these is the SHAKE method [322]. The original derivation of the SHAKE method began from the position-only, two-step form of the Stormer rule for q = F(q)... 

As we have seen in the Chap. 3, when a symplectic method is applied to a molecular dynamics problem it induces a perturbed Hamiltonian (energy) function. For the Verlet scheme the modified Hamiltonian is... 

---