4.2 Case Study 1: An MD Code for the Lennard-Jones Fluid
4.2.1 Introduction
Let us consider a few of the important elements of any MD program:
- Force evaluation;
- Integration; and
- 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). Recall the Lennard-Jones pair potential: \begin{equation} \label {eq:md:lj} u\left (r\right ) = 4\epsilon \left [\left (\frac {\sigma }{r}\right )^{12}- \left (\frac {\sigma }{r}\right )^{6}\right ] \end{equation} 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: \begin{equation} t [=] \sigma \sqrt {m/\epsilon } \end{equation} We also measure reduced temperature in units of \(\epsilon /k_B\); so for a system of identical Lennard-Jones particles: \begin{equation} \label {eq:md:T} \frac {3}{2}NT = \mathscr {K} = \frac {1}{2}\sum _{i=1}^{N}\left |{\bf v}_i\right |^2 \end{equation} (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: \begin{equation} {\bf f}_{ij}\left (r_{ij}\right ) = \frac {{\bf r}_{ij}}{r_{ij}^2}\left \{48\epsilon \left [\left (\frac {\sigma }{r_{ij}}\right )^{12} - \frac {1}{2}\left (\frac {\sigma }{r_{ij}}\right )^6\right ]\right \} \equiv {\bf r}_{ij}f. \end{equation} 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:
0 double forces ( double * rx, double * ry, double * rz,
1 double * fx, double * fy, double * fz, int n ) {
2 int i,j;
3 double dx, dy, dz, r2, r6, r12;
4 double e = 0.0, f = 0.0;
5
6 for (i=0;i<n;i++) {
7 fx[i] = 0.0;
8 fy[i] = 0.0;
9 fz[i] = 0.0;
10 }
11 for (i=0;i<(n-1);i++) {
12 for (j=i+1;j<n;j++) {
13 dx = (rx[i]-rx[j]);
14 dy = (ry[i]-ry[j]);
15 dz = (rz[i]-rz[j]);
16 r2 = dx*dx + dy*dy + dz*dz;
17 r6i = 1.0/(r2*r2*r2);
18 r12i = r6i*r6i;
19 e += 4*(r12i - r6i);
20 f = 48/r2*(r6i*r6i-0.5*r6i);
21 fx[i] += dx*f;
22 fx[j] -= dx*f;
23 fy[i] += dy*f;
24 fy[j] -= dy*f;
25 fz[i] += dz*f;
26 fz[j] -= dz*f;
27 }
28 }
29 return e;
30 }
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 incrementation 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 [8], 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:
1 for (i=0;i<N;i++) {
2 rx[i]+=vx[i]*dt+0.5*dt2*fx[i];
3 ry[i]+=vy[i]*dt+0.5*dt2*fy[i];
4 rz[i]+=vz[i]*dt+0.5*dt2*fz[i];
5 vx[i]+=0.5*dt*fx[i];
6 vy[i]+=0.5*dt*fy[i];
7 vz[i]+=0.5*dt*fz[i];
8 }
9
10 PE = total_e(rx,ry,rz,fx,fy,fz,N,L,rc2,ecor,ecut,&vir);
11
12 KE = 0.0;
13 for (i=0;i<N;i++) {
14 vx[i]+=0.5*dt*fx[i];
15 vy[i]+=0.5*dt*fy[i];
16 vz[i]+=0.5*dt*fz[i];
17 KE+=vx[i]*vx[i]+vy[i]*vy[i]+vz[i]*vz[i];
18 }
19 KE*=0.5;
You will notice the update of positions in lines 1-8 (Eq. 104), 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. 105). 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. 106). Also note that the kinetic energy, \(\mathscr {K}\), is computed in this loop.
4.2.2 The Code, mdlj.c
The complete C program, mdlj.c, contains a complete implementation of the Lennard-Jones force routine and the velocity-Verlet integrator. Compilation instructions appear in the header comments. Let us now consider some sample results from mdlj.c. First, mdlj.c includes an option that outputs a brief summary of the command line options available:
cfa@abrams01:/home/cfa>mdlj -h mdlj usage: mdlj [options] Options: -N [integer] Number of particles -rho [real] Number density -dt [real] Time step -rc [real] Cutoff radius -ns [real] Number of integration steps -so Short-form output (unused) -T0 [real] Initial temperature -fs [integer] Sample frequency -sf [a|w] Append or write config output file -icf [string] Initial configuration file -seed [integer] Random number generator seed -h Print this info. cfa@abrams01:/home/cfa>
Let us run mdlj.c for 512 particles and 1000 time-steps at a density of 0.85 and an initial temperature of 2.5. We will pick a relatively conservative (small) time-step of 0.001. We will not specify an input configuration, instead allowing the code to create initial positions on a cubic lattice. Here is what we see in the terminal:
cfa@abrams01:/home/cfa/dxu/che800-002/md1/T0=2.0>../../mdlj -N 512 -fs 10 -ns 1000 -sf w -T0 2.5 -rho 0.85 -rc 2.5 # NVE MD Simulation of a Lennard-Jones fluid # L = 8.44534; rho = 0.85000; N = 512; rc = 2.50000 # nSteps 1000, seed 23410981, dt 0.00100 # step PE KE TE drift T P 0 -2662.98197 1919.36147 -743.62050 3.52596e-07 2.49917 4.28921 1 -2661.06414 1917.44274 -743.62140 1.56869e-06 2.49667 4.30335 2 -2657.85692 1914.23397 -743.62296 3.66100e-06 2.49249 4.32697 3 -2653.34343 1909.71825 -743.62517 6.64471e-06 2.48661 4.36017 4 -2647.49960 1903.87154 -743.62807 1.05355e-05 2.47900 4.40307 5 -2640.29423 1896.66259 -743.63164 1.53466e-05 2.46961 4.45584 6 -2631.68894 1888.05303 -743.63591 2.10836e-05 2.45840 4.51868 7 -2621.63837 1877.99751 -743.64086 2.77377e-05 2.44531 4.59184 8 -2610.09114 1866.44448 -743.64666 3.55346e-05 2.43027 4.67548 9 -2596.98921 1853.33668 -743.65253 4.34344e-05 2.41320 4.76991 10 -2582.27118 1838.61192 -743.65927 5.24918e-05 2.39403 4.87544 11 -2565.87073 1822.20481 -743.66592 6.14413e-05 2.37266 4.99241 12 -2547.72317 1804.05002 -743.67316 7.11682e-05 2.34902 5.12119 ^C
Each line of output after the header information corresponds to one time-step. The first column is the time-step, the second the potential energy, the third the kinetic energy, the fourth the total energy, the fifth the “drift,” the sixth the instantaneous temperature (Eq. 125), and seventh the instantaneous pressure (Eq. 82).
The drift is output to assess the stability of the explicit integration. As a rule of thumb, we would like to keep the drift to below 0.01% of the total energy. The drift reported by mdlj.c is computed as \begin{equation}\label {eq:md:edrift} \Delta \mathscr {T}\left (t\right ) = \frac {\mathscr {T}\left (t\right )-\mathscr {T}\left (0\right )}{\mathscr {T}\left (0\right )} \end{equation} where \(\mathscr {T}\) is the total energy. The plots below show the output trace for the full 1,000 time-steps.
Now, this invokation of mdlj.c produces a series of configuration files. Look at the listing of the directory in which the program is run:
cfa@abrams01:/home/cfa/dxu/che800-002/md1/T0=2.0>ls *.xyz 0.xyz 20.xyz 300.xyz 410.xyz 520.xyz 630.xyz 740.xyz 850.xyz 960.xyz 10.xyz 200.xyz 310.xyz 420.xyz 530.xyz 640.xyz 750.xyz 860.xyz 970.xyz 100.xyz 210.xyz 320.xyz 430.xyz 540.xyz 650.xyz 760.xyz 870.xyz 980.xyz 110.xyz 220.xyz 330.xyz 440.xyz 550.xyz 660.xyz 770.xyz 880.xyz 990.xyz 120.xyz 230.xyz 340.xyz 450.xyz 560.xyz 670.xyz 780.xyz 890.xyz 130.xyz 240.xyz 350.xyz 460.xyz 570.xyz 680.xyz 790.xyz 90.xyz 140.xyz 250.xyz 360.xyz 470.xyz 580.xyz 690.xyz 80.xyz 900.xyz 150.xyz 260.xyz 370.xyz 480.xyz 590.xyz 70.xyz 800.xyz 910.xyz 160.xyz 270.xyz 380.xyz 490.xyz 60.xyz 700.xyz 810.xyz 920.xyz 170.xyz 280.xyz 390.xyz 50.xyz 600.xyz 710.xyz 820.xyz 930.xyz 180.xyz 290.xyz 40.xyz 500.xyz 610.xyz 720.xyz 830.xyz 940.xyz 190.xyz 30.xyz 400.xyz 510.xyz 620.xyz 730.xyz 840.xyz 950.xyz cfa@abrams01:/home/cfa/dxu/che800-002/md1/T0=2.0>
Each one of these files is a configuration snapshot, and contains the \((x,y,z)\) position of each particle. The format of each data file is special: it is called the “XYZ” format (this is a standard format for many simulation programs). The filenames indicate which snapshot a file is; the number is the time-step value. Look inside one of the files, 690.xyz:
cfa@abrams01:/home/cfa/dxu/che800-002/md1/T0=2.0>more 690.xyz 512 1 16 0.29604446 0.56828714 0.37752815 -1.30786300 1.84221403 -0.89413789 16 1.21401795 8.19177706 0.33515333 3.52032041 -0.81317122 1.23241312 16 2.64726066 0.82254281 0.26994972 -0.50630797 0.93443275 1.81549173 16 3.75202252 0.06731821 0.39268752 1.14238124 0.02215183 0.43795406 16 4.62876356 0.70002120 0.97873282 -1.09076342 0.93155434 0.67508513 16 5.32461628 0.72950028 0.14052347 1.70998632 -0.86727858 0.03644588 16 6.18785137 1.40475117 8.26593130 1.25736693 0.75750213 -1.78288755 16 7.56095527 0.53604317 0.71204718 2.35604054 0.69772173 0.63651689 16 0.84782080 1.61955526 8.37085344 0.73714470 -0.05463901 0.04286659 16 1.68195240 1.15978775 0.40224323 -2.07249584 -0.98665932 -0.17256677 16 2.90316975 1.71194098 0.82270713 -1.44619783 2.35862108 -0.62228759 16 3.64445686 1.16196197 0.92857556 -1.17480938 -1.88716489 -0.75592199 16 4.63278427 1.47318092 0.01705540 1.40835198 -0.48737538 -0.29628750 16 6.14017777 2.46476252 0.04160396 0.31157778 3.21760722 1.10011638 16 7.36957085 1.92449651 0.87918692 -0.13392746 -0.59350997 0.47063528 16 8.19214702 1.21748389 1.23876905 1.17032561 -0.99116310 -3.36639536 16 1.31101056 2.32117258 0.60909154 0.23149412 1.08341509 1.35597275 16 2.08284899 2.56390254 0.07651616 2.75209335 0.31391146 -0.48580440 16 3.05603350 2.77944216 0.43924763 -2.05687178 -1.56404114 0.58830945 16 3.89647173 2.24282063 0.44115397 -0.77140624 -1.99094522 0.41455603 ^C cfa@abrams01:/home/cfa/dxu/che800-002/md1/T0=2.0>
The first line contains two numbers: The number of particles (512) and a “flag” which indicates whether velocities are included. (Note: this is actually not the standard XYZ format as originally defined; most programs that use the XYZ format don’t care about velocities.) Then a blank line appears – this is actually a line that is reserved for a descriptive title of the configuration. I am just too lazy to put one. Then, the actual configuration data begins: We have \(x\), \(y\), and \(z\) components of the position and velocity vectors for each particle, one per line. The number “16” at the beginning designates the “type” as an atomic number; for simplicity, I have decided to call all of my particles sulfur. (XYZ format is often used for atomically-specific configurations.) The functions xyz_out() and xyz_in() write and read this format, respectively, in mdlj.c. We will use these functions in other programs as well, typically those which analyze configuration data. Examples of such analysis codes are the subjects of the next two sections.
At this point, you can’t do much with all this data, except appreciate just how much data an MD code can produce. It is not unusual nowadays for researchers to use MD to produce hundreds of gigabytes of configuration data in order to write a single paper. It leads one to think that perhaps a lesson on handling large amounts of data is appropriate for a course on Molecular Simulation; however, I’ll forego that for now by trying to keep our sample exercises small.
One thing we can do with this data is make nice pictures using VMD. (I showed you this in class – a tutorial will appear soon.) Below are two renderings, one of the initial snapshot, and the other at time \(t\) = 450. Notice how the initially perfect crystalline lattice has been wiped out.
