The Newtonian equations of motion can be expressed as \begin{equation} m \ddot {\bf r}_i + \nabla _i\mathscr {U} = 0 \end{equation} 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: \begin{equation} {\bf f}_i = -\nabla _i\mathscr {U} = -\frac {\partial \mathscr {U}}{\partial {\bf r}_i} \end{equation} 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: \begin{equation} \label {eq:md:tot} \mathscr {U} = \sum _j\sum _{k>j} u_{jk}\left (r_{jk}\right ) \end{equation} 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. 92 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\): \begin{equation} \label {eq:md:frci} {\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} \end{equation} 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: \begin{equation}\label {eq:md:fij} \begin {aligned} {\bf f}_{ij}\left (r_{ij}\right ) &= -\frac {\partial U_{LJ}\left (r_{ij}\right )}{\partial {\bf r}_i} \\ &= -\frac {{\bf r}_{ij}}{r_{ij}}\frac {\partial U_{LJ}\left (r_{ij}\right )}{\partial r_{ij}} \end {aligned} \end{equation} Eq. 94 illustrates that, because \({\bf r}_{ij} = -{\bf r}_{ji}\), \begin{equation} {\bf f}_{ij} = -{\bf f}_{ji} \end{equation} This leads us to the comforting result that \begin{equation} {\bf F} = \sum _i {\bf f}_i = 0 \end{equation} 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)\): \begin{equation} 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 + \frac {\left (\Delta t\right )^3}{3!}\stackrel {\dots }{r} + \mathscr {O}\left [\left (\Delta t\right )^4\right ], \end{equation} and, letting \(\Delta t \rightarrow -\Delta t\), \begin{equation} 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 - \frac {\left (\Delta t\right )^3}{3!}\stackrel {\dots }{r} + \mathscr {O}\left [\left (\Delta t\right )^4\right ]. \end{equation} When we add these together, we obtain Eq. 4.2.3 in F&S: \begin{equation} \label {eq:verlet1} 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. \end{equation} Eq. 99 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: \begin{equation} \mathscr {K} = \frac {1}{2}\sum _i m_i\left |{\bf v}_i\right |^2 \end{equation} 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: \begin{equation} \label {eq:verlet_vel} 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 ] \end{equation} Eq. 101 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\): \begin{equation} \frac {3}{2}Nk_BT = \mathscr {K} = \frac {1}{2}\sum _{i=1}^{N}m_i\left |{\bf v}_i\right |^2 \end{equation} 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: \begin{equation} \begin {aligned} 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 \\ v\left (t+\Delta t\right ) & = v\left (t\right ) + \frac {f\left (t+\Delta t\right )+f\left (t\right )}{2m}\Delta t \end {aligned} \end{equation} 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: \begin{equation}\label {eq:vv:pos} 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 \quad \text {(Update positions)} \end{equation} \begin{equation}\label {eq:vv:first_half} v\left (t+\frac {1}{2}\Delta t\right ) = v\left (t\right ) + \frac {f\left (t\right )}{2m}\Delta t \quad \text {(Half-update velocities)} \end{equation} \[ r\left (t+\Delta t\right ) \rightarrow f\left (t+\Delta t\right ) \quad \text {(Compute forces)} \] \begin{equation}\label {eq:vv:second_half} 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 \quad \text {(Half-update velocities)} \end{equation}
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.
In this section, we present an elegant formalism for deriving MD integrators, as discussed by Tuckerman et al. [9]. What we present here is essentially the first two parts of the second section of Reference [9], including some of my own elaboration and some of that presented in section 4.3 of F&S.
Imagine a quantity \(f\) which is a function of particle positions \({\bf r}^N\) and momenta \({\bf p}^N\). Its time derivative is given by \begin{equation} \label {eq:liouville_1} \dot {f} = \dot {\bf r}\frac {\partial f}{\partial {\bf r}} + \dot {\bf p}\frac {\partial f}{\partial {\bf p}} \end{equation} We can write down a formal solution to this equation. First, define the Liouville operator as \begin{equation} iL = \dot {\bf r}\frac {\partial }{\partial {\bf r}} + \dot {\bf p}\frac {\partial }{\partial {\bf p}} \end{equation} As Tuckerman points out in his coursenotes, the \(i\) is there by convention and ensures that the operator is Hermitian. We can re-express Eq. 107 as \begin{equation} \dot {f} = iLf \end{equation} which we solve directly to yield \begin{equation} f\left (t\right ) = \exp \left (iLt\right )f\left (0\right ). \end{equation} If \(f\) is itself a vector quantity identical to the set of positions and momenta, \(\Gamma \), we have a way to express, formally, the evolution of the system: \begin{equation} \label {eq:propagator} \Gamma \left (t\right ) = U\left (t\right )\Gamma \left (0\right ) \end{equation} where \(U(t) = \exp \left (iLt\right )\) is the classical propagator. The idea with numerical integration is that we find a way to represent the propagator as a discrete algorithm for constructing the system at some time \(t+ \Delta t\) given the system at time \(t\).
Let’s build our discrete integrator by decomposing the operator: \begin{equation} iL = iL_1 + iL_2 \end{equation} This does not necessarily lead to two independent propagators, because the two components do not commute; that is: \begin{equation} \exp \left [\left ( iL_1 + iL_2\right )t\right ] \ne \exp \left (iL_1t\right )\exp \left (iL_2t\right ) \end{equation}
Consider the action of the partial Liouville operator \begin{equation} iL_1 \equiv \dot {\bf r}\left (0\right )\frac {\partial }{\partial {\bf r}}, \end{equation} which gives \begin{equation} \begin {aligned} f\left [{\bf p}^N\left (t\right ),{\bf r}^N\left (t\right )\right ] & = \exp \left [t\dot {\bf r}\left (0\right )\frac {\partial }{\partial {\bf r}}\right ] f\left [{\bf p}^N\left (0\right ),{\bf r}^N\left (0\right )\right ] \\ & = \sum _{n=0}^{\infty } \frac {\left (\dot {\bf r}\left (0\right )t\right )^n}{n!}\frac {\partial ^n}{\partial {\bf r}^n}f\left [{\bf p}^N\left (0\right ),{\bf r}^N\left (0\right )\right ] \\ & = f\left [{\bf p}^N\left (0\right ),{\bf r}^N\left (0\right )+\dot {\bf r}\left (0\right )\right ] \end {aligned} \end{equation} The last line is the collapse of the Taylor expansion of the line immediately above it. So, the effect of this operator fragment is a simple shift of coordinates given some initial velocities. This is an interesting fact: we can consider first-order integration as a Taylor expansion.
The next step of Tuckerman was to apply the Trotter identity: \begin{equation} \exp \left [\left ( iL_1 + iL_2\right )t\right ] = \lim _{P\rightarrow \infty } \left [\exp \left (iL_1t/2P\right )\exp \left (iL_t2/P\right )\exp \left (iL_1t/2P\right )\right ]^P \end{equation} When \(P\) is large but finite: \begin{equation} \exp \left [\left (iL_1 + iL_2\right )t\right ] = \left [\exp \left (iL_1t/2P\right )\exp \left (iL_2t/P\right )\exp \left (iL_1t/2P\right )\right ]^P\exp \left [\mathscr {O}\left (1/P^2\right )\right ] \end{equation} Now, we define a finite timestep as \(\Delta t = t/P\) and we have then a discrete operator that, when applied to a configuration at time \(t\), will produce the configuration at time \(t + \Delta t\): \begin{equation} \exp \left (iL_1\Delta t/2\right )\exp \left (iL_2\Delta t\right )\exp \left (iL_1\Delta t/2\right )\Gamma \left (t\right ) = \Gamma \left (t+\Delta t\right ) \end{equation} By performing this operation sequentially \(P\) times, we recover a discretized version of the formal solution to generate \(\Gamma \left (t\right )\) given \(\Gamma \left (0\right )\).
Now we explicitly consider the decomposition: \begin{equation} \begin {aligned} iL_1 & = \dot {\bf p}\left (0\right )\frac {\partial }{\partial {\bf p}} \\ iL_2 & = \dot {\bf r}\left (0\right )\frac {\partial }{\partial {\bf r}} \end {aligned} \end{equation}
We can perform one \(\Delta t\)’s worth of update using the following operation on \(f\): \[ \exp \left (\dot {\bf p}\left (0\right )\left (\frac {\partial }{\partial {\bf p}}\right )\frac {\Delta t}{2}\right ) \exp \left (\dot {\bf r}\left (0\right )\left (\frac {\partial }{\partial {\bf r}}\right )\Delta t\right ) \exp \left (\dot {\bf p}\left (0\right )\left (\frac {\partial }{\partial {\bf p}}\right )\frac {\Delta t}{2}\right ) f\left [{\bf p}^N\left (t\right ),{\bf r}^N\left (t\right )\right ] \] The action of the rightmost operator, \(\exp \left (\dot {\bf p}\left (0\right )\left (\frac {\partial }{\partial {\bf p}}\right )\frac {\Delta t}{2}\right )\): \[ f\left [{\bf p}^N\left (t\right ),{\bf r}^N\left (t\right )\right ] \rightarrow f\left \{\left [{\bf p}\left (t\right )+\frac {\Delta t}{2}\dot {\bf p}\left (\Delta t\right )\right ]^N,\left [{\bf r}\left (t\right )\right ]^N\right \} \] The action of the next rightmost, \(\exp \left (\dot {\bf r}\left (0\right )\left (\frac {\partial }{\partial {\bf r}}\right )\Delta t\right )\): \[ f\left \{\left [{\bf p}\left (t\right )+\frac {\Delta t}{2}\dot {\bf p}\left (\Delta t\right )\right ]^N,\left [{\bf r}\left (t\right )\right ]^N\right \} \rightarrow f\left \{\left [{\bf p}\left (t\right )+\frac {\Delta t}{2}\dot {\bf p}\left (0\right )\right ]^N, \left [{\bf r}\left (t\right )+ \Delta t \dot {\bf r}\left (\frac {\Delta t}{2}\right )\right ]^N\right \} \] Then, the action of the final operator: \[ f\left \{\left [{\bf p}\left (t\right )+\frac {\Delta t}{2}\dot {\bf p}\left (0\right )\right ]^N, \left [{\bf r}\left (t\right )+ \Delta t \dot {\bf r}\left (\frac {\Delta t}{2}\right )\right ]^N\right \} \rightarrow f\left \{\left [{\bf p}\left (t\right )+\frac {\Delta t}{2}\dot {\bf p}\left (0\right )+ \frac {\Delta t}{2}\dot {\bf p}\left (\Delta t\right )\right ]^N, \left [{\bf r}\left (t\right )+ \Delta t \dot {\bf r}\left (\frac {\Delta t}{2}\right )\right ]^N\right \}\\ \]
Noting that \({\bf p} = m\dot {\bf r}\) and \({\bf F} = m{\bf a} = \dot {\bf p}\), we can summarize the effect of this three-step update of the positions and velocities as \begin{equation} \begin {aligned} {\bf r}\left (\Delta t\right ) & = {\bf r}\left (0\right ) + \Delta t \dot {\bf r}\left (0\right ) + \frac {\left (\Delta t\right )^2}{2} \frac {{\bf F}\left [{\bf r}\left (0\right )\right ]}{m}, \\ \dot {\bf r}\left (\Delta t\right ) & = \dot {\bf r}\left (0\right ) + \frac {\Delta t}{2m} \left \{{\bf F}\left [{\bf r}\left (0\right )\right ] + {\bf F}\left [{\bf r}\left (\Delta t\right )\right ]\right \} \end {aligned} \end{equation}
This is the velocity-Verlet algorithm, seen previously in Eqs 104-106. Interestingly, we can also reverse the order of the decomposition; i.e., \begin{equation} \begin {aligned} iL_1 & = \dot {\bf r}\left (0\right )\frac {\partial }{\partial {\bf r}} \\ iL_2 & = \dot {\bf p}\left (0\right )\frac {\partial }{\partial {\bf p}} \end {aligned} \end{equation} The update algorithm that arises is \begin{equation} \begin {aligned} \dot {\bf r}\left (\Delta t\right ) & = \dot {\bf r}\left (0\right ) + \Delta t {\bf F}\left [{\bf r}\left (0\right ) + \frac {\Delta t}{2m}\dot {\bf r}\left (0\right )\right ] \\ {\bf r}\left (\Delta t\right ) & = {\bf r}\left (0\right ) + \frac {\Delta t}{2}\left [\dot {\bf x}\left (0\right ) +\dot {\bf x}\left (\Delta t\right )\right ]. \end {aligned} \end{equation} This is termed the position Verlet algorithm [9]. Tuckerman et al. showed that this new algorithm results in a slightly lower drift in total energy in MD simulation of a simple Lennard-Jones fluid, when the time-step is greater than about 0.004.