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 system. Trajectory files containing a time series of atomic positions, velocities, and forces 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 draw 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 Newtons equation. A numerical approximation is employed, wherein one treats the velocities of the atoms as constant during a small 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 updated 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 atomic velocities are useful to estimate the kinetic energy of the system. 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.) "