Introduction¶
In molecular dynamics (MD) simulations, the program needs to compute the evolution of
the atoms' coordinates, velocities, and accelerations. The atoms' coordinates contain the information
about the structure of the molecular or material systems. Trajectory files containing a time series of atomic position, position, and force data are typically outputted from an MD simulation. Researchers can use visualization programs like
Visual Moleuclar Dyanmics to view the evolution
of the system and make scientific conclusions based on the results. Hence, it is critical
to calculate the atoms' position, velocity, and acceleration accurately.
In classic molecular dynamics simulations, the molecular motions are computed by
solving the Newton equation. A numerical approximation is employed, wherein one treats the
velocities of the atoms as constant during a tiny period $\Delta t$. Here, we used a two-body
system (sun and earth) to demonstrate this numerical method.
Verlet algorithm
Most molecular dynamics programs employ the algorithm adopted
by Loup Verlet
[Phys. Rev. 159, 98, (1967)]. There are many variants of the Verlet
algorithm. One of the most revealing methods is the so-called
velocity Verlet algorithm, which was proposed by William C. Swope
and Hans C. Andersen
[J. Chem. Phys. 76, 637 (1982)]. There are three equations to demonstrate
the velocity Verlet algorithm.
$$\normalsize V(t+\frac{1}{2} \Delta t)=v(t)+\frac{1}{2}\Delta a(t) \quad (1)$$
$$\normalsize r(t+\Delta t)=r(t)+\Delta tv(t+\frac{1}{2}\Delta t) \quad (2)$$
$$\normalsize v(t+\Delta t)=v(t+\frac{1}{2}\Delta t)+\frac{1}{2}\Delta t a(t+\Delta t) \quad (3)$$
Eq. (1) shows that the velocities of the atoms advance half of the time
step $\Delta t$. The new positions of the atoms are updated from previous
positions with constant velocities $v(t+\frac{1}{2}\Delta t)$ as shown
in Eq. (2). The velocities are updating from the middle of the time step
$v(t+\frac{1}{2}\Delta t)$.
In the original implementation of the Verlet algorithm, the positions of
the atoms are expressed via the Taylor expansion:
$$\normalsize r(t+\Delta t)=r(t) + \Delta tv(t) + \frac{1}{2}\Delta t^2a(t) + \dots \quad (4)$$
$$\normalsize r(t-\Delta t)=r(t) - \Delta tv(t) + \frac{1}{2}\Delta t^2a(t) - \dots \quad (5)$$
Adding Eq. (4) to Eq. (5):
$$\normalsize r(t+\Delta t)=2r(t)-r(t-\Delta t) + \Delta t^2a(t) + O(\Delta t^4)$$
The velocities are not needed to compute the trajectories of the atoms.
However, the velocities are useful to estimate the kinetic energy. They
can be approximated as:
$$v(t) = \frac{r(t+\Delta t)-r(t-\Delta t)}{2\Delta t}$$
A more in-depth description of the Verlet time integration scheme is given in "Understanding molecular simulation: from algorithms to applications (Frenkel, Daan, and Berend Smit., Elsevier, 2001.) "