introduction to Molecular dynamics, part 2: The MD Algorithm; why ∆t equal to 1 femtosecond?

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.