Our simulation system is typically a parallelepiped supercell in which particles are placed either in a very regular manner, as in modeling a crystal lattice, or in some random manner, as in modeling a gas or liquid. For the simulation of perfect crystals, the number of particles in the simulation cell can be quite small, and only certain discrete values, such as 256, 500, and 864, should be specified. These numbers pertain to a face — centered-cubic crystal that has four atoms in each unit cell. If our simulation cell has l unit cells along each side, then the number of particles in the cube will be 4l3. The above numbers then correspond to cubes with 4, 5, and 6 cells along each side, respectively.

Once we have chosen the number of particles we want to simulate, the next step is to choose the system density we want to study. Choosing the density is equivalent to choosing the system volume since density p = N/O, where N is the number of particles and O is the supercell volume. An advantage of the Lennard-Jones potential is that one can work in dimensionless reduced units. The reduced density pa3 has typical values of about 0.9—1.2 for solids and 0.6-0.85 for liquids. For reduced temperature kB T/e, the values are 0.4-0.8 for solids and 0.8-1.3 for liquids. Notice that assigning particle velocities according to the Maxwellian velocity distribution probability = (m/2nkB T)3/2exp[-m( v2 + vj + ^)/2кв T]dvx dvy dvz is tantamount to setting the system temperature T.

For simulation of bulk properties (system with no free surfaces), it is conventional to use the periodic boundary condition (PBC). This means that the cubical simulation cell is surrounded by 26 identical image cells. For every particle in the simulation cell, there is a corresponding image particle in each image cell. The 26 image particles move in exactly the same manner as the actual particle, so if the actual particle should happen to move out of the simulation cell, the image particle in the image cell opposite to the exit side will move in and become the actual particle in the simulation cell. The net effect is that par­ticles cannot be lost or created. It follows then that the particle number is conserved, and if the simulation cell volume is not allowed to change, the system density remains constant.

Since in the pair potential approximation, the particles interact two at a time, a procedure is needed to decide which pair to consider among the pairs between actual particles and between actual and image particles. The minimum image convention is a procedure in which one takes the nearest neighbor to an actual particle as the interaction partner, regardless of whether this neighbor is an actual parti­cle or an image particle. Another approximation that is useful in keeping the computations to a manage­able level is the introduction of a force cutoff distance beyond which particle pairs simply do not see each other (indicated as rc in Figure 4). In order to avoid a particle interacting with its own image, it is necessary to set the cutoff distance to be less than half of the simulation cell dimension.

Another book-keeping device often used in MD simulation is a neighbor list to keep track of who are the nearest, second nearest, … neighbors of each particle. This is to save time from checking every particle in the system every time a force calculation is made. The list can be used for several time steps before updating. In low-temperature solids where the particles do not move very much, it is possible to do an entire simulation without, or with only a few, updating, whereas in simu­lation ofliquids, updation every 5 or 10 steps is common.

If one uses a naive approach in updating the neighbor list (an indiscriminate double loop over all particles), then it will get expensive for more than a few thousand particles because it involves N x N operations for an N-particle system. For short-range interactions, where the interatomic potential can be safely taken to be zero outside of a cutoff rc, acceler­ated approaches exist that can reduce the number of operations from order-N2 to order-N. For example, in the so-called ‘cell lists’ approach,18 one partitions the supercell into many smaller cells, and each cell maintains a registry of the atoms inside (order-N operation). The cell dimension is chosen to be greater than rc, so an atom cannot possibly interact with more than one neighbor atom. This will reduce the number of operations in updating the neighbor list to order-N.

With the so-called Parrinello-Rahman method,1 the supercell size and shape can change dynamically during a MD simulation to equilibrate the internal stress with the externally applied constant stress. In these simulations, the supercell is generally non­orthogonal, and it becomes much easier to use the so-called scaled coordinates sy to represent particle positions. The scaled coordinates Sj are related to the real coordinates ry through the relation, ry = H • sy, when both rj and sj are written as column vectors. H is a 3 x 3 matrix whose columns are the three repeat vectors of the simulation cell. Regardless of the shape of the simulation cell, the scaled coordi­nates of atoms can always be mapped into a unit cube, [0, 1) x [0, 1) x [0, 1). The shape change of the simulation cell with time can be accounted for by including the matrix H into the equation of motion. A ‘cell lists’ algorithm can still be worked out for a dynamically changing H, which minimizes the number of updates.13

For modeling ionic crystals, the long-range elec­trostatic interactions must be treated differently from short-ranged interactions (covalent, metallic, van der Waals, etc.). This is because a brute-force evaluation of the electrostatic interaction energies involves computation between all ionic pairs, which is of the order N2, and becomes very time­consuming for large N. The so-called Ewald sum — mation20,21 decomposes the electrostatic interaction into a short-ranged component, plus a long-ranged component, which, however, can be efficiently summed in the reciprocal space. It reduces the computational time to order N3/2. The particle mesh Ewald22-24 method further reduces the computational time to order N log N.

with t taken to be as long as possible. In terms of discrete time steps, eqn [10] becomes

1 1

(A )=yE A(tk) [11]

L k=1

where L is the number of time steps in the trajectory. The second is a time-dependent quantity of the form

1L0

(A(0)B(t))= A(tk)B(tk + t) [12]

L k=1

where B is in general another dynamical variable, and L’ is the number of time origins. Equation [12] is called a correlation function oftwo-dynamical vari­ables; since it is manifestly time dependent, it is able to represent dynamical information of the system.

We give examples of both types of averages by considering the properties commonly calculated in MD simulation.

N

U = V(ry) potential energy [13]

i<j

N

m Vi • Vi

i=1

@V (rj’) h @rj ij,

1N

g (r )=rwN(|=: j=:d(r

radial distribution function

1N

MSD(t) =Njri(t) ^ri (0)l2

N i = 1

mean squared displacement

In eqn [19], va is the average volume of one atom, va is the a-component of vector Vi, and rja is the a-component of vector ri — ry. The interest in writing the stress tensor in the present form is to suggest that the macroscopic tensor can be decomposed into individual atomic contributions, and thus is known as the atomic level stress25 at atom i. Although this interpretation is quite appealing, one should be aware that such a decomposition makes sense only in a nearly homogeneous system where every atom ‘owns’ almost the same volume as every other atom. In an inhomogeneous system, such as in the vicinity of a surface, it is not appropriate to consider such decomposition. Both eqns [15] and [19] are written for pair potential models only. A slightly different expression is required for potentials that contain many-body terms.26