###
2.2 The MD Algorithm

####
2.2.1 Verlet’s Algorithm and Velocity Verlet’s Algorithm

The governing equations of motions are all ordinary differential equations and a

standard method for the solutions is the finite difference method. Of a large number of

finite-difference approaches devised, the most commonly used one is the Runge-Kutta

(RK) method. [33,34] The beauty of the RK method is that the order of a finite

difference method can often be increased by using positions and velocities from

several points in time instead of only the current time. However, if a fourth-order RK

method is applied to MD simulation, four interatomic force evaluations per atom

would be required per step, which will make the simulation too slow. Thus, reduction

of order is desired; especially, the time expense of calculating intermolecular force

more than once per atom per step should be also avoided. However, the algorithm

must be stable.

The most widely used method of integrating the equations of motions in molecular

dynamics is Verlet‘s method or algorithm. It uses positions and velocities from

previously evaluated steps to provide a direct solution of equation (2.4). In this study,

the simulation code LAMMPS is used and its default algorithm is the velocity Verlet.

[35] Hence, the basic idea of Verlet‘s algorithm, and specifically the velocity Verlet,

will be presented here. Information on other integration schemes can be found in texts

on atomic scale simulation. [32]

This standard, or basic, Verlet algorithm is a combination of two Taylor‘s

expansions as following. The Taylor series for position from time t forward to 𝑡 + ∆𝑡

And the series form t backward to 𝑡 − ∆𝑡 is

Adding equations (2.6) and (2.7) with some simplifying produces

(2.8)

This provides Verlet‘s algorithm for positions. The truncation error of the algorithm

when evolving the system by ∆t is of the order of ∆t^4. The acceleration needed in
equation (2.8) is obtained from the intermolecular forces and Newton‘s second law.
Velocities are not explicitly solved in equation (2.8), and they are calculated from
first-order central difference.

A more commonly used, related algorithm is the velocity Verlet algorithm. The

approach is similar with the standard one but starts with position and velocity

expansions explicitly

The implement scheme of this algorithm is following:

1. Calculate velocities at mid-step,

Calculate the positions at next step 𝑡 + ∆t

2. Calculate the accelerations at 𝑡 + ∆𝑡, 𝑎(𝑡 + ∆𝑡) from the interaction
potential using 𝑥(𝑡 + ∆𝑡).

3. Calculate the velocities at 𝑡 + ∆t

The advantage of this algorithm compared to the standard one is that we never need to simultaneously store the values at two different times for any of these

quantities. The velocity Verlet algorithm gives good feature of time reversibility,

simplicity and stability for even larger time steps compared to standard algorithm.

Ref:

http://preserve.lehigh.edu/cgi/viewcontent.cgi?article=2257&context=etd
*Molecular Dynamics Simulation of Gold Nanoparticles and Surface Stress Effect*
*Siming Zhang*
*Lehigh University*
*2011*

*---------------------------------------------------*
Verlet is an efficient algorithm, and its benefits are stability, time-reversibility and

energy conservation.

**The time step, ∆t, should be chosen in such manner that the simulation system**
**is stable and energy is conserved, that is to say that the provided ∆t should be**
**short enough. Usually, in classical MD simulations time steps of order 1–2 fs **are

being used (especially when using constraints), but in some cases time steps can

be even 4–5 fs (when using virtual site algorithms). The length of the time

step is dominated by the fastest degrees of motion in the system, and usually theses

are bond vibrations. Most often constraints are used in simulations to reduce the

degrees of freedom, removing high-frequency vibrations by restricting the movement

of covalent bonds. This is implemented by fixing the covalent bonds to constant

angles and lengths. In addition to constraints, also restraints can be used to fix

atomic positions, bond lengths, bond angles and dihedrals.