Verlet equation

Algorithms in this family are simple, accurate, and, as we will see below, time reversible. They are the most widely used methods for integrating the classical equations of motion. The initial form of the Verlet equations [3] is obtained by utilizing a Taylor expansion at times t — dt and t + dt... 

This is the Verlet equation for integrating the coordinates and, by renaming h as f and reassigning t At t, one obtains the usual expression ... 

A straightforward derivation (not reproduced here) shows that the effect of the diree successive steps embodied in equation (b3.3.7), with the above choice of operators, is precisely the velocity Verlet algorithm. This approach is particularly usefiil for generating multiple time-step methods. 

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

Verlet, L. Computer experiments on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules. Phys. Rev. 165 (1967) 98-103. Ryckaert, J.-P., Ciccotti,G., Berendsen, H.J.C. Numerical integration of the cartesian equations of motion of a system with constraints Molecular dynamics of n-alkanes. Comput. Phys. 23 (1977) 327-341. 

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. 

There is a number of algorithms to solve equations (1) and (2) that differ appreciably in their properties which are beyond the scope of the present article. In the discussion below we use the velocity Verlet algorithm. However, better approaches can be employed [2-5]. We define a rule - F X t), At) that modifies X t) to X t + At) and repeat the application of this rule as desired. For example the velocity Verlet algorithm ( rule ) is ... 

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

The constrained equations of motion in cartesian eoordinates can be solved by the SHAKE or (the essentially equivalent) RATTLE method (see [8]) which requires the solution of a non-linear system of equations in the Lagrange multiplier funetion A. The equivalent formulation in local coordinates ean still be integrated by using the explicit Verlet method. 

The idea is illustrated by Fig. 1. These equations constitute a readily understandable and concise representation of the widely used Verlet-I/r-RESPA impulse MTS method. The method was described first in [8, 9] but tested... 

The simplest of the numerical techniques for the integration of equations of motion is leapfrog-Verlet algorithm (LFV), which is known to be symplectic and of second order. The name leapfrog steams from the fact that coordinates and velocities are calculated at different times. 

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

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

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

The velocity Verlet method is actually implemented as a three-stage procedure because, as can be seen from Equation (7.15), to calculate the new velocities requires the accelerations at both t and t + 8t. Thus in the first step the positions at f I- are calculated according to Equation (7.14) using the velocities and the accelerations at time t. The velocities at time t + 6t are then determined using ... 

The equations of motion are integrated using a modified velocity Verlet algorithm. The modification is required because the force depends upon the velocity the extra step involves... 

Dissipative particle dynamics (DPD) is a technique for simulating the motion of mesoscale beads. The technique is superficially similar to a Brownian dynamics simulation in that it incorporates equations of motion, a dissipative (random) force, and a viscous drag between moving beads. However, the simulation uses a modified velocity Verlet algorithm to ensure that total momentum and force symmetries are conserved. This results in a simulation that obeys the Navier-Stokes equations and can thus predict flow. In order to set up these equations, there must be parameters to describe the interaction between beads, dissipative force, and drag. 

If the magnitudes of the dissipative force, random noise, or the time step are too large, the modified velocity Verlet algorithm will not correctly integrate the equations of motion and thus give incorrect results. The values that are valid depend on the particle sizes being used. A system of reduced units can be defined in which these limits remain constant. 

Again, elimination of the velocities from these equations recovers the Verlet algorithm. In practice, the velocity Verlet algorithm consists of the following steps ... 

The g r) that results from the modified Verlet (MV) closure is very close to the simulation results in Figs. 2 and 3. The MV results for g d), or equivalently, y d), are plotted in Fig. 4(a). The resulting equation of state is similar to the CS expression. An even more demanding test is an examination of the MV results for y r) for r < d. As is seen in Fig. 4(b), the MV results for y(0) are quite good [25], and are better than the PY and HNC results. Some results have also indicated that the MV closure gives quite accurate results for a mixture of hard spheres [26]. 

MD runs for polymers typically exceed the stability Umits of a micro-canonical simulation, so using the fluctuation-dissipation theorem one can define a canonical ensemble and stabilize the runs. For the noise term one can use equally distributed random numbers which have the mean value and the second moment required by Eq. (13). In most cases the equations of motion are then solved using a third- or fifth-order predictor-corrector or Verlet s algorithms. 

Adding these two equations gives the Verlet algorithm, which is used to advance the position vector r from its value at time t to time f + Sf... 

This is the Verlet algorithm for solving Newton s equation numerically. Notice that the term involving the change in acceleration (b) disappears, i.e. the equation is correct to third order in At. At the initial point the previous positions are not available, but may be estimated from a first-order approximation of eq. (16.29). 

In our calculations we make use of several standard techniques of molecular dynamics simulations. The integration of the equations of motions is done by the velocity form of the Verlet-algorithm with a time step of 1.5 The temperature is controlled... 

Thus, positions of the particles at the time t + At can be computed from the positions at times t and t - At, and the second derivative, (frildfi) that corresponds to the acceleration. The latter can be computed via Eqs. (35.1) and (35.2). Equation (35.3) is known as the Verlet algorithm. A number of methods are discussed in specihc textbooks (e.g., Allen and Tildesley, 1992), and the choice must be taken by making a balance between accuracy and cost in execution speed. 

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. 

---