The Verlet Method

The Verlet method (also known as leapfrog or Stormer-Verlet) is a second order method that is popular for molecular simulation. It is specialized to problems that can be expressed in the former = v,Mv = F( ), with even dimensional phase space which includes constant energy molecular dynamics. Some generalizations exist for other classes of Hamiltonian systems. 

The Verlet method is a numerical method that respects certain conservation principles associated to the continuous time ordinary differential equations, i.e. it is a geometric integrator. Maintaining these conservation properties is essential in molecular simulation as they play a key role in maintaining the physical environment. As a prelude to a more general discussion of this topic, we demonstrate here that it is possible to derive the Verlet method from the variational principle. This is not the case for every convergent numerical method. The Verlet method is thus a special type of numerical method that provides a discrete model for classical mechanics. 

Mr = F(r) = -V,V(r) where M is a diagonal mass matrix with diagonal 

Remews in Computational Chemistry, Volume 26 edited by Kenny B. Lipkowitz and Thomas R. Cundari Copyright 2009 John Wiley Sons, Inc. 

Computer simulation of the system modeled by Eq. [1] requires some sort of time discretization scheme. The method proposed by Verlet propagates positions by 

Elere the superscripts denote the indices of time steps, each of which is of size h, so 

In contrast to the constant energy regime described above, it is sometimes desirable to perform simulations at a fixed temperature. This can be accomplished by the Langevin dynamics model  

---

In this paper, we discuss semi-implicit/implicit integration methods for highly oscillatory Hamiltonian systems. Such systems arise, for example, in molecular dynamics [1] and in the finite dimensional truncation of Hamiltonian partial differential equations. Classical discretization methods, such as the Verlet method [19], require step-sizes k smaller than the period e of the fast oscillations. Then these methods find pointwise accurate approximate solutions. But the time-step restriction implies an enormous computational burden. Furthermore, in many cases the high-frequency responses are of little or no interest. Consequently, various researchers have considered the use of scini-implicit/implicit methods, e.g. [6, 11, 9, 16, 18, 12, 13, 8, 17, 3]. 

The standard discretization for the equations (9) in molecular dynamics is the (explicit) Verlet method. Stability considerations imply that the Verlet method must be applied with a step-size restriction k < e = j2jK,. Various methods have been suggested to avoid this step-size barrier. The most popular is to replace the stiff spring by a holonomic constraint, as in (4). For our first model problem, this leads to the equations d... 

A timestep of size At with the Verlet method ( velocity Verlet ) takes iqo,Po) to (51,pi) and can be divided into three steps (1) a kick ... 

The standard numerical integrators for the constrained system (12) are the SHAKE scheme [23], which extends the Verlet method (2),... 

Letting m/M 0 in the numerical method, it can be shown that the solution given by (21) tends to a small perturbation of the Verlet method formally applied to that equation ... 

Beeman s algorithm [Beeman 1976] is also related to the Verlet method ... 

In the Verlet method, this equation is written by using central finite differences (see Interpolation and Finite Differences ). Note that the accelerations do not depend upon the velocities. 

The method (2.4) is commonly referred to as StOrmer s rule. It was used by the mathematician Stormer for calculations in the first decade of the 1900s. In molecular dynamics this method is referred to as the Verlet method since it was used by Verlet in his important 1967 paper [387]. 

The most straightforward rewriting of the Verlet method is to put the equations and the discretization in Hamiltonian form, i.e. introducing momenta/> = Mv, and thus pn = Mv , which results in the flow map approximation (taking us from any point in phase space (q,p) to a new point (Q, P) ... 

As another illustration of the performance of the Verlet method, we mention that the dynamical trajectories given in the previous chapter (in particular those given Examples 1.8 and 1.9) were computed using this method (with stepsize h = 0.00001). In the case of the calculation of the exponential separation of... 

All of the methods mentioned so far are said to be explicit discretizations, since producing the next approximation point on a trajectory does not require solving any implicit equations defined in terms of the previous one. The Verlet method is not implicit even though Q appears on the right in the second equation, since it is given in an explicit way in terms of q and/ . Implicitness adds another layer in both analysis and numerical implementation, which, however, in certain applications, is readily justified. 

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

We have mentioned previously that it is possible to reduce the Verlet method to a scheme involving positions only ... 

This method requires that the positions (and forces) be known at two successive points h apart in time in order to initialize the iteration. These might be generated by using the Verlet method or some other self-starting scheme. Beeman s algorithm is explicit since, given q , q - andp , one directly obtains q + and then, q i, and thus p +i, with only one new force evaluation. Because it is a partitioned multistep method, its analysis is more involved than for the one-step methods, and, in particular its qualitative features are difficult to relate to those of the flow map. The order of accuracy of the scheme above can be shown to be three. 

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. 

One might recognize the modification as being proportional to one part of the commutator expansion in the Verlet method, in fact... 

The modified energy for the Verlet method for a single degree of freedom system with energy H = pV2 -h U q) is... 

We now set about comparing various methods. With the Verlet method, it was found that one million timesteps of length about h = 0.01 gave reasonable accuracy for the radial distribution function (see Fig. 3.6). 

Fig. 3.6 The radial distribution computed using the Verlet method is well resolved with a stepsize ofh = 0.01 (even the h = 0.02 is a reasonable approximation)...

Fig. 3.9 The radial distribution computed using the explicit 4th order Runge-Kutta method with stepsizEs of ft = 0.01 and ft = 0.02. Due to energy drift, the ft = 0.02 solution is completely inaccurate, whereas even the ft = 0.01 curve is much worse than that obtained using Verlet with ft = 0.02, despite the fact that the Runge-Kutta method uses four times the number of force evaluations per timestep. These calculations were performed using IM timesteps if a longer simulation were used, the RK4 errors in distribution would be expected to increase (due to the increasing energy drift) whereas the errors reported in Figure 3.6 for the Verlet method would not be expected to depend appreciably on the time interval...

Fig. 3.10 Eneigy as a function of time for the Verlet method with two different stepsizes. The drift is controlled for h = 0.025 but much more dramatic ath = 0.03...

---