Symplectic

The velocity Verlet algorithm may be derived by considering a standard approximate decomposition of the Liouville operator which preserves reversibility and is symplectic (which implies that volume in phase space is conserved). This approach [47] has had several beneficial consequences. 

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]. 

In practice modifications are made to incorporate thermostats or barostats that may destroy the time-reversible and symplectic properties. While extended-system algorithms such as Nose dynamics [41] can be designed on the principles of the reversible operators, methods that use proportional velocity or coordinate scaling [42] cannot. Such methods arc very... 

R. D. Skeel, G. H. Zhang, and T. Schlick. A family of symplectic integrators Stability, accuracy, and molecular dynamics applications. SIAM J. Scient. COMP., 18 203-222, 1997. 

Robert D. Skeel, Jeffrey J. Biesiadecki, and Daniel Okunbor. Symplectic integration for macromolecular dynamics. In Proceedings of the International Conference Computation of Differential Equations and Dynamical Systems. World Scientific Publishing Co., 1992. in press. 

Robert D. Skeel and Jeffrey J. Biesiadecki. Symplectic integration with variable stepsize. Ann. Num. Math., 1 191-198, 1994. 

Backward Analysis In this type of analysis, the discrete solution is regarded as an exact solution of a perturbed problem. In particular, backward analysis of symplectic discretizations of Hamiltonian systems (such as the popular Verlet scheme) has recently achieved a considerable amount of attention (see [17, 8, 3]). Such discretizations give rise to the following feature the discrete solution of a Hamiltonian system is exponentially close to the exact solution of a perturbed Hamiltonian system, in which, for consistency order p and stepsize r, the perturbed Hamiltonian has the form [11, 3]... 

Summarizing, from a mathematical point of view, both forward and backward analysis lead to the insight that long term trajectory simulation should be avoided even with symplectic discretizations. Rather, in the spirit of multiple as opposed to single shooting (cf. Bulirsch [4, 18]), only short term trajectories should be used to obtain reliable information. 

We restrict our attention to symplectic one-step discretizations of (1), which leads to discrete dynamical systems of the form... 

This section deals with the question of how to approximate the essential features of the flow for given energy E. Recall that the flow conserves energy, i.e., it maps the energy surface Pq E) = x e P H x) = E onto itself. In the language of statistical physics, we want to approximate the microcanonical ensemble. However, even for a symplectic discretization, the discrete flow / = (i/i ) does not conserve energy exactly, but only on... 

Inefficiency of Direct Simulation Suppose we want to compute the corresponding invariant measure p by direct simulation. Direct long term simulation by symplectic discretization of (1) yields the discrete solution ( ) a )i, ... 

G. Benettin and A. Giorgilli. On the Hamiltonian interpolation of near to the identity symplectic mappings with applications to symplectic integration algorithms. J. Stat. Phys. 74 (1994)... 

E. Hairer. Backward analysis of numerical integrators and symplectic methods. Annals of Numerical Mathematics 1 (1994)... 

When the integrator used is reversible and symplectic (preserves the phase space volume) the acceptance probability will exactly satisfy detailed balance and the walk will sample the equilibrium distribution... 

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.

This discussion suggests that even the reference trajectories used by symplectic integrators such as Verlet may not be sufficiently accurate in this more rigorous sense. They are quite reasonable, however, if one requires, for example, that trajectories capture the spectral densities associated with the fastest motions in accord to the governing model [13, 15]. Furthermore, other approaches, including nonsymplectic integrators and trajectories based on stochastic differential equations, can also be suitable in this case when carefully formulated. 

The implicit-midpoint (IM) scheme differs from IE above in that it is symmetric and symplectic. It is also special in the sense that the transformation matrix for the model linear problem is unitary, partitioning kinetic and potential-energy components identically. Like IE, IM is also A-stable. IM is (herefore a more reasonable candidate for integration of conservative systems, and several researchers have explored such applications [58, 59, 60, 61]. 

Note that the method is implicit unless q = 0 symplecticness is proven in [66, 67]. For these schemes, we obtain... 

The IE and IM methods described above turn out to be quite special in that IE s damping is extreme and IM s resonance patterns are quite severe relative to related symplectic methods. However, success was not much greater with a symplectic implicit Runge-Kutta integrator examined by Janezic and coworkers [40],... 

In LN, the bonded interactions are treated by the approximate linearization, and the local nonbonded interactions, as well as the nonlocal interactions, are treated by constant extrapolation over longer intervals Atm and At, respectively). We define the integers fci,fc2 > 1 by their relation to the different timesteps as Atm — At and At = 2 Atm- This extrapolation as used in LN contrasts the modern impulse MTS methods which only add the contribution of the slow forces at the time of their evaluation. The impulse treatment makes the methods symplectic, but limits the outermost timestep due to resonance (see figures comparing LN to impulse-MTS behavior as the outer timestep is increased in [88]). In fact, the early versions of MTS methods for MD relied on extrapolation and were abandoned because of a notable energy drift. This drift is avoided by the phenomenological, stochastic terms in LN. 

The heightened appreciation of resonance problems, in particular, has been quite recent [63, 62], and contrasts the more systematic error associated with numerical stability that grows systematically with the discretization size. Ironically, resonance artifacts are worse in the modern impulse multiple-timestep methods, formulated to be symplectic and reversible the earlier extrapolative variants were abandoned due to energy drifts. 

B. Leimkuhler and R. D. Skeel. Symplectic numerical integrators in constrained Hamiltonian systems. J. Comp. Phys., 112 117-125, 1994. 

J. C. Simo, N. Tarnow, and K. K. Wang. Exact energy-momentum conserving algorithms and symplectic schemes for nonlinear dynamics. Computer Methods in Applied Mechanics and Engineering, 100 63-116, 1994. 

O. Gonzales and J. C. Simo. On the stability of symplectic and energy-momentum conserving algorithms for nonlinear Hamiltonian systems with symmetry. Comp. Meth. App. Mech. Engin., 134 197, 1994. 

M. Zhang and R. D. Skeel. Cheap implicit symplectic integrators. Applied Numerical Mathematics, 25 297-302, 1997. 

D. Janezid and F. Merzel. An efficient symplectic integration algorithm for molecular dynamics simulations. J. Chem. Info. Comp. Set, 35 321-326, 1995. 

S. Reich. Preservation of adiabatic invariants under symplectic discretization. i4pp. Numer. Math., to appear, 1998. 

J.C. Simo and O. Gonzales. Assessment of energy-monentum and symplectic schemes for stiff dynamical systems. The American Society of Mechanical Engineering, 93-WA/PVP-4, 1993. 

One of the advantages of the Verlet integrator is that it is time reversible and symplectic[30, 31, 32]. Reversibility means that in the absence of numerical round off error, if the trajectory is run for many time steps, say nAt, and the velocities are then reversed, the trajectory will retrace its path and after nAt more time steps it will land back where it started. An integrator can be viewed as a mapping from one point in phase apace to another. If this mapping is applied to a measurable point set of states at on(> time, it will... 

Thus these integrators are measure preserving and give trajectories that satisfy the Liouville theorem. [12] This is an important property of symplectic integrators, and, as mentioned before, it is this property that makes these integrators more stable than non-symplectic integrators. [30, 33]... 

Since many systems of interest in chemistry have intrinsic multiple time scales it is important to use integrators that deal efficiently with the multiple time scale problem. Since our multiple time step algorithm, the so-called reversible Reference System Propagator Algorithm (r-RESPA) [17, 24, 18, 26] is time reversible and symplectic, they are very useful in combination with HMC for constant temperature simulations of large protein systems. 

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. 

---