In the simplest physical terms, MD may be charac­terized as a method of ‘particle tracking.’ Operation­ally, it is a method for generating the trajectories of a system of N particles by direct numerical integration

of Newton’s equations of motion, with appropriate specification of an interatomic potential and suitable initial and boundary conditions. MD is an atomistic modeling and simulation method when the particles in question are the atoms that constitute the material of interest. The underlying assumption is that one can treat the ions and electrons as a single, classical entity. When this is no longer a reasonable approxi­mation, one needs to consider both ion and electron motions. One can then distinguish two versions of MD, classical and ab initio, the former for treating atoms as classical entities (position and momentum) and the latter for treating separately the electronic and ionic degrees of freedom, where a wave func­tion description is used for the electrons. In this chapter, we are concerned only with classical MD. The use of ab initio methods in nuclear materials research is addressed elsewhere (Chapter 1.08, Ab Initio Electronic Structure Calculations for Nuclear Materials). Figure 2 illustrates the MD sim­ulation system as a collection of N particles contained in a volume O. At any instant of time t, the particle coordinates are labeled as a 3N-dimensional vector, r3N(t) = {ri(t), r2(t), …, rN(t)}, where r, repre­sents the three coordinates of atom i. The simula­tion proceeds with the system in a prescribed initial configuration, r3N(t0), and velocity, r3N(t0), at time t = tu As the simulation proceeds, the particles evolve through a sequence of time steps, r3N (t0) ! r3N (t1) ! r3N(t2) ! ••• ! r3N(tL), where tk = t0 + kAt,

k = 1,2,…, L, and At is the time step ofMD simulation. The simulation runs for L number of steps and covers a time interval of LAt. Typical values of L can range from 104 to 108 and At ~ 10~15s. Thus, nominal MD simulations follow the system evolution over time inter­vals not more than ~ 1-10 ns.

The simulation system has a certain energy E, the sum of the kinetic and potential energies of the particles, E = K + U where K is the sum of individual kinetic energies

1 ^

K = 2m v • v [1]

y=1

and U = U(r3N) is a prescribed interatomic interac­tion potential. Here, for simplicity, we assume that all particles have the same mass m. In principle, the potential U is a function of all the particle coordi­nates in the system if we allow each particle to interact with all the others without restriction. Thus, the dependence of U on the particle coordinates can be as complicated as the system under study demands. However, for the present discussion we introduce an approximation, the assumption of a two-body or pair-wise additive interaction, which is sufficient to illustrate the essence of MD simulation.

To find the atomic trajectories in the classical version of MD, one solves the equations governing the particle coordinates, Newton’s equations of motion in mechanics. For our N-particle system with potential energy U, the equations are

d2 r

m =-Vr;U (r3N ); j = 1; …; N [2]

where m is the particle mass. Equation [2] may look deceptively simple; actually, it is as complicated as the famous N-body problem that one generally cannot solve exactly when N is >2. As a system of coupled second-order, nonlinear ordinary differential equations, eqn [2] can be solved numerically, which is what is carried out in MD simulation.

Equation [2] describes how the system (particle coordinates) evolves over a time period from a given initial state. Suppose we divide the time period ofinter — est into many small segments, each being a time step of size At. Given the system conditions at some initial time t0, r3N(t0), and r3N(t0), integration means we advance the system successively by increments of At,

r3N (to) ! r3N (tl) ! r3N (t2) ! • • • ! r3N (tL) [3]

where L is the number of time steps making up the interval of integration.

How do we numerically integrate eqn [3] for a given U? A simple way is to write a Taylor series expansion,

r, (to + At) = r, (to) + v, (to)At

2 [4]

+ l/2ay (to)(At) + •••

and a similar expansion for ry (t0 — At). Adding the two expansions gives

r, (to + At)= — rj (to — At) + 2ry (to)

+ ay (to)(At) + •••

Notice that the left-hand side of eqn [5] is what we want, namely, the position of particle j at the next time step to + At. We already know the positions at to and the time step before, so to use eqn [5] we need the acceleration of particle y at time to. For this we substitute Fy(r3N(to))/m in place of acceleration ay(to), where Fy is just the right-hand side of eqn [2]. Thus, the integration of Newton’s equations of motion is accomplished in successive time incre­ments by applying eqn [5]. In this sense, MD can be regarded as a method of particle tracking where one follows the system evolution in discrete time steps. Although there are more elaborate, and therefore more accurate, integration procedures, it is impor­tant to note that MD results are as rigorous as classical mechanics based on the prescribed interatomic poten­tial. The particular procedure just described is called the Verlet (leapfrog)lci method. It is a symplectic integrator that respects the symplectic symmetry of the Hamiltonian dynamics; that is, in the absence of floating-point round-off errors, the discrete mapping rigorously preserves the phase space volume.1 , Symplectic integrators have the advantage of long­term stability and usually allow the use of larger time steps than nonsymplectic integrators. However, this advantage may disappear when the dynamics is not strictly Hamiltonian, such as when some thermostat — ing procedure is applied. A popular time integrator used in many early MD codes is the Gear predictor- corrector method13 (nonsymplectic) of order 5. Higher accuracy of integration allows one to take a larger value of At so as to cover a longer time interval for the same number of time steps. On the other hand, the trade-off is that one needs more computer mem­ory relative to the simpler method.

A typical flowchart for an MD code11 would look something like Figure 3. Among these steps, the part that is the most computationally demanding is the force calculation. The efficiency of an MD simulation therefore depends on performing the force calculation as simply as possible without compromising the phys­ical description (simulation fidelity). Since the force is calculated by taking the gradient of the potential U, the specification of U essentially determines the com­promise between physical fidelity and computational efficiency.

Figure 3 Flow chart of MD simulation.

1.09.2 The Interatomic Potential

This is a large and open-ended topic with an exten­sive literature.14 It is clear from eqn [2] that the interaction potential is the most critical quantity in MD modeling and simulation; it essentially controls the numerical and algorithmic simplicity (or complex­ity) of MD simulation and, therefore, the physical fidelity of the simulation results. Since Chapter 1.10, Interatomic Potential Development is devoted to interatomic potential development, we limit our dis­cussion only to simple classical approximations to U(rb r2, . . . , rN).

Practically, all atomistic simulations are based on the Born-Oppenheimer adiabatic approximation, which separates the electronic and nuclear motions.15 Since electrons move much more quickly because of their smaller mass, during their motion one can treat the nuclei as fixed in instantaneous positions, or equivalently the electron wave functions follow the nuclear motion adiabatically. As a result, the electrons are treated as always in their ground state as the nuclei move.

For the nuclear motions, we consider an expansion of U in terms of one-body, two-body, … N-body interactions:

NN

U (r3N )= V1(rj )+ V2(ri; rj )

}=1 i<J

[6]

N

+ V3(ri; rj ; rkH———-

i<j <k

The first term, the sum of one-body interactions, is usually absent unless an external field is present to couple with each atom individually. The second sum is the contribution ofpure two-body interactions (pairwise additive). For some problems, this term alone is sufficient to be an approximation to U. The third sum represents pure three-body interactions, and so on.