Finite Difference Method

The finite difference method can generate MD trajectories with continuous potential models. The essential idea is that the integration is divided into many small steps, each separated by a fixed time interval $\delta t$. The total interactions on each particle at time $t$ can be calculated from the sum of interactions from other particles. The force is assumed to be constant during the time step $t$ and $t+\delta t$. The forces and accelerations of the particles in new positions can be determined, and so on.

There are many algorithms for integrating the equations of motion using finite difference methods, which are commonly used in MD simulations. All algorithms assume that the positions and dynamics properties (velocities, accelerations, etc.) can be approximated in Taylor series expansions:

$$r(r+\delta t)=r(t)+\delta tv(t)+\frac{1}{2}\delta t^{2}a(t)+\frac{1}{6}\delta t^{3}b(t)+\frac{1}{24}\delta t^{4}c(t)+...$$ (9)

$$v(t+\delta t)=v(t)+\delta ta(t)+\frac{1}{2}\delta t^{2}b(t)+\frac{1}{6}\delta t^{3}c(t)+...$$ (10)

$$a(t+\delta t)=a(t)+\delta b(t)+\frac{1}{2}\delta t^{2}c(t)...$$ (11)

where $v$ is the velocity (the first derivative of the position with respect to time), $a$ is the acceleration (the second derivative), $b$ is the third derivative, and so on. The Verlet algorithm [30] is the most broadly used method for integrating the trajectories of motion in MD simulations. The Verlet algorithm uses the positions and accelerations at time $t$ and the positions from the previous step, $r(t-\delta t)$, to calculate the new positions $r(t+\delta t)$. We can write down the following equations between these quantities and the velocities at time $t$:

$$r(t+\delta t)=r(t)+\delta tv(t)+\frac{1}{2}\delta t^{2}a(t)+...$$ (12)

$$r(t-\delta t)=r(t)-\delta tv(t)+\frac{1}{2}\delta t^{2}a(t)-....$$ (13)

Combining these two equation gives

$$r(t+\delta t)=2r(t)-r(t-\delta t)+\delta t^{2}a(t).$$ (14)

The velocities do not explicitly appear in the Verlet integration algorithm. The velocities can be calculated by dividing the difference in positions at time $t+\delta t$ and $t-\delta t$ by $2\delta t$ :

$$v(t)=[r(t+\delta t)-r(t-\delta t)]/2\delta t.$$ (15)

Alternatively, the velocities can be estimated at the half-step, $t+\frac{1}{2}\delta t$ :

$$v(t+\frac{1}{2}\delta t)=[r(t+\delta t)-r(t)]/\delta t.$$ (16)

The execution of the Verlet algorithm is straightforward and comprises $r(t)$, $r(t-\delta t)$, and the accelerations $a(t)$. One of its disadvantage is that the positions $r(t+\delta t)$ are obtained by adding a small term, $\delta t^{2}a(t)$, to the difference of two much larger terms, $2r(t)$ and $r(t-\delta t$). This may lead to a loss of precision. The lack of an explicit velocity term in the equations makes it difficult to obtain the velocities. The velocities are not available until the positions have been computed at the next step. In addition, it is not a self-starting algorithm; the new positions are obtained from the current positions $r(t)$ and the previous position $r(t-\delta t)$. At $t=0$ there is obviously only one set of positions and so it is necessary to employ some other means to obtain positions at $t-\delta t$. One way to obtain $r(t-\delta t)$ is to use the Taylor series, Eq. (9), truncated after the first term. Thus, $r(-\delta t)=r(0)-\delta tv(0)$.

The velocity Verlet method gives positions, velocities and accelerations at the same time and does not influence the precision:

$$r(t+\delta t)=r(t)+\delta t\, v(t)+\frac{1}{2}\delta t^{2}a(t)$$ (17)

$$v(t+\delta t)=v(t)+\frac{1}{2}\delta t[a(t)+a(t+\delta t)].$$ (18)

The velocity Verlet method [26] is actually implemented as a three-stage procedure because, as can be seen from Eq. (18), calculating the new velocities requires the accelerations at both $t$ and $t+\delta t$. Thus in the first step the positions at $t+\delta t$ are calculated according to Eq. (17) using the velocities and the accelerations at time $t$. The velocities at time $t+\frac{1}{2}\delta t$ are then determined using

$$v(t+\frac{1}{2}\delta t)=v(t)+\frac{1}{2}\delta t\, a(t).$$ (19)

New forces are computed from the current positions, giving $a(t+\delta t$) . In the final step, the velocities at time $t+\delta t$ are determined using

$$v(t+\delta t)=v(t+\frac{1}{2}\delta t)+\frac{1}{2}\delta ta(t+\delta t).$$ (20)