Introduction

Let us consider a few of the important elements of any MD program:

1. Force evaluation;
2. Integration; and
3. Configuration output.

We will understand these elements by manipulating an existing simulation program that implements the Lennard-Jones fluid (which you may recall was analyzed using Metropolis Monte-Carlo simulation in Sec. 3). The program is called mdlj.c in the instructional codes repository.

Recall the Lennard-Jones pair potential:

$$u\left(r\right) = 4\epsilon\left[\left(\frac{\sigma}{r}\right)^{12}- \left(\frac{\sigma}{r}\right)^{6}\right]$$ (150)

When implementing this in an MD code, similar to its implementation in MC, we adopt a reduced unit system, in which we measure length in $\sigma$, and energy in $\epsilon$. Additionally, for MD, we need to measure mass, and we adopt units of particle mass, $m$. This convention makes the force on a particle numerically equivalent to its acceleration. With these conventions, time is a derived unit:

$$t [=] \sigma\sqrt{m/\epsilon}$$ (151)

We also measure reduced temperature in units of $\epsilon/k_B$; so for a system of identical Lennard-Jones particles:

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

(Recall that the mass $m$ is 1 in Lennard-Jones reduced units.) These conventions obviate the need to perform unit conversions in the code.

We have already encountered interparticle forces for the Lennard-Jones pair potential in the context of computing the pressure from the virial in the MC simulation of the LJ fluid (Sec. 3.6). Briefly, the force exerted on particle $i$ by virtue of its Lennard-Jones interaction with particle $j$, ${\bf f}_{ij}$, is given by:

$${\bf f}_{ij}\left(r_{ij}\right) = \frac{{\bf r}_{ij}}{r_{ij}^2}\l... ...{1}{2}\left(\frac{\sigma}{r_{ij}}\right)^6\right]\right\} \equiv {\bf r}_{ij}f.$$ (153)

And, as shown in the previous section, once we have computed the vector ${\bf f}_{ij}$, we automatically have ${\bf f}_{ji}$, because ${\bf f}_{ij} = -{\bf f}_{ji}$. The scalar $f$ is called a “force factor.” If $f$ is negative, the force vector ${\bf f}_{ij}$ points from $i$ to $j$, meaning $i$ is attracted to $j$. Likewise, if $f$ is positive, the force vector ${\bf f}_{ij}$ points from $j$ to $i$, meaning that $i$ is being forced away from $j$.

Below is a C-code fragment for computing both the total potential energy and interparticle forces:

 double forces ( double * rx, double * ry, double * rz, 
 	             double * fx, double * fy, double * fz, 
 	             int n ) {
  int i,j;
  double dx, dy, dz, r2, r6, r12;
  double e = 0.0, f = 0.0;
 
  for (i=0;i<n;i++) {
    fx[i] = 0.0;
    fy[i] = 0.0;
    fz[i] = 0.0;
  }
  for (i=0;i<(n-1);i++) {
    for (j=i+1;j<n;j++) {
      dx     = rx[i]-rx[j];
      dy     = ry[i]-ry[j];
      dz     = rz[i]-rz[j];
      r2     = dx*dx + dy*dy + dz*dz;
      r6i    = 1.0/(r2*r2*r2);
      r12i   = r6i*r6i;
      e     += 4*(r12i - r6i);
      f      = 48/r2*(r6i*r6i-0.5*r6i);
      fx[i] += dx*f;
      fx[j] -= dx*f;
      fy[i] += dy*f;
      fy[j] -= dy*f;
      fz[i] += dz*f;
      fz[j] -= dz*f;
    }
  }
  return e;
}

Notice that the argument list now includes arrays for the forces, and because force is a vector quantity, we have three parallel arrays for a three-dimensional system. These forces must of course be initialized, shown in lines 6-10. The $N^2$ loop for visiting all unique pairs of particles is opened on lines 11-12. The inside of this loop is very similar to the evaluation of potential first presented in the MC simulation of the Lennard-Jones fluid; the only real difference is the computation of the “force factor,” $f$, on line 20, and the subsequent increment of force vector components on lines 21-26. Notice as well that there is no implementation of periodic boundary conditions in this code fragment; it was left out for simplicity. What would this “missing” code do? (Hint: look at the code mdlj.c for the answer.)

The second major aspect of MD is the integrator. As discussed in class, we will primarily use Verlet-style (explicit) integrators. The most common version is the velocity-Verlet algorithm [5], first presented in Sec. 4.1.1. Below is a fragment of C-code for executing one time step of integration for a system of $N$ particles:

for (i=0;i<N;i++) {
  rx[i]+=vx[i]*dt+0.5*dt2*fx[i];
  ry[i]+=vy[i]*dt+0.5*dt2*fy[i];
  rz[i]+=vz[i]*dt+0.5*dt2*fz[i];
  vx[i]+=0.5*dt*fx[i];
  vy[i]+=0.5*dt*fy[i];
  vz[i]+=0.5*dt*fz[i];
}

PE = total_e(rx,ry,rz,fx,fy,fz,N,L,rc2,ecor,ecut,&vir);
      
KE = 0.0;
for (i=0;i<N;i++) {
  vx[i]+=0.5*dt*fx[i];
  vy[i]+=0.5*dt*fy[i];
  vz[i]+=0.5*dt*fz[i];
  KE+=vx[i]*vx[i]+vy[i]*vy[i]+vz[i]*vz[i];
}
KE*=0.5;

Notice the update of positions (Eq. 125), where vx[i] is the $x$-component of velocity, fx[i] is the $x$-component of force, dt and dt2 are the time-step and squared time-step, respectively. Notice that there is no implementation of periodic boundaries in this code fragment; what would this “missing code” look like? (Hint: see mdlj.c for the answer!) Lines 5-7 are the first half-update of velocities (Eq. 126). The force routine computes the new forces on the currently updated configuration on line 10. Then, lines 12-18 perform the second-half of the velocity update (Eq. 127). Also note that the kinetic energy, $\mathscr{K}$, is computed in this loop.