Newtonian Mechanics and Numerical Integration

The Newtonian equations of motion can be expressed as

$$m \ddot{\bf r}_i + \nabla_i\mathscr{U} = 0$$ (109)

where $\ddot{\bf r}_i$ is the acceleration of particle $i$, and the force acting on particle $i$ is given by the negative gradient of the total potential, $\mathscr{U}$, with respect to its position:

$${\bf f}_i = -\nabla_i\mathscr{U} = -\frac{\partial\mathscr{U}}{\partial{\bf r}_i}$$ (110)

Whereas in a typical MC simulation, in which all we really need is the ability to evaluate the potential energy of a configuration, in MD we actually need to evaluate all interparticle forces for a configuration.

We first encountered interparticle forces in Sec. 3.6 in a discussion of the virial in computing pressure in a standard Metropolis Monte Carlo simulation of the Lennard-Jones liquid. At this point, it suffices to consider a system with generic pairwise interactions, for which the total potential is given by:

$$\mathscr{U} = \sum_j\sum_{k>j} u_{jk}\left(r_{jk}\right)$$ (111)

where $r_{jk}$ is the scalar distance between particles $j$ and $k$, and $u_{jk}$ is the pair potential specific to pair $(j,k)$. For a system of $N$ identical particles, Eq. 111 is a summation of $\frac{1}{2}N\left(N-1\right)$ terms. So, the force on any particular particle, $i$, selects $N$ terms from the above summation; that is, those terms involving particle $i$:

$${\bf f}_i = - \sum_{j=1}^{N} \frac{\partial u_{ij}\left(r_{ij}\right)}{\partial {\bf r}_i} = \sum_{j=1}^{N} {\bf f}_{ij}$$ (112)

where we can define the quantity ${\bf f}_{ij}$ is the force exerted on particle $i$ by virtue of the fact that it interacts with particle $j$. Because $u_{ij}$ is a function of a scalar quantity, we can break the derivative up:

$${\bf f}_{ij}\left(r_{ij}\right) = -\frac{\partial U_{LJ}\left(r_{ij}\right)}{\partial {\bf r}_i}$$ (113)

$$= -\frac{{\bf r}_{ij}}{r_{ij}}\frac{\partial U_{LJ}\left(r_{ij}\right)}{\partial r_{ij}}$$ (114)

Eq. 114 illustrates that, because ${\bf r}_{ij} = -{\bf r}_{ji}$,

$${\bf f}_{ij} = -{\bf f}_{ji}$$ (115)

This leads us to the comforting result that

$${\bf F} = \sum_i {\bf f}_i = 0$$ (116)

That is, the total force on the collection of particles is zero. (The same result holds identically for all potentials which are functions of relative interatomic positions only.) But the practical advantage of this result is that, when we visit the pair $(i,j)$ and compute the force on $i$ due to its interaction with $j$, ${\bf f}_{ij}$, we automatically have the force on $j$ due to its interaction with $i$, $-{\bf f}_{ij}$. Some refer to this as “Newton's Third Law.”

The other key aspect of a simple MD program is a means of numerical integration of the equations of motion of each particle. We first consider the “simple Verlet” algorithm, which is an explicit integration scheme. Let us consider a Taylor-expanded version of one coordinate of the position of a particular particle, $r(t)$:

$$r\left(t+\Delta t\right) = r\left(t\right) + v\left(t\right)\Delt... ...t)^3}{3!}\stackrel{\dots}{r} + \mathscr{O}\left[\left(\Delta t\right)^4\right],$$ (117)

and, letting $\Delta t \rightarrow -\Delta t$,

$$r\left(t-\Delta t\right) = r\left(t\right) - v\left(t\right)\Delt... ...t)^3}{3!}\stackrel{\dots}{r} + \mathscr{O}\left[\left(\Delta t\right)^4\right].$$ (118)

When we add these together, we obtain:

$$r\left(t+\Delta t\right) \approx 2r\left(t\right) - r\left(t-\Delta t\right) + \frac{f\left(t\right)}{2m}\left(\Delta t\right)^2.$$ (119)

Eq. 119 is termed the “Verlet” algorithm (going back to Verlet's simulations of liquid argon [4]). Notice that, when one chooses a small $\Delta t$, one can predict the position of a particle at time $t + \Delta t$ given its position at time $t$ and the force acting on it at time $t$. We see that the new position coordinate has an error of order $\left(\Delta t\right)^4$. $\Delta t$ is called the “time step” in a molecular dynamics simulation.

A system obeying Newtonian mechanics conserves total energy. For a dynamical system (i.e., a system of interacting particles) obeying Newtonian mechanics, the configurations generated by integration are members of the microcanonical ensemble; that is, the ensemble of configurations for which $NVE$ is constant, constrained to a subvolume $\Omega$ in phase space. The “natural” ensemble for Metropolis Monte Carlo, you will recall, is canonical; for MD, it is microcanonical. Later, we will consider techniques for conducting MD simulations in other ensembles (at constant temperature and/or pressure, for example).

When the Verlet algorithm is used to integrate Newtonian equations of motion, the total energy of the system is conserved to within a finite error, so long as $\Delta t$ is “small enough.” How does one determine a reasonable value for $\Delta t$? Basically the same way we determined reasonable maximum displacements in continuous-space MC simulation: trial and error. We will play with time-step values in the next section, in which we consider MD simulation of the Lennard-Jones liquid.

In saying that the total energy is conserved, we realize that total energy is the sum of potential and kinetic energy. To integrate the equations of motion, we need to compute neither the potential or kinetic energy, so we have to take extra steps in an MD program to make sure total energy is being conserved. Potential energy is easily accumulated during the calculation of forces, but kinetic energy has to be computed using particle velocities:

$$\mathscr{K} = \frac{1}{2}\sum_i m_i\left\vert{\bf v}_i\right\vert^2$$ (120)

But where are velocities in the Verlet algorithm? They are not necessary for updating positions, but can be easily “generated” provided one stores previous, current, and next-time-step positions in implementing the algorithm:

$$v\left(t\right) = \frac{r\left(t+\Delta t\right) - r\left(t-\Delta t\right)}{2\Delta t} + \mathscr{O}\left[\left(\Delta t\right)^2\right]$$ (121)

Eq. 121 is used in Algorithm 6 in F&S (Integration of the Equations of Motion) simply in order to compute $\mathscr{K}$. Energy consservation can be checked by tracking the total energy, $\mathscr{U}+\mathscr{K}$.

While we are considering the instantaneous kinetic energy, $\mathscr{K}$, it is useful to recognize a working definition of instantaneous temperature, $T$:

$$\frac{3}{2}Nk_BT = \mathscr{K} = \frac{1}{2}\sum_{i=1}^{N}m_i\left\vert{\bf v}_i\right\vert^2$$ (122)

Because $\mathscr{K}$ fluctuates in time as the system evolves, so does the temperature. So, the actual temperature of a system in a microcanonical MD simulation is a time-average.

Perhaps the most popular integrator is the “velocity Verlet” algorithm [5]. Every MD code I have ever written or used (this totals a dozen or so) has used the velocity Verlet algorithm, so I feel at least it is worth explaining here. The velocity Verlet algorithm requires updates of both positions and velocities:

$$r\left(t+\Delta t\right) = r\left(t\right) + v\left(t\right)\Delta t + \frac{f\left(t\right)}{2m}\left(\Delta t\right)^2$$ (123)

$$v\left(t+\Delta t\right) = v\left(t\right) + \frac{f\left(t+\Delta t\right)+f\left(t\right)}{2m}\Delta t$$ (124)

The update of velocities uses an arithmetic average of the force at time $t$ and $t + \Delta t$. This results in a slightly more stable integrator compared to the simple Verlet algorithm, in that one may use slightly larger time-steps to achieve the same level of energy conservation. (Aside: a nice project idea is to quantify this statement.) This might imply that one has to maintain two parallel force arrays. In practice, this is not necessary, because the velocity update can be split to either side of the force computation, forming a so-called “leapfrog” algorithm:

$$r\left(t+\Delta t\right) = r\left(t\right) + v\left(t\right)\Delta t + \frac{f\left(t\right)}{2m}\left(\Delta t\right)^2$$ (125)

$$v\left(t+\frac{1}{2}\Delta t\right) = v\left(t\right) + \frac{f\left(t\right)}{2m}\Delta t$$ (126)

$$r\left(t+\Delta t\right) \rightarrow f\left(t+\Delta t\right)$$

$$v\left(t+\Delta t\right) = v\left(t+\frac{1}{2}\Delta t\right) + \frac{f\left(t+\Delta t\right)}{2m}\Delta t$$ (127)

In the next section (Sec. 4.2), we will consider the velocity Verlet algorithm in the context of an MD simulation of the Lennard-Jones fluid.

As a final tidbit, we must consider periodic boundaries applied in a molecular dynamics simulation. Consider modes of a system. Think of a mode as a concerted vibration of collections of particles with a characteristic wavelength. A dense system will have short wavelength (local) modes, and long wavelength modes, like large-scale concerted “sloshing” of the particles in the system. These modes exist naturally in matter, and the partitioning of energy among these various modes is important to understand in describing some transport properties. The key caveat is that modes with wavelengths that are incommensurate with the box size are not permitted in a periodic system because they cancel themselves. A mode is commensurate with the box so long as an integer multiple of its wavelength is the box length. This can actually be very restrictive in systems with a wide span of wavelengths, like amorphous unstructured solids, but is not that important for amorphous liquids.