Particle Trajectories

Note that the classical equation of motion

\[\vec{F}_i = m_i \frac{\partial \vec{v}_i}{\partial t} = m_i \frac{\partial ^2 r_i}{\partial t^2}\]

is a second order differential equation where the forces are assumed to be constant over the timestep. The usual way to solve this kind of differential equation numerically is by finite difference techniques, which involve the use of a fixed timestep \(\delta t\).

The general idea is to obtain the dynamic information at a later time \(t + \delta t\) to a significant degree of accuracy, based on the given atomic positions and velocities at time \(t\). This is repeatedly solved on a step-by-step basis and the solution is then propagated through time by using this approximation. The degree of accuracy can be obtained by expanding a Taylor series about \(\vec{r}(t)\) for an atom \(I\):

\[\vec{r}_i(t+\delta t) = \vec{r}_i(t) + \vec{v}_i(t)\delta t + \frac{1}{2}\vec{a}(t)\delta t^2 + ...\]

where \(\vec{a}\) is the acceleration, equal to \(\frac{\vec{f}}{m}\). The diagram below illustrates the general idea for the integration of motion. Essentially, we are interested in determining the net displacement of an atom over the timestep :math:delta t` due to the force experienced by the atom, to produce a new set of dynamic information.

One of the simplest and fastest methods to compute this integration is the leapfrog Verlet integration scheme. The scheme requires values of position and force at time \(t\) but the velocity at half a timestep behind, \(t - \frac{\delta t}{2}\).

First of all, the force at time \(t\) is computed. Then the velocity is advanced to \(t + \frac{\delta t}{2}\) using the kinematic equation v = u + ft/m, where m is the mass of the atom:

\[\vec{v}_i \left(t+\frac{\delta t}{2}\right) = \vec{v}_i \left(t-\frac{\delta t}{2}\right) + \frac{\vec{f}(\vec{r}_i(t)) \delta t}{m}\]

where the initial velocity (half a step behind) is obtained as

\[\vec{v}_i \left(t-\frac{\delta t}{2}\right) = \frac{\vec{r}_i(t) - \vec{r}_i(t-\delta t)}{\delta t}\]

After that, the new position is advanced, or ‘leaps forward’, a full step ahead to \(t + \delta t\) using the half-step velocity calculated above:

\[\vec{r}_i\left(t+\delta t\right) = \vec{r}_i\left(t\right) + \vec{v}_i \left(t+\frac{\delta t}{2}\right)\delta t\]

However, the current velocity at time \(t\) would be needed to calculate properties such as the kinetic energy and pressure:

\[\vec{v}_i(t) = \frac{1}{2} \left[\vec{v}_i \left(t-\frac{\delta t}{2}\right) + \vec{v}_i\left(t+\frac{\delta t}{2}\right)\right]\]

After that, the force at the new position will be calculated and the whole integration process is repeated. The diagram below shows a schematic illustration of the integration algorithm.

Note that, although leapfrog Verlet is simple and faster than most other schemes, it is not time reversible and can become unstable over long time scales.

The use of a symplectic integrator provides time reversibility and long-term stability. One commonly-used integration scheme for MD simulations is Velocity Verlet (VV), which is second-order accurate (i.e. errors are proportional to the square of the timestep size \(\delta t\)) and updates particle positions and velocities without requiring additional memory to do so.

Starting with positions, velocities and forces all at time \(t\), the velocity for atom \(i\) is advanced to \(t + \frac{\delta t}{2}\) using the force:

\[\vec{v}_i \left(t+\frac{\delta t}{2}\right) = \vec{v}_i \left(t\right) + \frac{\vec{f}(\vec{r}_i(t)) \delta t}{2 m}\]

before its position is advanced a full step ahead to \(t + \delta t\) using the half-step velocity calculated above:

\[\vec{r}_i\left(t+\delta t\right) = \vec{r}_i\left(t\right) + \vec{v}_i \left(t+\frac{\delta t}{2}\right)\delta t\]

The force at the new position will be calculated, which is then used to advance the velocity by another half step:

\[\vec{v}_i \left(t+\delta t\right) = \vec{v}_i \left(t+\frac{\delta t}{2}\right) + \frac{\vec{f}(\vec{r}_i(t+\delta t)) \delta t}{2 m}\]

While this integration algorithm requires two passes through the atoms per timestep, it keeps the velocities synchronised with the positions and forces, avoiding the need to store previous values for property calculations.

Note

Older versions of DL_POLY offer VV by default, while only VV is available in DL_POLY_5. It generates trajectories in the microcanonical (NVE) ensemble in which the total energy is conserved. If this drifts or fluctuates excessively in the course of a simulation, this may mean the timestep is too large or the potential cutoffs are too small.

^ GO TO TOP ^