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
