Newtonian Mechanics and Numerical Integration

The Newtonian equations of motion can be expressed as

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

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}$$ (104)

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)$$ (105)

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. 105 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}$$ (106)

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}$$ (107)

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

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

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

This leads us to the comforting result that

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

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. The first algorithm considered in F&S is 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)... ...l{\dots}{r} + \mathscr{O}\left[\left(\Delta t\right)^4\right],$$ (111)

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

$$r\left(t-\Delta t\right) = r\left(t\right) - v\left(t\right)... ...l{\dots}{r} + \mathscr{O}\left[\left(\Delta t\right)^4\right].$$ (112)

When we add these together, we obtain Eq. 4.2.3 in F&S:

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

Eq. 113 is termed the ``Verlet'' algorithm (going back to Verlet's simulations of liquid argon [7]). 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$$ (114)

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-... ...}{2\Delta t} + \mathscr{O}\left[\left(\Delta t\right)^2\right]$$ (115)

Eq. 115 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$$ (116)

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.

Sec. 4.3.1 in F&S details a few other integration algorithms. Among them is the most popular integrator, the ``Velocity Verlet'' algorithm [8]. Every MD code I have every 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$$ (117)

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

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 standard 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 \mbox{Update positions}$$ (119)

$$v\left(t+\frac{1}{2}\Delta t\right) = v\left(t\right) + \frac{f\left(t\right)}{2m}\Delta t \mbox{Half-update velocities}$$ (120)

$$r\left(t+\Delta t\right) \rightarrow f\left(t+\Delta t\right) \mbox{Compute forces.}$$

$$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 \mbox{Half-update velocities}$$ (121)

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. This is an aspect not explicitly mentioned in F&S (at least at the point where periodic boundaries are introduced). 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 longer than a box size are excluded in systems with periodic boundaries because they cancel themselves.