Computer simulation methods have been applied to solve many-body problems in

physics, physical chemistry and biophysics; these methods are essentially experiments

in statistical mechanics. [30] The systems which statistical mechanics addresses are

composed of a large (i.e. statistical) number of particles; this branch of science allows

one to relate the collective properties of interesting atoms and molecules to the

macroscopic bulk (i.e. thermodynamic) properties of materials. In many cases, it is

impossible to study specific aspects of how atomic interactions manifest continuum

properties unless via computer simulation. Among simulation techniques, it is natural

to classify them according to the length scale at which the model is resolved. For

atomic scale models, there are two main families of simulation technique, Molecular

Dynamics (MD) method and Monte Carlo (MC) method [31, 32].

Molecular Dynamics simulation is a computational tool used to describe how

positions, velocities, and orientations of interacting particles change over time.

Particles may represent atoms, groups of atoms, or even larger entities; herein, focus

is on the use of MD as an atomic scale simulation technique. A set of models (i.e.

mathematic functions) that describe atomic scale interaction are foundations of any

MD simulation. Dynamics are governed by the system‘s Hamiltonian, and

Hamilton‘s equations of motion:

where H is the Hamiltonian, T is the kinetic energy which is a function of momenta, T=p*p/(2m), and U is the potential energy which is a function of position in the atomic

ensemble, U=U(q). The simulation model for atomic interactions is defined by an

interatomic potential energy function, or potential for brevity. [33] The potential can

be represented as 𝑈 = 𝑈(𝑟⃗ 𝑁), where 𝑟⃗ 𝑁 is the set of location vectors of an ensemble

of N atoms. For simplicity, in many simulations only pair interactions among atoms

are taken into account, so the potential energy of a system of N atoms is:

where 𝑢(𝑟𝑖𝑗) is a pair potential function having a known form, and 𝑟𝑖𝑗 is the

absolute value of distance between atom i and j. Since the interatomic forces are

derived from potential, they are all conservative, and can be related to the potential as

the gradient of U, or

In Molecular Dynamics, positions of atoms are attained by seeking the solution of

Newton‘s equation of motion:

where 𝐹𝑖

⃗⃗ is the force induced by the interactions between atom i and the other N-1

atoms, and the double dots mean second order time derivative while m is the mass of

each atom. Substitution of (2.3) into (2.4) relates the force to the potential function,

which is a known input to a given MD simulation. Integration once yields the

momenta of atoms, and the atomic position can be found by integrating an additional

time.

Repeating the integration procedure many times will reveal the trajectories of

atoms, based on which time average of some desired thermodynamic value,

, can
be computed

When a state of equilibrium is reached, the average will be independent from the

initial time 𝑡0. Because there are many times of integration involved to achieve

equilibrium, it is of great importance that an appropriate integration scheme and a

time steps to the algorithm are selected. Poor integration schemes or too large a time

step will introduce systematic errors in the thermodynamic ensemble and therefore in

the computed value of

.
MD simulation can be used to realize different sorts of thermodynamic ensembles,

which dictate the state thermodynamic quantities that will remain constant throughout

the simulation. Once the ensemble is established, time averages may be obtained, say,

for the thermodynamic state variables that are not held constant. For example, a

commonly employed MD ensemble is the microcanonical, [30] where

**number of**
**atoms N, volume V, and total energy E (sum of kinetic and potential energy) are held**
**constant (i.e. NVE ensemble)**; in such an ensemble, one may seek, for instance, time

averages of the system pressure P or temperature T. Other often employed ensembles

include the canonical, i.e. constant-NVT, ensemble, and the isothermal/isobaric, i.e. constant-NPT, ensemble. In order to implement a given ensemble in a MD simulation,

certain thermodynamic quantities must be held constant (either absolutely or in a time

averaged sense). Various numerical algorithms have been formulated to fix the

corresponding properties to a desired value. More theoretically detailed information

on ensembles can be found in [32] [33] and texts of statistical mechanics. The key

point is that properly implemental MD algorithms maintain thermodynamic

ensembles, form which methods of statistical mechanics reveal desired

thermodynamic quantities (i.e. material properties).

In short, the mechanism of the MD simulation process can be summarized as

following

1. Select the potential energy function that will be used to describe the

interatomic interactions

2. Set the initial position and velocity of each atom in a system and the

dimensions of the simulation space

3. Calculate all the forces on total N atoms via the governing interaction model

4. Find accelerations of each atom

5. Obtain each atom‘s velocity and position after a given timestep through

numerical integration by a chosen algorithm

6. Update all atomic positions

7. Repeat steps 3 - 6

8. Post process, or analyze, generated thermodynamic and atomic trajectory data

MD simulation was first introduced by B. Alder and T. Wainright in 1957 to study

the phase transition of hard spheres system, where particle propagation is determined

by each successive binary collision, rather than an integration of the equations of

motion. In 1964 a model system of liquid Argon was simulated by A. Rahman. In this

simulation not only the collision events between two particles but also interactions

modeled by Lennard-Jones potential were taken into account. This was the first time

that MD simulation was applied to atoms interacting through a continuous potential.

During 1970s, MD method continued to mature. Rahman and Stillinger (1971)

performed a simulation of 216 water molecules using an effective pair potential and

Ewald summation to calculate long-ranged forces triggered by electrostatic (i.e.

Coulomb) interactions. Several years later, the first protein simulation was carried out

by McCammon for bovine pancreatic tyrpsin inhibitor. [33] After that, MD

simulations began to be operated under different conditions and ensembles and

became more sophisticated. For instance, using MD to study non equilibrium

processes was first advanced in 1983 by Gillan and Dixon to calculate the

autocorrelation function of the microscopic heat current and hence the thermal

conductivity of a liquid. The Nose-Hoover thermostat algorithm, which is used to

obtain isothermal ensembles, was developed in 1984. The Car-Parrinello method to

carry out ab initio quantum mechanical MD (in which quantum mechanical effects for

the electronic degrees of freedom are taken into account) was developed in 1985. The

progress of available computer facilities is a major factor driving the development of

MD simulations. [31,32] Today MD methods are applied to a huge class of subjects,

e.g. surface properties, molecular clusters, biomolecules, and properties of liquids and

solids, and many different simulation codes have been developed.

Another well-developed simulation methodology is Monte Carlo (MC) method.

MC simulation is a computer simulation method that explores the phase space of a

system in a time-independent and stochastic manner. It uses a sequence of random

numbers to render and sample the phase space, from which diverse physical

properties of the system can be extracted. The main advantage of Monte Carlo

simulation over many analytical techniques is that it is a simple and straightforward

method that can provide an approximate solution to the problems encountered in

complex systems. Specifically, MC is able to reveal long time equilibrium properties.

The main advantage of MD method, with respect to MC method, is that MD

probes equilibrium configurations as well as dynamical properties of the system; this

makes MD well suited to compute time correlation functions and associated transport

coefficients. MD is therefore more efficient to evaluate some thermodynamic

properties like heat capacity, compressibility, and interfacial transport properties.

There are also computational advantages of MD method, particularly as compared to

ab initio methods. For example, the dynamic behavior of the atomic system is

described empirically without having to solve Schrodinger‘s equation at each time

step. Furthermore, it has been extensively discussed that many MD models are

inherently able to be executed in parallel. On the other hand, parallel implementations

of MC or ab initio based MD are highly non-trivial to realize.

Computing every force on each atom is a demanding part of a MD simulation. In

order to limit the number of force calculations, the smallest possible system – in terms

of number of atoms – should be simulated. The minimum system size necessary to

compute a given desired property is not always obvious and a common practice is to

analyze how increasing system size influences obtained results. A number of

algorithmic approaches exist that help minimize system size such as applying periodic

boundary conditions; these will be discussed in greater detail in the following section.

In addition to how large a system must be, one must determine how long a simulation

must be run to obtain reliable answers (i.e. with sufficient accuracy and precision).

Related to this is that the MD simulation time step should be small compared to the

period associated with the highest frequency atomic event. A typical time step for MD

simulations of metals is of order 1 fs (i.e. 1 E-15 s); this is roughly 1/1000 of the

period for an atomic oscillation. The very small time step is the main drawback of

MD method compared to other commonly used numerical simulations. Processes with

time scale of order microseconds and larger are outside the reasonable realm for even

highly parallel, efficient MD codes. Conversely, MC techniques cannot probe

dynamics but they can probe long time equilibrium properties.

Figure 2.1 Schematic comparison of time- and length-scales, accessible to different

types of simulation techniques (quantum simulations (QM), molecular dynamics

(MD), Brownian Dynamics (BD) and hydrodynamics/fluid dynamics (HD))

Accessible length- and time- scales of microscopic simulation should be

considered. Figure 2.1 shows a schematic representation for different types of

simulations in a length-time-diagram. It is clear that the more detailed a simulation

technique addresses, the smaller is the accessibility of times and length scales. Hence,

quantum simulations (QS), where fast motion of electrons are considered, are located

in the lower left corner of the diagram and typical length time scales are of order of Å

(Å) and ps. MD method puts fixed partial charges on interaction sites or adds a

dynamics model for polarization effects to simulate electronic distribution. [31]

Hence, the time scale of the system is determined by the time of collision between

atoms and atomic vibration frequencies. Consequently, the accessible time scale is of

order ns and corresponding accessible length scale becomes of order 10-10000 Å.

Note Figure 2.1 was originally presented in 2002; accessible length scales for MD

have increased dramatically. Brownian Dynamics (BD) can be applied to deal with

particles in a solvent medium, where one is not interested in a detailed description of

the solvent. In BD the effect of the solvent is hidden in average quantities. If one is

not interested in a molecular picture but in macroscopic quantities, the concepts of

hydrodynamics (HD) may be applied, where the system properties are hidden in

effective numbers, like density and sound velocity. These example methods, or levels

of theory, are geared towards liquids but analogous examples exist for solids.

Ref:

http://preserve.lehigh.edu/cgi/viewcontent.cgi?article=2257&context=etd
*Molecular Dynamics Simulation of Gold Nanoparticles and Surface Stress Effect*
*Siming Zhang*
*Lehigh University*
*2011*