A molecular dynamics primer

Furio Ercolessi
International School for Advanced Studies (SISSA-ISAS)
Trieste, Italy

Email: furio@sissa.it
WWW: http://www.sissa.it/furio/

Spring College in Computational Physics, ICTP, Trieste, June 1997.

Contents
Preface 1
1 Introduction 2
1.1 The role of computer experiments . . . . . . . . . . . . . . . . . . 2
1.1.1 But is it theory or experiment? . . . . . . . . . . . . . . . 3
1.2 What is molecular dynamics? . . . . . . . . . . . . . . . . . . . . 4
1.3 Historical notes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.4 Today's role of molecular dynamics . . . . . . . . . . . . . . . . . 6
1.5 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.5.1 Use of classical forces . . . . . . . . . . . . . . . . . . . . 7
1.5.2 Realism of forces . . . . . . . . . . . . . . . . . . . . . . . 8
1.5.3 Time and size limitations . . . . . . . . . . . . . . . . . . 8
2 The basic machinery 10
2.1 Modeling the physical system . . . . . . . . . . . . . . . . . . . . 10
2.1.1 The Lennard-Jones potential . . . . . . . . . . . . . . . . 10
2.1.2 Potential truncation and long-range corrections . . . . . . 11
2.2 Periodic boundary conditions . . . . . . . . . . . . . . . . . . . . 12
2.2.1 The minimum image criterion . . . . . . . . . . . . . . . . 13
2.2.2 Geometries with surfaces . . . . . . . . . . . . . . . . . . 14
2.3 Time integration algorithm . . . . . . . . . . . . . . . . . . . . . 14
2.3.1 The Verlet algorithm . . . . . . . . . . . . . . . . . . . . . 15
2.3.2 Predictor-corrector algorithm . . . . . . . . . . . . . . . . 16
3 Running, measuring, analyzing 18
3.1 Starting a simulation . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.1.1 Starting from scratch . . . . . . . . . . . . . . . . . . . . 18
3.1.2 Continuing a simulation . . . . . . . . . . . . . . . . . . . 19
3.2 Controlling the system . . . . . . . . . . . . . . . . . . . . . . . . 19
3.3 Equilibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.4 Looking at the atoms . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.5 Simple statistical quantities to measure . . . . . . . . . . . . . . 21
3.5.1 Potential energy . . . . . . . . . . . . . . . . . . . . . . . 22
3.5.2 Kinetic energy . . . . . . . . . . . . . . . . . . . . . . . . 22
3.5.3 Total energy . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.5.4 Temperature . . . . . . . . . . . . . . . . . . . . . . . . . 23

i
 3.5.5 The caloric curve . . . . . . . . . . . . . . . . . . . . . . . 23
3.5.6 Mean square displacement . . . . . . . . . . . . . . . . . . 24
3.5.7 Pressure . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.6 Measuring the melting temperature . . . . . . . . . . . . . . . . . 25
3.7 Real space correlations . . . . . . . . . . . . . . . . . . . . . . . . 26
3.8 Reciprocal space correlations . . . . . . . . . . . . . . . . . . . . 27
3.9 Dynamical analysis . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.10 Annealing and quenching: MD as an optimization tool . . . . . . 29
3.11 Other statistical ensembles . . . . . . . . . . . . . . . . . . . . . 30
4 Interatomic potentials 32
4.1 The Born-Oppenheimer approximation . . . . . . . . . . . . . . . 32
4.2 The design of potentials . . . . . . . . . . . . . . . . . . . . . . . 33
4.3 The problems with two-body potentials . . . . . . . . . . . . . . 34
4.4 Many-body potentials for metals . . . . . . . . . . . . . . . . . . 35
4.5 Many-body potentials for semiconductors . . . . . . . . . . . . . 37
4.5.1 The Stillinger-Weber potential . . . . . . . . . . . . . . . 37
4.5.2 The Terso potential . . . . . . . . . . . . . . . . . . . . . 38
4.6 Long-range forces . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
4.7 But do we really need potentials? . . . . . . . . . . . . . . . . . . 39
4.8 Fitting to ab initio data by \force matching" . . . . . . . . . . . 40
A Choosing hardware and software 42
A.1 What kind of computer? . . . . . . . . . . . . . . . . . . . . . . . 42
A.1.1 Selecting a computer for molecular dynamics . . . . . . . 42
A.1.2 Storage space . . . . . . . . . . . . . . . . . . . . . . . . . 44
A.2 What language? . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
A.2.1 High level languages . . . . . . . . . . . . . . . . . . . . . 44
A.2.2 Fortran or C? . . . . . . . . . . . . . . . . . . . . . . . . . 45
Bibliography 47

ii
Preface
These notes constitute the skeleton for a short hands-on course on molecular
dynamics. It is assumed that the reader has basic notions of classical mechanics,
statistical mechanics and thermodynamics, but knows very little|if anything
at all|on methods for atomistic computer simulation. The basics of molecular
dynamics are explained, together with a few tricks of the trade.
No attempt is made to address \advanced topics" in molecular dynamics,
such as calculation of free energies, non equilibrium simulations, etc., and many
important subjects are treated only super cially. This document is meant to be
just a starting point. Many excellent books devoted to computer simulations
in condensed matter physics are available [1, 2, 3, 4, 5, 6, 7], and the reader is
advised to make use of them to dig into more speci c subjects.
These note are accompanied by example programs written in
Fortran 90 and available on the World Wide Web at the URL
http://www.sissa.it/furio/md/f90. While still relatively uncommon
among scientists, Fortran 90 is bound to become one of the most important
languages for scienti c programming in the years to come. Fortran 90
incorporates Fortran 77|allowing it to inherit smoothly the legacy of the
immense body of existing codes|but also features elements now recognized to
be essential in a language for large projects in scienti c computing.
I hope that the examples contained here will be of help to young compu-
tational physicists to get familiar with Fortran 90, and to people grown with
Fortran 77 (or 66!) to migrate towards this new, more modern environment.
Any comment or suggestion to improve this document, which will also
be kept online at the URL http://www.sissa.it/furio/md/, will be wel-
come and appreciated. The best way to contact the author is by email at
furio@sissa.it.

1
Chapter 1
Introduction
1.1 The role of computer experiments
Computer experiments play a very important role in science today. In the past,
physical sciences were characterized by an interplay between experiment and
theory. In experiment, a system is subjected to measurements, and results,
expressed in numeric form, are obtained. In theory, a model of the system is
constructed, usually in the form of a set of mathematical equations. The model
is then validated by its ability to describe the system behavior in a few selected
cases, simple enough to allow a solution to be computed from the equations.
In many cases, this implies a considerable amount of simpli cation in order to
eliminate all the complexities invariably associated with real world problems,
and make the problem solvable.
In the past, theoretical models could be easily tested only in a few simple
\special circumstances". So, for instance, in condensed matter physics a model
for intermolecular forces in a speci c material could be veri ed in a diatomic
molecule, or in a perfect, in nite crystal. Even then, approximations were often
required to carry out the calculation. Unfortunately, many physical problems of
extreme interest (both academic and practical) fall outside the realm of these
\special circumstances". Among them, one could mention the physics and
chemistry of defects, surfaces, clusters of atoms, organic molecules, involving
a large amount of degrees of freedom; an accurate treatment of temperature
e ects, including anharmonicities and phase transitions; disordered systems in
general, where symmetry is of no help to simplify the treatment; and so on.
The advent of high speed computers|which started to be used in the 50s|
altered the picture by inserting a new element right in between experiment and
theory: the computer experiment. In a computer experiment, a model is still
provided by theorists, but the calculations are carried out by the machine by
following a \recipe" (the algorithm, implemented in a suitable programming
language). In this way, complexity can be introduced (still with caution!) and
more realistic systems can be investigated, opening a road towards a better
understanding of real experiments.
Needless to say, the development of computer experiments altered substan-
tially the traditional relationship between theory and experiment. On one side,

2
computer simulations increased the demand for accuracy of the models. For
instance, a molecular dynamics simulation allows to evaluate the melting tem-
perature of a material, modeled by means of a certain interaction law. This
is a dicult test for the theoretical model to pass|and a test which has not
been available in the past. Therefore, simulation \brings to life" the models,
disclosing critical areas and providing suggestions to improve them.
On the other side, simulation can often come very close to experimental
conditions, to the extent that computer results can sometimes be compared
directly with experimental results. When this happens, simulation becomes an
extremely powerful tool not only to understand and interpret the experiments
at the microscopic level, but also to study regions which are not accessible
experimentally, or which would imply very expensive experiments, such as under
extremely high pressure.
Last but not least, computer simulations allow thought experiments|things
which are just impossible to do in reality, but whose outcome greatly increases
our understanding of phenomena|to be realized. Fantasy and creativity are
important qualities for the computer simulator!
1.1.1 But is it theory or experiment?
Simulation is seen sometimes as theory, sometimes as experiment. On one
side, we are still dealing with models, not with the \real thing": this suggests
to classify simulation as belonging to theoretical methods without hesitation.
On the other side, the procedure of verifying a model by computer simulation
resembles an experiment quite closely: we perform a run, and then analyze the
results in pretty much the same way as experimental physicists do. So how
should we classify simulation?
There is no sharp answer to this question: both sides represent legitimate
point of views, and this is what makes computational science a branch on its
own. There is however another important consideration.
Theory is traditionally based on the reductionistic approach: we deal with
complexity by reducing a system to simpler subsystems, continuing until the
subsystems are simple enough to be represented with solvable models. When
we look at simulation as simply a practical tool to \verify and test" models
in situations which are too complex to handle analytically (for example, when
computing the phase diagram of a substance modeled by a certain force law),
we are implicitly assuming that the model represents the \theory level" where
the interest is focused.
But it is important to realize that simulation may play a more important
and interesting role. We can consider it not as an aid to reductionism but|
to some extent|as an alternative to it. Simulation increases the threshold of
complexity which separates \solvable" and \unsolvable" models. We can take
advantage of this threshold shift and move up one level in our description of
physical systems. Thanks to the presence of simulation, we do not need to work
with models as simple as those used in the past. This gives us an additional
degree of freedom to explore and opens entirely new possibilities.
One example of this point of view is represented by interatomic potentials.

3
In the past, interactions were obtained by two-body potentials with simple
analytical form, such as Morse or Lennard-Jones. Today, the most accurate
potentials contain many-body terms and are determined numerically by repro-
ducing as closely as possible forces predicted by rst-principle methods (this is
discussed in 4.8). We have thus moved up one level in the degree of reduction-
ism contained in the potential, now limited only to the selection of its analytical
form. The advantage is of course a much better realism, which in turn allows
investigation of physics problem which require a level of accuracy in the model
that was not achieved before. These new potentials could not exist without
simulation: simulation is not only a link between experiment and theory, it is
also a poweful tool to propel progress in new directions.
The reader interested in these \philosophical" aspects of computational sci-
ence can nd a very enjoyable discussion in chapter 1 of ref. [6].

1.2 What is molecular dynamics?

We call molecular dynamics (MD) a computer simulation technique where the
time evolution of a set of interacting atoms is followed by integrating their
equations of motion.
In molecular dynamics we follow the laws of classical mechanics, and most
notably Newton's law:
Fi = miai (1.1)
for each atom i in a system constituted by N atoms. Here, mi is the atom
mass, ai = d2 ri =dt2 its acceleration, and Fi the force acting upon it, due
to the interactions with other atoms. Therefore, in contrast with the Monte
Carlo method, molecular dynamics is a deterministic technique: given an initial
set of positions and velocities, the subsequent time evolution is in principle1
completely determined. In more pictorial terms, atoms will \move" into the
computer, bumping into each other, wandering around (if the system is uid),
oscillating in waves in concert with their neighbors, perhaps evaporating away
from the system if there is a free surface, and so on, in a way pretty similar to
what atoms in a real substance would do.
The computer calculates a trajectory in a 6N -dimensional phase space (3N
positions and 3N momenta). However, such trajectory is usually not particu-
larly relevant by itself. Molecular dynamics is a statistical mechanics method.
Like Monte Carlo, it is a way to obtain a set of con gurations distributed ac-
cording to some statistical distribution function, or statistical ensemble. An
example is the microcanonical ensemble, corresponding to a probability density
in phase space where the total energy is a constant E :
(H (;) ; E ):
Here, H (;) is the Hamiltonian, and ; represents the set of positions and mo-
menta.  is the Dirac function, selecting out only those states which have a
1
In practice, the niteness of the integration time step and arithmetic rounding errors will
eventually cause the computed trajectory to deviate from the true trajectory.

4
speci c energy E . Another example is the canonical ensemble, where the tem-
perature T is constant and the probability density is the Boltzmann function
exp(;H (;)=kB T ):
According to statistical physics, physical quantities are represented by aver-
ages over con gurations distributed according to a certain statistical ensemble.
A trajectory obtained by molecular dynamics provides such a set of con gura-
tions. Therefore, a measurements of a physical quantity by simulation is simply
obtained as an arithmetic average of the various instantaneous values assumed
by that quantity during the MD run.
Statistical physics is the link between the microscopic behavior and ther-
modynamics. In the limit of very long simulation times, one could expect the
phase space to be fully sampled, and in that limit this averaging process would
yield the thermodynamic properties. In practice, the runs are always of -
nite length, and one should exert caution to estimate when the sampling may
be good (\system at equilibrium") or not. In this way, MD simulations can
be used to measure thermodynamic properties and therefore evaluate, say, the
phase diagram of a speci c material.
Beyond this \traditional" use, MD is nowadays also used for other purposes,
such as studies of non-equilibrium processes, and as an ecient tool for opti-
mization of structures overcoming local energy minima (simulated annealing).

1.3 Historical notes

A full account of early developments of the molecular dynamics technique is
certainly beyond the scope of this document. I will simply mention below a few
key papers appeared in the 50s and in the 60s that can be regarded as milestones
in molecular dynamics, and whose reading is certainly recommended. I will not
mention developments in the Monte Carlo method, which also occurred during
the same period of time. The interested reader can nd a rst-hand account by
Wood in [9].
Reprints of all these papers, and of many other important works in the area
of computer simulation of liquids and solids published up to 1986, can be found
in ref. [8].
 The rst paper reporting a molecular dynamics simulation was written
by Alder and Wainwright [10] in 1957. The purpose of the paper was to
investigate the phase diagram of a hard sphere system, and in particular
the solid and liquid regions. In a hard sphere system, particles interact
via instantaneous collisions, and travel as free particles between collisions.
The calculations were performed on a UNIVAC and on an IBM 704.
 The article Dynamics of radiation damage by J. B. Gibson, A. N. Goland,
M. Milgram and G. H. Vineyard from Brookhaven National Laboratory,
appeared in 1960 [11], is probably the rst example of a molecular dy-
namics calculation with a continuous potential based on a nite di erence

5
 time integration method. The calculation for a 500-atoms system was
performed on an IBM 704, and took about a minute per time step. The
paper, dealing with creation of defects induced by radiation damage (a
theme appropriate to cold war days), is done exceedingly well, and is hard
to believe that it is almost 40 years old!
 Aneesur Rahman at Argonne National Laboratory has been a well known
pioneer of molecular dynamics. In his famous 1964 paper Correlations in
the motion of atoms in liquid argon [12], he studies a number of properties
of liquid Ar, using the Lennard-Jones potential on a system containing
864 atoms and a CDC 3600 computer. The legacy of Rahman's computer
codes is still carried by many molecular dynamics programs in operation
around the world, descendant of Rahman's.
 Loup Verlet calculated in 1967 [13] the phase diagram of argon using
the Lennard-Jones potential, and computed correlation functions to test
theories of the liquid state. The bookkeeping device which became known
as Verlet neighbor list was introduced in these papers. Moreover the
\Verlet time integration algorithm" (2.3.1) was used. Phase transitions
in the same system were investigated by Hansen and Verlet a couple of
years later [14].

1.4 Today's role of molecular dynamics

Reviewing even a tiny subset of the applications of molecular dynamics is far
beyond the purpose of these notes. I will only brie y mention a few areas of
current interest where MD has brought and/or could bring important contri-
butions. The list should not be considered as exhaustive, and certainly re ects
my background in condensed matter physics:
Liquids. As shown in 1.3, liquids are where everything started! Nevertheless,
they remain an important subject. Availability of new realistic interac-
tion models allows to study new systems, elemental and multicomponent.
Through non-equilibrium techniques, transport phenomena such as vis-
cosity and heat ow have been investigated [2].
Defects. Another subject whose investigation by MD started a long time ago
(see [8]), defects in crystals|crucial for their mechanical properties and
therefore of technological interest|remain a favoured topic. The focus
shifted perhaps from point defects (vacancies, interstitials) to linear (dislo-
cations) and planar (grain boundaries, stacking faults) defects [15]. Again,
improved realism thanks to better potentials constitutes a driving force.
Fracture. Under mechanical action, solids break into two or more pieces. The
fracture process can occur in di erent ways and with di erent speeds de-
pending of several parameters. The technological importance is obvious,
and simulation is providing insight.

6
Surfaces. Surface physics had a boom starting in the 80s, thanks to the avail-
ability of new wonderful experimental tools with microscopic resolution
(scanning tunneling microscopy, high resolution electron microscopy, sev-
eral scattering-based techniques). Simulation is still playing a big role in
understanding phenomena such as surface reconstructions, surface melt-
ing, faceting, surface di usion, roughening, etc, often requiring large sam-
ples and simulation times.
Friction. Even more recent are investigations of adhesion and friction between
two solids, propelled by the development of the atomic force microscope
(AFM). The body of \macroscopic" knowledge is being revised and ex-
panded on microscopic grounds.
Clusters. Clusters|conglomerates of a number of atoms ranging from a few
to several thousands|constitute a bridge between molecular systems and
solids, and exhibit challenging features. Frequently, an astonishingly large
number of di erent con gurations have very similar energies, making it
dicult to determine stable structures. Their melting properties can also
be signi cantly di erent from those of the solid, due to the nite size, the
presence of surfaces and the anisotropy. Metal clusters are extremely im-
portant from the technological point of view, due to their role as catalysts
in important chemical reactions (for instance, in catalytic exhaust pipes
of cars).
Biomolecules. MD allows to study the dynamics of large macromolecules,
including biological systems such as proteins, nucleic acids (DNA, RNA),
membranes [16]. Dynamical events may play a key role in controlling
processes which a ect functional properties of the biomolecule. Drug
design is commonly used in the pharmaceutical industry to test properties
of a molecule at the computer without the need to synthesize it (which is
far more expensive).
Electronic properties and dynamics. The development of the Car-Parri-
nello method (see 1.5.2), where the forces on atoms are obtained by solv-
ing the electronic structure problem instead of by an interatomic po-
tential, allows to study electronic properties of materials fully including
their dynamics (and, therefore, phase transitions and other temperature-
dependent phenomena). This important work gave rise to a very success-
ful research line during the last decade.

1.5 Limitations
Molecular dynamics is a very powerful technique but has|of course|limitations.
We quickly examine the most important of them.
1.5.1 Use of classical forces
One could immediately ask: how can we use Newton's law to move atoms, when
everybody knows that systems at the atomistic level obey quantum laws rather

7
than classical laws, and that Schrodinger's equation is the one to be followed?
A simple test of the validity of the classical approximation is based on the
de Broglie thermal wavelength [17], de ned as
s

 = Mk 2h 2 (1.2)
BT
where M is the atomic mass and T the temperature. The classical approxima-
tion is justi ed if   a, where a is the mean nearest neighbor separation. If
one considers for instance liquids at the triple point, =a is of the order of 0.1
for light elements such as Li and Ar, decreasing further for heavier elements.
The classical approximation is poor for very light systems such as H2 , He, Ne.
Moreover, quantum e ects become important in any system when T is suf-
ciently low. The drop in the speci c heat of crystals below the Debye temper-
ature [18], or the anomalous behavior of the thermal expansion coecient, are
well known examples of measurable quantum e ects in solids.
Molecular dynamics results should be interpreted with caution in these re-
gions.
1.5.2 Realism of forces
How realistic is a molecular dynamics simulation?
In molecular dynamics, atoms interact with each other. These interactions
originate forces which act upon atoms, and atoms move under the action of
these instantaneous forces. As the atoms move, their relative positions change
and forces change as well.
The essential ingredient containing the physics is therefore constituted by
the forces. A simulation is realistic|that is, it mimics the behavior of the real
system|only to the extent that interatomic forces are similar to those that real
atoms (or, more exactly, nuclei) would experience when arranged in the same
con guration.
As described in 2.1, forces are usually obtained as the gradient of a potential
energy function, depending on the positions of the particles. The realism of the
simulation therefore depends on the ability of the potential chosen to reproduce
the behavior of the material under the conditions at which the simulation is run.
The problem of selecting|or constructing|potentials is addressed in more
detail in chapter 4.
1.5.3 Time and size limitations
Typical MD simulations can be performed on systems containing thousands|
or, perhaps, millions|of atoms, and for simulation times ranging from a few
picoseconds to hundreds of nanoseconds. While these numbers are certainly
respectable, it may happen to run into conditions where time and/or size limi-
tations become important.
A simulation is \safe" from the point of view of its duration when the
simulation time is much longer than the relaxation time of the quantities we
are interested in. However, di erent properties have di erent relaxation times.

8
In particular, systems tend to become slow and sluggish in the proximity of
phase transitions, and it is not uncommon to nd cases where the relaxation
time of a physical property is orders of magnitude larger than times achievable
by simulation.
A limited system size can also constitute a problem. In this case one has
to compare the size of the MD cell with the correlation lengths of the spatial
correlation functions of interest. Again, correlation lengths may increase or
even diverge in proximity of phase transitions, and the results are no longer
reliable when they become comparable with the box length.
This problem can be partially alleviated by a method known as nite size
scaling. This consist of computing a physical property A using several box with
di erent sizes L, and then tting the results on a relation
A(L) = A + Lcn
using A ; c; n as tting parameters. A then corresponds to limL!1 A(L), and
should therefore be taken as the most reliable estimate for the \true" physical
quantity.

9
Chapter 2
The basic machinery
2.1 Modeling the physical system
The main ingredient of a simulation is a model for the physical system. For
a molecular dynamics simulation this amounts to choosing the potential: a
function V (r1 ; : : : ; rN ) of the positions of the nuclei, representing the potential
energy of the system when the atoms are arranged in that speci c con gura-
tion. This function is translationally and rotationally invariant, and is usually
constructed from the relative positions of the atoms with respect to each other,
rather than from the absolute positions.
Forces are then derived as the gradients of the potential with respect to
atomic displacements:
Fi = ;rri V (r1 ; : : : ; rN ) (2.1)
This form implies the presence of a conservation law of the total energy E =
K + V , where K is the instantaneous kinetic energy.
The simplest choice for V is to write it as a sum of pairwise interactions:
V (r1 ; : : : ; rN ) = (jri ; rj j)
XX
(2.2)
i j>i
The clause j > i in the second summation has the purpose of considering each
atom pair only once. In the past most potentials were constituted by pairwise
interactions, but this is no longer the case. It has been recognized that the two-
body approximation is very poor for many relevant systems, such as metals and
semiconductors. Various kinds of many-body potentials are now of common use
in condensed matter simulation, and will be brie y reviewed in x4.4 and x4.5.
The development of accurate potentials represents an important research
line. A quick overview of the current situation will be presented in Chapter 4.
For now, let us mention only the most commonly used interaction model: the
Lennard-Jones pair potential.
2.1.1 The Lennard-Jones potential
The Lennard-Jones 12-6 potential is given by the expression
" #

LJ (r) = 4"  12 ;   6 (2.3)

r r

10
for the interaction potential between a pair of atoms. The total potential of a
system containing many atoms is then given by 2.2.
This potential has an attractive tail at large r, it reaches a minimum around
1.122 , and it is strongly repulsive at shorter distance, passing through 0 at
r =  and increasing steeply as r is decreased further.
The term  1=r12 , dominating at short distance, models the repulsion be-
tween atoms when they are brought very close to each other. Its physical origin
is related to the Pauli principle: when the electronic clouds surrounding the
atoms starts to overlap, the energy of the system increases abruptly. The expo-
nent 12 was chosen exclusively on a practical basis: equation (2.3) is particularly
easy to compute. In fact, on physical grounds an exponential behavior would
be more appropriate.
The term  1=r6 , dominating at large distance, constitute the attractive
part. This is the term which gives cohesion to the system. A 1=r6 attraction
is originated by van der Waals dispersion forces, originated by dipole-dipole
interactions in turn due to uctuating dipoles. These are rather weak interac-
tions, which however dominate the bonding character of closed-shell systems,
that is, rare gases such as Ar or Kr. Therefore, these are the materials that a
LJ potential could mimic fairly well1. The parameters " and  are chosen to t
the physical properties of the material.
On the other hand, a LJ potential is not at all adequate to model situa-
tions with open shells, where strong localized bonds may form (as in covalent
systems), or where there is a delocalized \electron sea" where the ions sit (as
in metals). In these systems the two-body interactions scheme itself fails very
badly. Potentials for these systems will be discussed in chapter 4.
However, regardless of how well it is able to model actual materials, the
LJ 12-6 potential constitutes nowadays an extremely important model system.
There is a vast body of papers who investigated the behavior of atoms interact-
ing via LJ on a variety of di erent geometries (solids, liquids, surfaces, clusters,
two-dimensional systems, etc). One could say that LJ is the standard potential
to use for all the investigations where the focus is on fundamental issues, rather
than studying the properties of a speci c material. The simulation work done
on LJ systems helped us (and still does) to understand basic points in many
areas of condensed matter physics, and for this reason the importance of LJ
cannot be underestimated.
When using the LJ potentials in simulation, it is customary to work in a
system of units where  = 1 and " = 1. The example codes accompanying these
notes follow this convention.
2.1.2 Potential truncation and long-range corrections
The potential in eq. (2.3) has an in nite range. In practical applications, it
is customary to establish a cuto radius Rc and disregard the interactions be-
tween atoms separated by more than Rc . This results in simpler programs and
1
Many-body e ects in the interaction are however present also in these systems and, in all
cases, pair potentials more accurate than LJ have been developed for rare gases.

11
enormous savings of computer resources, because the number of atomic pairs
separated by a distance r grows as r2 and becomes quickly huge.
A simple truncation of the potential creates a new problem though: when-
ever a particle pair \crosses" the cuto distance, the energy makes a little jump.
A large number of these events is likely to spoil energy conservation in a simu-
lation. To avoid this problem, the potential is often shifted in order to vanish
at the cuto radius 2 :
(

V (r) = 0 LJ(r) ; LJ(Rc ) ifif rr >

 Rc
Rc (2.4)

Physical quantities are of course a ected by the potential truncation. The

e ects of truncating a full-ranged potential can be approximately estimated by
treating the system as a uniform (constant density) continuum beyond Rc. For
a bulk system (periodic boundary conditions along each direction, see x2.2), this
usually amounts to a constant additive correction. For example, the potential
tail (attractive) brings a small additional contribution to the cohesive energy,
and to the total pressure. Truncation e ects are not so easy to estimate for
geometries with free surfaces, due to the lower symmetry, and can be rather
large for quantities like surface energy.
Commonly used truncation radii for the Lennard-Jones potential are 2:5 
and 3:2 . It should be mentioned that these truncated Lennard-Jones models
are so popular that they acquired a value on their own as reference models for
generic two-body systems (as also discussed at the end of the previous section).
In many cases, there is no interest in evaluating truncation corrections because
the truncated model itself is the subject of the investigation.
Potentials for metals and semiconductors are usually designed from the start
with a cuto radius in mind, and go usually to zero at Rc together with their
rst two derivatives at least. Such potentials do not therefore contain true van
der Waals forces. Since such forces are much weaker than those associated with
metallic or covalent bonding, this is usually not a serious problem except for
geometries with two or more separate bodies.

2.2 Periodic boundary conditions

What should we do at the boundaries of our simulated system? One possibility
is doing nothing special: the system simply terminates, and atoms near the
boundary would have less neighbors than atoms inside. In other words, the
sample would be surrounded by surfaces.
Unless we really want to simulate a cluster of atoms, this situation is not
realistic. No matter how large the simulated system is, its number of atoms
N would be negligible compared with the number of atoms contained in a
macroscopic piece of matter (of the order of 1023 ), and the ratio between the
number of surface atoms and the total number of atoms would be much larger
2
Shifting eliminates the energy discontinuity, but not the force discontinuity. Some re-
searchers altered the form of the potential near Rc to obtain a smoother junction, but there
is no standard way to do that.

12
than in reality, causing surface e ects to be much more important than what
they should.
A solution to this problem is to use periodic boundary conditions (PBC).
When using PBC, particles are enclosed in a box, and we can imagine that
this box is replicated to in nity by rigid translation in all the three cartesian
directions, completely lling the space. In other words, if one of our particles is
located at position r in the box, we assume that this particle really represents
an in nite set of particles located at
r + `a + mb + nc ; (`; m; n = ;1; 1);
where `; m; n are integer numbers, and I have indicated with a; b; c the vec-
tors corresponding to the edges of the box. All these \image" particles move
together, and in fact only one of them is represented in the computer program.
The key point is that now each particle i in the box should be thought as
interacting not only with other particles j in the box, but also with their images
in nearby boxes. That is, interactions can \go through" box boundaries. In
fact, one can easily see that (a) we have virtually eliminated surface e ects from
our system, and (b) the position of the box boundaries has no e ect (that is, a
translation of the box with respect to the particles leaves the forces unchanged).
Apparently, the number of interacting pairs increases enormously as an ef-
fect of PBC. In practice, this is not true because potentials usually have a
short interaction range. The minimum image criterion discussed next simpli-
es things further, reducing to a minimum the level of additional complexity
introduced in a program by the use of PBC.
2.2.1 The minimum image criterion
Let us suppose that we are using a potential with a nite range: when separated
by a distance equal or larger than a cuto distance Rc , two particles do not
interact with each other. Let us also suppose that we are using a box whose
size is larger than 2Rc along each cartesian direction.
When these conditions are satis ed, it is obvious that at most one among
all the pairs formed by a particle i in the box and the set of all the periodic
images of another particle j will interact.
To demonstrate this, let us suppose that i interacts with two images j1 and
j2 of j. Two images must be separated by the translation vector bringing one
box into another, and whose length is at least 2Rc by hypothesis. In order to
interact with both j1 and j2, i should be within a distance Rc from each of
them, which is impossible since they are separated by more than 2Rc .
When we are in these conditions, we can safely use is the minimum image
criterion: among all possible images of a particle j, select the closest, and throw
away all the others. In fact only the closest is a candidate to interact; all the
others certainly do not.
This operating conditions greatly simplify the set up of a MD program, and
are commonly used. Of course, one must always make sure that the box size is
at least 2Rc along all the directions where PBCs are in e ect.

13
2.2.2 Geometries with surfaces
The purpose of PBCs is to eliminate surface e ects. However, we may also be
interested in situations where we want to have surfaces. Obviously, all surface
physics problems belong to this class!
For a surface simulation, the model usually adopted is that of the slab: a
thick slice of material, delimited by two free surfaces. This is simply obtained
by removing PBCs along one direction (usually taken to be z) while retaining
them in the orthogonal plane. Therefore, a slab must be thought as replicated
to in nity in the xy plane, but there is no replication along the slab normal z.
If the slab is thick enough, its inner part is expected to be quite similar
to the bulk of the material. The two surfaces (at the top and at the bottom
of the slab) can then be thought of as two decoupled, independent surfaces.
With this assumption, the system behavior would be close to that of a single
surface at the top of a semiin nite system|a condition which would be closer
to actual surface experiments but which is unfortunately impossible to realize
in simulation.
There are also circumstances where researchers nd it preferable to \freeze"
a few layers of material at one side of the slab|with the atoms constrained to
sit at perfect bulk-like crystal positions|leaving only the other side free. This
is done when it is believed that spurious e ects induced by the frozen side into
the \good" side of the slab are of smaller entity than the e ects induced by
another free surface at the same separation distance. This is to be preferred in
those cases where massive perturbations occur at the surface, such as atomic
rearrangements (surface reconstructions) or local melting. If a slab with a frozen
side is chosen, care must be taken to freeze a number of layers corresponding
to a thickness of at least Rc, to guarantee that no mobile atom would be able
to \see" the surface \through" the xed atoms.
One can also leave PBCs along only one direction only, thus obtaining a
wire geometry with the wire axis along the PBC direction, or remove them
altogether. This latter condition correspond to a cluster of atoms. Clusters are
important systems and simulation of clusters have been performed since the
early days of molecular dynamics and Monte Carlo.

2.3 Time integration algorithm

The engine of a molecular dynamics program is its time integration algorithm,
required to integrate the equation of motion of the interacting particles and
follow their trajectory.
Time integration algorithms are based on nite di erence methods, where
time is discretized on a nite grid, the time step
t being the distance between
consecutive points on the grid. Knowing the positions and some of their time
derivatives at time t (the exact details depend on the type of algorithm), the
integration scheme gives the same quantities at a later time t +
t. By iterating
the procedure, the time evolution of the system can be followed for long times.
Of course, these schemes are approximate and there are errors associated
with them. In particular, one can distinguish between

14
  Truncation errors, related to the accuracy of the nite di erence method
with respect to the true solution. Finite di erence methods are usually
based on a Taylor expansion truncated at some term, hence the name.
These errors do not depend on the implementation: they are intrinsic to
the algorithm.
 Round-o errors, related to errors associated to a particular implementa-
tion of the algorithm. For instance, to the nite number of digits used in
computer arithmetics.
Both errors can be reduced by decreasing
t. For large
t, truncation errors
dominate, but they decrease quickly as
t is decreased. For instance, the Verlet
algorithm discussed in x2.3.1 has a truncation error proportional to
t4 for each
integration time step. Round-o errors decrease more slowly with decreasing

t, and dominate in the small
t limit. Using 64-bit precision (corresponding to
\double precision" when using Fortran on the majority of today's workstations)
helps to keep round-o errors at a minimum.
Two popular integration methods for MD calculations are the Verlet algo-
rithm and predictor-corrector algorithms. They are quickly presented in the
sections below. For more detailed informations on time integration algorithms,
the reader is referred to refs. [3, 6] for a general survey, and to ref. [19] for a
deeper analysis.
2.3.1 The Verlet algorithm
In molecular dynamics, the most commonly used time integration algorithm is
probably the so-called Verlet algorithm [13]. The basic idea is to write two third-
order Taylor expansions for the positions r(t), one forward and one backward in
time. Calling v the velocities, a the accelerations, and b the third derivatives
of r with respect to t, one has:
r(t +
t) = r(t) + v(t)
t + (1=2)a(t)
t + (1=6)b(t)
t + O(
t )
2 3 4

r(t ;
t) = r(t) ; v(t)
t + (1=2)a(t)
t ; (1=6)b(t)
t + O(
t )

2 3 4

Adding the two expressions gives

r(t +
t) = 2r(t) ; r(t ;
t) + a(t)
t + O(
t )
2 4
(2.5)
This is the basic form of the Verlet algorithm. Since we are integrating Newton's
equations, a(t) is just the force divided by the mass, and the force is in turn a
function of the positions r(t):
a(t) = ;(1=m)rV (r(t)) (2.6)
As one can immediately see, the truncation error of the algorithm when
evolving the system by
t is of the order of
t4 , even if third derivatives do
not appear explicitly. This algorithm is at the same time simple to implement,
accurate and stable, explaining its large popularity among molecular dynamics
simulators.

15
 A problem with this version of the Verlet algorithm is that velocities are
not directly generated. While they are not needed for the time evolution, their
knowledge is sometimes necessary. Moreover, they are required to compute the
kinetic energy K , whose evaluation is necessary to test the conservation of the
total energy E = K + V . This is one of the most important tests to verify that a
MD simulation is proceeding correctly. One could compute the velocities from
the positions by using
v(t) = r(t +
t)2
; r(t ;
t) :
t
However, the error associated to this expression is of order
t2 rather than
t4 .
To overcome this diculty, some variants of the Verlet algorithm have been
developed. They give rise to exactly the same trajectory, and di er in what
variables are stored in memory and at what times. The leap-frog algorithm, not
reported here, is one of such variants [20] where velocities are handled somewhat
better.
An even better implementation of the same basic algorithm is the so-called
velocity Verlet scheme, where positions, velocities and accelerations at time
t +
t are obtained from the same quantities at time t in the following way:
r(t +
t) = r(t) + v(t)
t + (1=2)a(t)
t2
v(t +
t=2) = v(t) + (1=2)a(t)
t
a(t +
t) = ;(1=m)rV (r(t +
t))
v(t +
t) = v(t +
t=2) + (1=2)a(t +
t)
t
Note how we need 9N memory locations to save the 3N positions, velocities
and accelerations, but we never need to have simultaneously stored the values
at two di erent times for any one of these quantities.
2.3.2 Predictor-corrector algorithm
Predictor-corrector algorithms constitute another commonly used class of meth-
ods to integrate the equations of motion. Those more often used in molecular
dynamics (see e.g. [12]) are due to Gear, and consists of three steps:
1. Predictor. From the positions and their time derivatives up to a certain
order q, all known at time t, one \predicts" the same quantities at time
t +
t by means of a Taylor expansion. Among these quantities are, of
course, accelerations a.
2. Force evaluation. The force is computed taking the gradient of the poten-
tial at the predicted positions. The resulting acceleration will be in general
di erent from the \predicted acceleration". The di erence between the
two constitutes an \error signal".
3. Corrector. This error signal is used to \correct" positions and their deriva-
tives. All the corrections are proportional to the error signal, the coe-
cient of proportionality being a \magic number" determined to maximize
the stability of the algorithm.

16
 For more details the reader is referred to [3], x3.2, or [6], x4.4. A detailed
comparison between the Verlet and the Gear scheme can be found in [19].

17
Chapter 3
Running, measuring,
analyzing
In this chapter I discuss how we interact with a molecular dynamics program:
 how we start a simulation
 how we control a simulation, for instance to explore a phase diagram of a
substance
 how we extract results from a simulation
plus a couple of extra topics.

3.1 Starting a simulation

To start a simulation we need to de ne the MD box and the set of positions
and velocities to assign initially to the particles. There are two common ways
of doing this:
3.1.1 Starting from scratch
If we are starting from scratch, we have to \create" a set of initial positions
and velocities.
Positions are usually de ned on a lattice, assuming a certain crystal struc-
ture. This structure is typically the most stable at T=0 with the given potential.
Initial velocities may be taken to be zero, or from a Maxwellian distribution as
described below.
Such initial state will not of course correspond to an equilibrium condition.
However, once the run is started equilibrium is usually reached within a time
of the order of 100 time steps.
Some randomization must be introduced in the starting sample. If we do
not do so, then all atoms are symmetrically equivalent, and the equation of
motions cannot do anything di erent that evolving them in exactly the same
way. In the case of a perfect lattice, there would be no net force on the atoms
for symmetry reasons, and therefore atoms would sit idle inde nitely.

18
 Typical ways to introduce a random component are
1. small random displacements are added to the lattice position. The am-
plitude of these displacements must not be too large, in order to avoid
overlap of atomic cores. A few percent of the lattice spacing is usually
more than adequate.
2. the initial velocities are assigned taking them from a Maxwell distribution
at a certain temperature T. When doing this, the system will have a small
total linear momentum, corresponding to a translational motion of the
whole system. Since this is somewhat inconvenient to have, it is common
practice to subtract this component from the velocity of each particle in
order to operate in a zero total momentum condition.
Initial randomization is usually the only place where chance enters a molec-
ular dynamics simulation. The subsequent time evolution is completely deter-
ministic.
3.1.2 Continuing a simulation
Another possibility is to take the initial positions and velocities to be the nal
positions and velocities of a previous MD run.
This is in fact the most commonly used method in actual production. For
instance, if one has to take \measurements" at di erent temperatures, the stan-
dard procedure is to set up a chain of runs, where the starting point at each T
is taken to be the nal point at the preceding T (lower if we are heating, higher
if we are cooling).

3.2 Controlling the system

In a molecular dynamics simulation a system could be in a state characterized
by a certain density, temperature, pressure (the calculation of these quantities
is described below): the phase diagram of the simulated material can be inves-
tigated. However, how can we bring the system in the desired region? In other
words, how can we control the system?
In a standard calculation, the density is controlled by the choice of the box
volume V . Programs usually have provisions for rescaling the volume upon
request at the beginning of the calculation. Volume changes should be modest,
of the order of at most a few percent. The pressure will be measured during
the run, following the method described in x3.5.7. In other, more advanced
simulation schemes, the user chooses the pressure, and the volume of the box
is a variable to be measured (see x3.11).
Temperature changes are usually achieved by enabling a device in the code
which brings the system to a desired temperature by rescaling the velocities.
In the velocity Verlet algorithm discussed at the end of x2.3.1, this may be
accomplished by replacing the equation
v(t +
t=2) = v(t) + (1=2)a(t)
t

19
with s

v(t +
t=2) = TT(t) v(t) + (1=2)a(t)
t

where T is the desired temperature, and T (t) the \instantaneous temperature".
Such a modi cation means that we are no longer following Newton's equations,
and the total energy is no longer conserved. Important data should not be
collected in this stage: these \controlled temperature" simulations should be
used only to bring the system from one state to another. One should always
wait for the system to reach equilibrium under clean constant energy conditions
before collecting data.

3.3 Equilibration
Every time the state of the system changes, the system will be \out of equi-
librium" for a while. We are referring here to thermodynamic equilibrium. By
this, it is meant that the indicators of the system state|of which many are
described below|are not stationary (that is, uctuating around a xed value),
but relaxing towards a new value (that is, uctuating around a value which is
slowly drifting with time).
The state change may be induced by us or spontaneous. It is induced by us
when we change a parameter of the simulation|such as the temperature, or the
density|thereby perturbing the system, and then wait for a new equilibrium
to be reached. It is spontaneous when, for instance, the system undergoes a
phase transition, thereby moving from one equilibrium state to another.
In all cases, we usually want equilibrium to be reached before starting per-
forming measurements on the system. A physical quantity A generally ap-
proaches its equilibrium value exponentially with time:
A(t) = A + C exp(;t= ) (3.1)
where A(t) indicates here a physical quantities averaged over a short time to get
rid of instantaneous uctuations, but not of its long-term drift. The relevant
variable is here the relaxation time  . We may have cases where  is of the order
of hundreds of time steps, allowing us to see A(t) converge to A and make a
direct measurement of the equilibrium value. In the opposite case,  could be
much larger than our overall simulation time scale; for instance, of the order
of one second. In this case, we do not see any relaxation occurring during the
run, and molecular dynamics hardly seems a valid technique to use in such a
situation. In intermediate situations, we may see the drift but we cannot wait
long enough to observe convergency of A(t) to A . In these cases, one can often
obtain an estimate for A by applying (3.1) on the available data, even if the
nal point is still relatively far from it.

3.4 Looking at the atoms

The simplest|at least conceptually|kind of \probe" that we can apply to our
system to understand what is going on is \looking at it". One can assign a

20
radius to the atoms (for instance corresponding to half the distance at which
the pair potential between two atoms becomes strongly repulsive), represent
the atoms as balls having that radius, and have a plotting program construct a
\photograph" of the system from a certain point of observation. The resulting
image can be observed on a graphic screen or printed.
Programs capable of doing this are commonly available. One of them, called
BallRoom1 , has been written by the author and its source code is freely avail-
able. Among other things, this program allows to map properties attached to
individual atoms (for instance, di usion, potential energy, coordination, etc)
to colors, providing a vivid representation of the spatial distribution of these
properties.
In some cases, such as in the analysis of complex structural arrangements
of atoms, visual inspection is an invaluable tool. In other cases, however, it
gives little information, and useful data must be extracted by processing atomic
coordinates and/or velocities, obtaining for instance correlation functions, as
discussed in the following.

3.5 Simple statistical quantities to measure

Measuring quantities in molecular dynamics usually means performing time
averages of physical properties over the system trajectory. Physical properties
are usually a function of the particle coordinates and velocities. So, for instance,
one can de ne the instantaneous value of a generic physical property A at time
t:
A(t) = f (r1 (t); : : : ; rN (t); v1 (t); : : : ; vN (t))
and then obtain its average
1 NT
hAi = N
X
A(t)
T t=1
where t is an index which runs over the time steps from 1 to the total number
of steps NT .
There are two equivalent ways to do this in practice:
1. A(t) is calculated at each time step by the MD program while running.
The sum t A(t) is also updated at each step. At the end of the run the
P

average is immediately obtained by dividing by the number of steps. This

is the preferred method when the quantity is simple to compute and/or
particularly important. An example is the system temperature.
2. Positions (and possibly velocities) are periodically dumped in a \trajec-
tory le" while the program is running. A separate program, running
after the simulation program, processes the trajectory and computes the
desired quantities. This approach can be very demanding in terms of disk
space: dumping positions and velocities using 64-bit precision takes 48
1
http://www.sissa.it/furio/ballroom.html

21
 bytes per step and per particle. However, it is often used when the quan-
tities to compute are complex, or combine di erent times as in dynamical
correlations, or when they are dependent on other additional parameters
that may be unknown when the MD program is run. In these cases, the
simulation is run once, but the resulting trajectory can be processed over
and over.
The most commonly measured physical properties are discussed in the fol-
lowing.
3.5.1 Potential energy
The average potential energy V is obtained by averaging its instantaneous value,
which is usually obtained straightforwardly at the same time as the force com-
putation is made. For instance, in the case of two-body interactions
 (jri (t) ; rj (t)j)
XX
V (t) =
i j>i
Even if not strictly necessary to perform the time integration (forces are all is
needed), knowledge of V is required to verify energy conservation. This is an
important check to do in any MD simulation.
3.5.2 Kinetic energy
The instantaneous kinetic energy is of course given by
K (t) = 12 mi [vi (t)]2
X

i
and is therefore extremely easy to compute. We will call K its average on the
run.
3.5.3 Total energy
The total energy E = K + V is a conserved quantity in Newtonian dynamics.
However, it is common practice to compute it at each time step in order to
check that it is indeed constant with time. In other words, during the run
energy ows back and forth between kinetic and potential, causing K(t) and
V(t) to uctuate while their sum remains xed.
In practice there could be small uctuations in the total energy, in a typical
amount of, say, one part in 104 or less. These uctuations are usually caused
by errors in the time integration (see 2.3), and can be reduced in magnitude by
reducing the time step if considered excessive. Ref. [19] contains an in-depth
analysis of total energy uctuations using various time integration algorithms.
Slow drifts of the total energy are also sometimes observed in very long
runs. Such drifts could also be originated by an excessive
t. Drifts are more
disturbing than uctuations because the thermodynamic state of the system is
also changing together with the energy, and therefore time averages over the

22
run do not refer to a single state. If drifts over long runs tend to occur, they
can be prevented, for instance by breaking the long run into smaller pieces and
restoring the energy to the nominal value between one piece and the next. A
common mechanism to adjust the energy is to modify the kinetic energy via
rescaling of velocities.
A nal word of caution: while one may be tempted to achieve \perfect"
energy conservation by reducing the time step as much as desired, working
with an excessively small time step may result in waste of computer time. A
practical compromise would probably allow for small energy uctuations and
perhaps slow energy drifts, as a price to pay to work with a reasonably large

t. See also refs. [3], x3.5, and [6], x4.4.4.
3.5.4 Temperature
The temperature T is directly related to the kinetic energy by the well-known
equipartition formula, assigning an average kinetic energy kB T=2 per degree of
freedom:
K = 23 NkB T (3.2)
An estimate of the temperature is therefore directly obtained from the aver-
age kinetic energy K (see x3.5.2). For practical purposes, it is also common
practice to de ne an \instantaneous temperature" T(t), proportional to the
instantaneous kinetic energy K(t) by a relation analogous to the one above.
3.5.5 The caloric curve
Once the total energy E (also called internal energy) and the temperature T
are measured in di erent runs corresponding to di erent thermodynamic states,
one can construct the caloric curve E(T). This is a useful tool to monitor the
occurrence of phase transitions, which imply a jump of E (T ) if rst-order, or
in a derivative of E (T ) if second- or higher order. This strategy, however, is
not suitable for a reliable estimate of the melting temperature Tm , as discussed
later in x3.6.
The most common rst-order transition is melting, easily observable by
simulation. When the system abandons the crystalline structure and becomes
a disordered liquid, we observe a jump of E (T ), corresponding to the latent heat
of fusion. This usually happens at a temperature somewhat higher than the true
melting temperature Tm of the model, due to hysteresis e ects associated with
the necessity to wait for a seed of the liquid phase to appear by a spontaneous
uctuation. Only when a liquid seed containing a few atoms is formed, the
liquid phase can start to grow at the expense of the solid one. On the typical
molecular dynamics time scale, one has to \overshoot" and set the temperature
20-30% above Tm to see melting happening.

23
3.5.6 Mean square displacement
The mean square displacement of atoms in a simulation can be easily computed
by its de nition
MSD = hjr(t) ; r(0)j2 i (3.3)
where h: : :i denotes here averaging over all the atoms (or all the atoms in a
given subclass). Care must be taken to avoid considering the \jumps" of par-
ticles to refold them into the box when using periodic boundary conditions as
contributing to di usion.
The MSD contains information on the atomic di usivity. If the system is
solid, MSD saturates to a nite value, while if the system is liquid, MSD grows
linearly with time. In this case it is useful to characterize the system behavior
in terms of the slope, which is the di usion coecient D:
D = tlim 1 hjr(t) ; r(0)j2 i (3.4)
!1 6t
The 6 in the above formula must be replaced with 4 in two-dimensional systems.
3.5.7 Pressure
The measurement of the pressure in a molecular dynamics simulation is based
on the Clausius virial function
N
W (r1 ; : : : ; rN ) = ri
Fi
X
TOT
(3.5)
i=1
where FTOT
i is the total force acting on atom i. Its statistical average hW i will
be obtained, as usual, as an average over the molecular dynamics trajectory:
1 Z t N
hW i = tlim ri( )
miri( ):
X

!1 t
d
0
i=1
where use has been made of Newton's law. By integrating by parts,
1 Z t N
hW i = ; tlim mijr_ i ( )j2 :
X

!1 t d
0
i=1
This is twice the average kinetic energy, therefore by the equipartition law of
statistical mechanics
hW i = ;DNkB T (3.6)
where D is the dimensionality of the system (2 or 3), N the number of particles,
kB the Boltzmann constant.
Now, one may think of the total force acting on a particle as composed of
two contributions:
FTOT
i = Fi + FEXT
i (3.7)
where Fi is the internal force (arising from the interatomic interactions), and
FEXT
i is the external force exerted by the container's walls. If the particles

24
are enclosed in a parallelepipedic container of sides Lx , Ly , Lz , volume V =
Lx Ly Lz , and with the coordinates origin on one of its corners, the part hW EXT i
due to the container can be evaluated using the de nition (3.5):
hW i = Lx(;PLy Lz ) + Ly (;PLxLz ) + Lz (;PLxLy ) = ;DPV
EXT

where ;PLy Lz is, for instance, the external force Fx applied by the yz wall
EXT

along the x directions to particles located at x = Lx , etc. Eq. (3.6) can then
be written *
N +

ri
Fi ; DPV = ;DNkB T
X

i=1
or *
N +
1
PV = NkB T + D
X
ri
Fi (3.8)
i=1
This important result is known as the virial equation. All the quantities except
the pressure P are easily accessible in a simulation, and therefore (3.8) consti-
tutes a way to measure P . Note how (3.8) reduces to the well-known equation
of state of the perfect gas if the particles are non-interacting.
In the case of pairwise interactions via a potential (r), it is left as an
exercise to the reader to verify that (3.8) becomes
* +

PV = NkB T ; D1 rij d
XX

i j>i dr rij : (3.9)

This expression has the additional advantage over (3.8) to be naturally suited
to be used when periodic boundary conditions are present: it is sucient to
take them into account in the de nition of rij .

3.6 Measuring the melting temperature

The behavior of the mean square displacement as a function of time easily
allows to discriminate between solid and liquid. One might be then tempted
to locate the melting temperature Tm of the simulated substance by increasing
the temperature of a crystalline system until di usion appears, and the caloric
curve exhibits a jump indicating absorption of latent heat.
While these are indeed indicators of a transition from solid to liquid, the
temperature at which this happen in a MD simulation is invariably higher
than the melting temperature. In fact, the melting point is by de nition the
temperature at which the solid and the liquid phase coexist (they have the
same free energy). However, lacking a liquid seed from where the liquid could
nucleate and grow, overheating above melting commonly occurs. In this region
the system is in a thermodynamically metastable state, nevertheless it appears
stable within the simulation time. An overheated bulk crystal breaks down
when its mechanical instability point is reached. This point may correspond to
the vanishing of one of the shear moduli of the material or to similar instabilities,
and is typically larger than Tm by an amount of the order of 20-30%.

25
 While mechanical instability is certainly useful to estimate roughly the lo-
cation of Tm , more precise methods exist. One of them consists of setting up
a large sample consisting approximately of 50% solid and 50% liquid. Such a
sample could be initially constructed arti cially2, and it will contain interface
regions between solid and liquid. One then tries to establish an equilibrium
state where solid and liquid can coexist. When this goal is achieved, the tem-
perature of the state has to be Tm by de nition.
To this end, the system is evolved microcanonically at an energy E believed
to be in proximity of the melting jump in the caloric curve. Let us call T the
initial temperature assumed by the system. This temperature is an increasing
but unknown function of the energy E, and it will in general di er from Tm .
Suppose that T > Tm . Then, the liquid is favored over the solid, and the solid-
liquid interfaces present in the system will start to move in order to increase
the fraction of liquid at the expense of the solid. In other words, some material
melts. As this happens, the corresponding latent heat of melting is absorbed by
the system. Since the total energy is conserved, absorption of latent heat auto-
matically implies a decrease of the kinetic energy. Therefore, the temperature
goes down. The force driving the motion of the solid-liquid interface will also
decrease, and as time goes on, the position of the interface and the temperature
will exponentially reach an equilibrium state. The temperature will then be Tm .
If T < Tm , the inverse reasoning applies: latent heat is released and con-
verted into kinetic energy, and temperature again converges exponentially to
Tm , this time from below.
The melting temperature depends on the pressure of the system. A pressure
measurement on the nal state is also required to characterize it thermodynam-
ically. More often, the above procedure is carried out using constant pressure
techniques (see x3.11). In this case, the volume of the box automatically changes
as material melts or crystallizes, to accommodate the di erence in density be-
tween the two phases at constant pressure.

3.7 Real space correlations

Real space correlation functions are typically of the form
hA(r)A(0)i
and are straightforward to obtain by molecular dynamics: one has to compute
the quantity of interest A(r) starting from the atomic positions and velocities for
several con gurations, construct the correlation function for each con guration,
and average over the available con gurations.
The simplest example is the pair correlation function g(r), which is essen-
tially a density-density correlation. g(r) is the probability to nd a pair a
distance r apart, relative to what expected for a uniform random distribution
A possible way to do this is to \freeze" half of the particles (no motion allowed), raise the
2

temperature in order to melt the other half, and then release the frozen particles.

26
of particles at the same density:
*
N N +

g(r) = N1 (r ; rij )

X X
(3.10)
i=1 j =1
(j 6=i)

This function carries information on the structure of the system. For a crystal,
it exhibits a sequence of peaks at positions corresponding to shells around a
given atom. The positions and magnitude of the peaks are a \signature" of the
crystal structure (fcc, hcp, bcc, . . . ) of the system. For a liquid, g(r) exhibits
its major peak close to the average atomic separation of neighboring atoms,
and oscillates with less pronounced peaks at larger distances. The magnitude
of the peaks decays exponentially with distance as g(r) approaches 1. In all
cases, g(r) vanishes below a certain distance, where atomic repulsion is strong
enough to prevent pairs of atoms from getting so close.
One quantity often computed from g(r) is the average number of atoms
located between r1 and r2 from a given atom,
Z r2
 g(r)4r2 dr:
r1
This allows to de ne coordination numbers also in situations where disorder is
present.
The calculation of g(r) is intrinsically a O(N 2 ) operation, and therefore
it can slow down considerably an optimized molecular dynamics program. If
the behavior at large r is not important, it might be convenient to de ne a
cuto distance, and use techniques borrowed from fast force calculations using
short-ranged potentials to decrease the computational burden. It should also be
noted that periodic boundary conditions impose a natural cuto at L=2, where
L is the minimum between the box sizes Lx, Ly , Lz in the three directions. For
larger distance the results are spoiled by size e ects.

3.8 Reciprocal space correlations

Structural informations can also be collected by working in reciprocal space.
The most basic quantity is the space Fourier transform of the density (r) =
i r ; ri ):
 (
P

N
(k) = eik
ri
X
(3.11)
i=1
This quantity, easily and quickly obtained for a given con guration, is the build-
ing block to construct the static structure factor
S (k) = N1 hD(k)(;k)i E
= N1 D ij eik
(ri ;rj ) E
P
(3.12)
= N ij cos(k
rij )
P
1

This quantity is extremely useful for comparisons with results from scattering
experiments. In liquids, which are isotropic, S (k) depends only on the modulus

27
of the wavevector, and is essentially a Fourier transform of the pair correlation
function:
S (k) = 1 + g(k):
The computation of S(k) via (3.12) is however much faster than that of g(r).
By replacing (k) with some other generic observable A(k), any other cor-
relation function in reciprocal space can be computed.
It should be noted that, working with periodic boundary conditions, results
could not depend on the choice of a particular particle image. That is, one
must impose that (choosing for instance direction x) eikx (xi +Lx) = eikx xi where
Lx is the length of the repeated box. This implies a restriction in the choice of
the allowable wavevectors:
k = 2 n
x Lx x
where nx is an integer, and similarly along y and z . Therefore, the allowed
wavevectors are quantized as a consequence of adopting periodic boundary con-
ditions.

3.9 Dynamical analysis

One strength of molecular dynamics with respect to other methods such as
Monte Carlo is the fact that the time evolution of the system is followed and|
therefore|information on the dynamics of the system is fully present. Experi-
mental data (for instance, phonon dispersion curves in the case of crystals) are
usually available in the wavevector-frequency space. Information obtained by
MD in real space and real time will therefore be Fourier-transformed with re-
spect to both space and time for useful comparisons with experimental results.
To perform dynamical analysis, one typically instructs the simulation program
to dump periodically on a mass storage medium the positions and possibly the
velocities of all the particles. These can easily produce very large les: the cal-
culation of dynamical quantities requires plenty of room for data. These data
are then analyzed later by other programs to extract the quantities of interest.
In some cases, dynamical analysis are directly performed by the MD program
while it is running, thus saving disk space.
One typical quantity is the Van Hove correlation function G(r; t), represent-
ing the probability to nd an atom at position r at time t if there was an atom
at position r = 0 at time t = 0. This quantity is connected by a double Fourier
transform to the dynamic structure factor S (k; !) reported by experiments:
Z Z

S (k ; ! ) = dr dt ei(k
r;!t) G(r; t): (3.13)

As for static quantities, the most convenient way to proceed in terms of com-
putational e ort is that of using the space Fourier transform of the density as
a building block:
(k; t) = eik
ri(t)
X
(3.14)
i
This is a single-particle quantity which is very fast to compute. From here, one
constructs a time-dependent density-density correlation function (also called

28
intermediate scattering function) by performing an average over a molecular
dynamics trajectory:
F (k; t) = 1 h(k; t + t )(;k; t )i :
N (3.15)
The average is made for all the available choices for t . The maximum t (con-
nected with the frequency resolution) is clearly limited by the length of the run,
while the periodicity with which data are collected|typically a few time steps|
controls the maximum frequency in the nal spectra. The dynamic structure
factor is nally obtained by a time Fourier transform:
Z

S (k; !) = dt e;i!t F (k; t) (3.16)

Again, the allowed k are quantized as a consequence of periodic boundary
conditions. S (k; !) will exhibit peaks in the (k; !) plane, corresponding to
the dispersion of propagating modes. The broadening of these peaks is related
with anharmonicity. By performing calculations at di erent temperatures, one
can for instance study the e ect of temperature on the phonon spectra, fully
including anharmonic e ects: in this sense, the approach is clearly superior to
standard lattice dynamics techniques. The main inconvenience of MD is that
it is not so straightforward to extract eigenvectors.
When dealing with vibrational modes and di erent polarizations, it is often
convenient to follow the approach outlined above, but using currents
j(k; t) = vi (t)eik
ri t
X
( )
(3.17)
i
in place of the density (k; t).

3.10 Annealing and quenching: MD as an optimiza-

tion tool
Molecular dynamics may also have a role as on optimization tool. Let us sup-
pose that a set of N particles has many possible equilibrium con gurations.
The energy of these con gurations is in general di erent, and one of them will
be the optimal one; all of them, however, correspond to local minima in the en-
ergy, and are each other separated by energy barriers. Such a situation occurs
very commonly, for instance, in cluster physics.
Finding the optimal structure within an approach based on traditional mini-
mization techniques (steepest descent method, conjugate gradient method, etc.)
is tricky, as these methods do not normally overcome energy barriers and tend
to fall into the nearest local minimum. One therefore would have to try out
several di erent starting points, corresponding to di erent \attraction basins"
in the energy landscape, and relax each of them to the bottom of the basin.
The optimal structure would then be the one with the lowest energy, provided
that we were sage enough to select it in the list of candidates to try.

29
 Temperature in a molecular dynamics (or Monte Carlo) calculation pro-
vides a way to \ y over the barriers": states with energy E are visited with a
probability exp(;E=kB T ). If T is suciently large, the system can \see" the
simultaneous existence of many di erent minima, spending more time in the
deeper ones. By decreasing slowly T to 0, there is a good chance that the system
will be able to pick up the best minimum and land into it. This consideration
is the base of simulated annealing methods, where the system is equilibrated
at a certain temperature and then slowly cooled down to T = 0. While this
procedure does not guarantee that the true minimum will be reached, it often
brings the system into a good minimum, and as no a priori assumptions are
made about the optimal structure, it often brings about structures which were
dicult to foresee by intuition alone.
This method is often used to optimize atomic structures, but its validity
is more general: given an objective function Z ( 1 ; : : : ; N ) depending on N
parameters, one can regard each of these parameters as a degree of freedom,
assign a \mass' to it, and move the system with a molecular dynamics or Monte
Carlo algorithm to perform simulated annealing. One of the early application
of this method can be found in a famous paper discussing an application to the
\traveling salesman" problem [21].

3.11 Other statistical ensembles

We have discussed so far the standard molecular dynamics scheme, based on
the time integration of Newton's equations and leading to the conservation of
the total energy. In the statistical mechanics parlance, these simulations are
performed in the microcanonical ensemble, or NVE ensemble: the number of
particles, the volume and the energy are constant quantities.
The microcanonical average of a physical quantity A is obtained as a time
average on the trajectory:
NT
hAi = N1
X
NVE A(;(t)) (3.18)
T t=1
where I denote with ;(t) the phase space coordinate of the system (3N positions
and 3N velocities).
There are other important alternatives to the NVE ensemble, that we will
mention here only brie y. The basic idea is that of integrating other equations
in place of Newton's equations, in such a way that sampling is performed in an-
other statistical ensemble. Averages of physical quantities in the new ensemble
will be again obtained as time averages, similarly to eq. (3.18).
A scheme for simulations in the isoenthalpic-isobaric ensemble (NPH) has
been developed by Andersen [22]. Here, an additional degree of freedom repre-
senting the volume of the box has been introduced, and all the particle coordi-
nates are given in units relative to the box. The volume V of the box becomes
a dynamical variable, with a kinetic energy and a potential energy which just
PV , where P is the external pressure. The enthalpy H = E + PV is a conserved
quantity.

30
 Parrinello and Rahman [23] developed a variant where the shape of the box
can vary as well as the volume. This is achieved by introducing 9 new degrees of
freedom instead of 1: the components of the three vectors spanning the MD box.
Each of them is a new dynamical variable, evolving accordingly to equation of
motion derived from an appropriate Lagrangian. This scheme allows to study
structural phase transitions as a function of pressure, where for example the
system abandons a certain crystal structure in favor of a more compact one.
Another very important ensemble is the canonical ensemble (NVT). In a
method developed by Nose and Hoover [24], this is achieved by introducing a
time-dependent frictional term, whose time evolution is driven by the imbal-
ance between the instantaneous kinetic energy and the average kinetic energy
(3N=2)kB T .

31
Chapter 4
Interatomic potentials
In molecular dynamics forces are derived from a potential energy function V,
which depend on the particle coordinates:
Fi = ;rV (r1 ; : : : ; rN ) (4.1)
The problem of modelling a material can therefore be restated as that of nding
a potential function V (r1 ; : : : ; rN ) for that material.
Can such a function exist? And if yes, how is it related with the behavior
of the real material, which we know to be controlled by the law of quantum
mechanics rather than classical mechanics, and where the electrons certainly
play the major role in determining the bonding properties of the system? And
perhaps the most important question: how to nd a potential in practice?

4.1 The Born-Oppenheimer approximation

What do we mean exactly when we talk about \interatomic potentials"? A
system of interacting atoms is really made up of nuclei and electrons which
interact with each other. The true Hamiltonian for this system may be written
as
H = 2PMi + 2pmn + 21 jRZi Z;j R e2 + 1 X e2 ;X Zi e2 (4.2)
X 2 X 2 X

i i n ij i j j 2 nn jrn ; rn j in jRi ; rn j
0
0

where indexes i, j run on nuclei, n and n0 on electrons, Ri and Pi are positions

and momenta of the nuclei, rn and pn of the electrons, Zi the atomic number
of nucleus i, Mi its mass and m the electron mass. One clearly recognizes the
kinetic energy terms and the various coulombic interactions. In principle, one
should solve a Schrodinger equation for the total wavefunction (Ri ; rn ) and
then everything about the system is known.
Of course, this is impossible to carry out in practice, and approximation
schemes have to be employed. In 1923, Born and Oppenheimer noted that
nuclei are much heavier than electrons, and move on a time scale which is
about two order of magnitude longer than that of the electrons:
s
!el M
!nuc  m  100 (4.3)

32
It is therefore sensible to regard the nuclei as xed as far as the electronic part
of the problem is concerned, and factorize the total wavefunction as
(Ri ; rn ) = (Ri )(rn ; Ri ) (4.4)
where (Ri ) describes the nuclei, and (rn ; Ri ) the electrons (depending para-
metrically on the positions of the nuclei). With this assumptions, the problem
is reformulated in terms of two separate Schrodinger equations:
Hel(rn; Ri ) = V (Ri )(rn; Ri ) (4.5)
where
Hel =
X p2n + 1 X ZiZj e2 + 1 X e2 ; X Zi e2 ; (4.6)
n 2m 2 ij jRi ; Rj j 2 nn jrn ; rn j in jRi ; rn j
0
0

and " #
X Pi2 + V (R ) (R ) = E (R ): (4.7)
2Mi i i i
i
Eq. (4.5) is the equation for the electronic problem, considering the nuclei as
xed. The eigenvalue of the energy V (Ri ) will depend parametrically on the
coordinates of the nuclei; we call this quantity the interatomic potential. Once
found, this quantity enters (4.7) which will give the motion of the nuclei. Note
how in this equation there are no electronic degrees of freedom: all the electronic
e ects are incorporated in V (Ri ). It is customary to replace this Schrodinger
equation with a Newton equation, that is, to move the nuclei classically. We
have already discussed in x1.5.1 the conditions under which this can be done
safely.

4.2 The design of potentials

The above section de nes rigorously a potential, but it is not of much practical
help: solving (4.5) is still a formidable task and constitutes in fact a large
fraction of the activity done in condensed matter physics.
One possibility is to adopt some approximation, but really solve this equa-
tion: we then have rst-principles molecular dynamics, brie y mentioned later
in x4.7. While this is feasible, it requires massive computer resources, and poses
severe limits on the maximum size of the system and on the simulation time.
In these notes I will rather focus on the traditional approach where we get rid
completely of electronic degrees of freedom and move our nuclei (atoms) ac-
cording to some potential function V (Ri ) whose analytical form we specify in
advance. Our goal is to select functional forms which mimic the behavior of the
\true" potential in realistic ways for speci c materials.
Constructing a potential involves two steps:
1. Selecting an analytical form for the potential. In the past this was almost
invariably a sum of pairwise terms, with the energy of a pair depending
on their relative distance. Today novel many-body forms are tried in

33
 the attempt to capture as much as possible the physics and chemistry
of the bonding. A typical analytical form is constituted by a number of
functions, depending on geometrical quantities such as distances or angles,
or on intermediate variables such as atom coordinations.
2. Finding an actual parametrization for the functions that constitute the
analytical form we have chosen. This step is very important and can be
technically elaborate.
A variety of techniques have been utilized over the years to this end. At one
extreme, some groups tried to start from a rst principle description (that is,
where the electronic structure is kept into account), and obtain an expression
for the energy as a function of the nuclei position by means of successive approx-
imations. At the other extreme, other groups chose to t the potential (that
is, to parametrize the functions which appear in its analytical form) on exper-
imental data, giving maximum priority to realism rather than to connections
with rst-principles.
In all cases, potentials are designed with a \range of applicability" in mind.
Due to the vast di erences in the electronic structure, it would probably be
too ambitious to try modeling a bulk metal and a diatomic molecule of the
same element with the same potential: the environment is dramatically di er-
ent. However, it could be feasible to model simultaneously a bulk and a surface
environment, where the environment di ers (due to the reduced coordination
of the atoms at the surface) but not as dramatically as in the previous ex-
ample. The ability of a potential to work properly in di erent environments
is called transferability. When using a potential, the simulator should always
be familiar with its transferability properties, and look critically at results ob-
tained in unusual conditions|for example, very low coordinations, or very high
temperature, or very high pressure.

4.3 The problems with two-body potentials

In x2.1.1 I presented the Lennard-Jones potential, which is a typical pairwise
potential and probably the most commonly used one. In spite of that, the class
of materials which can be realistically modeled using this approach is in practice
limited to rare gases, where no electrons are available for bonding and atoms
are attracted with each other only through the weak van der Waals forces.
Systems of more practical interest such as metals and semiconductors cannot
be modeled with pairwise forces. Considering for instance noble metals, one
can easily identify a few indicators of many-atom e ects, summarized in the
following table, where experimental data for a few metals are compared with
Lennard-Jones data (but other pair potentials yield similar results):
Property Cu Ag Pt Au LJ
Ec=kB Tm 30 28 33 33 13
Ev =Ec 0.33 0.36 0.26 0.25  1
C12 =C44 1.5 1.9 3.3 3.7 1

34
 Ec =kB Tm is the ratio between the cohesive energy and the melting temper-
ature. This ratio is about 30 in metals, and about 10 in two-body systems.
This result indicates that metals exhibit some \extra cohesion" with respect
to pairwise systems, which is less e ective than two-body forces in keeping the
system in the crystalline state.
Ev =Ec is the ratio between the vacancy formation energy and the cohesive
energy. This number is between 1/4 and 1/3 in metals, but about 1 in two-
body systems (exactly 1 if relaxations are neglected, actually). When forming a
vacancy in a crystal structure with coordination number Z , one has to pay the
energy Ev to decrease coordination from Z to Z ; 1 for Z atoms. In contrast,
Ec is the energy to pay to decrease the coordination of a single atom from Z
to 0. In a two-body model, where a xed energy contribution is attached to
bonds between pairs of atoms, these two energies are the same as they are both
associated to the breaking of Z bonds. But this is not what happens in metals!
C12 =C44 is the ratio between two elastic constants of a cubic crystal (all
the systems in the table are fcc). This ratio is exactly 1 in two-body systems:
this is a so-called Cauchy relation that can be demonstrated analytically. But
deviations in metals are very common. The high value in Au is to be related
with the well-known high ductility and malleability of gold.
If one consider semiconductors, deviations from a two-body behavior are
even worse. For instance, silicon undergoes a series of structural phase tran-
sitions (from tetrahedral to -tin to simple cubic to fcc) under pressure: this
indicates that the energy di erence between this structures is not too large.
In other words, the cohesive energy is nearly independent upon coordination,
while a two-body model should favor the more packed structures, which have
more bonds.
Due to these and other [25] shortcomings, since the mid '80s researchers
started to gure out how to improve realism by incorporating many-atom e ects
in potentials.

4.4 Many-body potentials for metals

A signi cant progress was made during the '80s by the development of many-
atom potentials for metals based on the concept of density, or coordination, as
the key variable to consider. The main physical point to model is that bonds
become weaker when the local environment becomes more crowded: another
consequence of the Pauli principle. So, a plot of the cohesive energy as a function
of coordination should not be decreasing linearly as in two-body systems, but
should exhibit a positive curvature: decreasing faster when coordination is low,
and more slowly as coordination increases.
A possible form for the attractive part of the potential (the repulsive one
being still conveniently treated by a pairwise law) can be qualitatively obtained
by working out a connection with the tight-binding formalism [26]. The result
of this (non-rigorous) argument suggests a form for the energy of an atom i
s

Ei / h2ij / Zi
X p
(4.8)
j

35
where hij = h(rij ) = hijH jj i are the overlap integrals between i and its neigh-
bors. In the tight-binding formalism the orbitals are localized and these func-
tions vanish beyond a certain cuto distance.
The key result contained in (4.8) is that the energy is proportional to the
square root of the coordination, rather than to the coordination itself as in
two-body models. One can easily verify that such a scheme goes in the right
direction in solving the problems associated with pairwise potentials in metals.
For instance, p p
Ev  Z [ Z ;p Z ; 1]  1 ; (4.9)
E c Z 2
and one can also verify that the Cauchy relation C12 = C44 no longer holds.
Following this and other considerations, several schemes to construct many-
atom potentials for metals were developed, all essentially based on an analytical
form N N
V = 21
X X
(rij ) + U (ni ) (4.10)
i;j =1 i=1
(j 6=i)

where (r) is a two-body part, and U (n) is a function giving the energy of an
atom as a function of a \generalized coordination" n. n for a given atom is in
turn constructed as a superpositions of contributions from neighboring atoms:
XN
ni = (rij ) (4.11)
j =1
j6 i
( = )

where (r) is a short-ranged, decreasing function of distance. Belonging to

this scheme are the glue model [25], the embedded atom method [27] and the
Finnis-Sinclair potentials [28]. Also similar is the e ective medium theory [29].
Even if sharing the analytical form, these scheme di er vastly in the procedure
used to build the three functions (r), (r), U (n) constituting the model, often
resulting in rather di erent parametrizations for the same material. Fitting is
in fact the true crux of the matter of building potentials, and some groups use
molecular dynamics simulation to test trial potentials (for example, to evaluate
thermal expansion and melting point) during the construction stage.
It is interesting to note that eq. (4.10) is invariant with respect to the
transformation
^(r) = (r) + 2(r) (4.12)
U^ (n) = U (n) ; n
for any value of . This implies that there is no unique choice for (r) and U (n),
and one can arbitrarily impose a condition such as U 0 (n ) = 0 (where n is,
for instance, the coordination of a bulk atom at equilibrium at T = 0). It also
shows that if U (n) is a linear function of n, then the whole scheme reduces to
a two-body potential: many-body e ects are related to the curvature of U (n).
In fact, U 00 (n) is invariant with respect to (4.12).

36
 To implement these potentials in a molecular dynamics program, one has
to evaluate the forces
Fi = ; 0 (rij ) + [U 0 (ni) + U 0 (nj )]0 (rij ) rrij
X 
(4.13)
j 6=i ij
This calculation is only slightly more complex than that required for a two-body
system. Energy and forces can still be obtained by using pair distances as the
only input: no angular term or other 3-body or 4-body terms are present. This
allows to write very fast MD programs. On the other hand, the lack of true
angular forces makes it dicult to use these schemes to model metals where
covalent e ects in the bonding are important, for instance transition metals.
New, more powerful schemes are being developed to handle these materials.
The interested reader is invited to check [26] for further details.

4.5 Many-body potentials for semiconductors

Semiconductors are even more challenging than metals. Take for instance sili-
con, and consider that
 the most stable phase in the absence of pressure has the diamond struc-
ture: very open, with a coordination number of only 4;
 when a pressure is applied, new structures with increasing coordinations
appear in sequence ( -tin, simple cubic, fcc), indicating that they are not
so far in energy;
 the liquid is a metal, and is more dense than the solid (anomalous).
These features alone guarantee that constructing a potential for Si is not a
trivial task. Nevertheless, due to its technological importance, many groups
have challenged themselves in this territory, with various degrees of success.
4.5.1 The Stillinger-Weber potential
The Stillinger-Weber potential [30] is one of the rst attempts to model a semi-
conductor with a classical model. It is based on a two-body term and a three-
body term:
1 X X


V = 2 (rij ) + g(rij )g(rik ) cos jik + 3 1 2

(4.14)
ij ijk
where jik is the angle formed by the ij bond and the ik bond, and g(r) is
a decaying function with a cuto between the rst- and the second-neighbor
shell. The intent is obvious: favor those con guration where cos jik = ;1=3,
that is, where angles are as close as possible to those found in the diamond-
like tetrahedral structure, and make this structure more stable than compact
structure, as it should.
This potential gives a fairly realistic description of crystalline silicon. How-
ever, its built-in tetrahedral bias creates problems in other situations: it cannot

37
predict the right energies of the non-tetrahedral polytypes found under pres-
sure; the coordination of the liquid is too low; surface structures are not correct.
These are called transferability problems: it is dicult to use the potential un-
der conditions di erent from those it was designed for. Recent research [26]
shows that, to build more realistic models, one should consider analytic forms
which take into account the concept of local environment, with bond strengths
depending on it.
Nevertheless, potentials based on similar concepts (using only geometrical
quantities such as distances and bond angles as variables) are also currently
used in simulations of organic systems, where they are usually called force
elds. In spite of the limited transferability, which makes it dicult to change
the topology of a molecule by breaking and reforming bonds, or to work at
high temperatures or pressures, they are capable of modeling with a very high
precision the structural and dynamical properties of a large variety of molecules,
and are therefore vastly used in research and industry.
4.5.2 The Terso potential
The family of potentials developed by Terso [31], are based on the concept
of bond order: the strength of a bond between two atoms is not constant,
but depends on the local environment. This idea is similar to that of the \glue
model" for metals, to use the coordination of an atom as the variable controlling
the energy. In semiconductors, however, the focus is on bonds rather than
atoms: that is where the electronic charge is sitting in covalent bonding.
At rst sight, a Terso potential has the appearance of a pair potential:
V = 21 R (rij ) + 12 Bij A (rij )
X X
(4.15)
ij ij
where R and A mean \repulsive" and \attractive". However, it is not a pair
potential because Bij is not a constant. In fact, it is the bond order for the
bond joining i and j , and it is a decreasing function of a \coordination" Gij
assigned to the bond:
Bij = B (Gij ): (4.16)
Gij is in turn de ned as
fc(rik )g(jik )f (rij ; rik )
X
Gij = (4.17)
k
where fc (r), f (r) and g() are suitable functions. The basic idea is that the
bond ij is weakened by the presence of other bonds ik involving atom i. The
amount of weakening is determined by where these other bonds are placed.
Angular terms appear necessary to construct a realistic model.
This scheme works in a broader spectrum of situations than the Stillinger-
Weber potential, however it is not exempt from problems. One of the biggest
problems is perhaps that the t is dicult: with 6 functions to t and angular
terms, nding a good parametrization is not an easy task.
These potentials have also been applied to hydrocarbons [32].

38
4.6 Long-range forces
So far it has always been assumed that interactions are short ranged, that is,
that the e ect of atom j on the force perceived by atom i vanishes when their
separation rij is larger than a cuto distance Rc . While this is a reasonable
assumption in many cases, there are circumstances where long-range forces are
crucial and must be taken into account. For instance:
 ions are charged, as, say, in a simulation of a salt like NaCl. In this case
there is a Coulomb interaction ( 1=r) between ions
 individual particles are neutral but polarized, such as in water. There will
be dipole-dipole interactions ( 1=r3 )
 van der Waals attractions between atoms ( 1=r6 ), arising from induced
dipole moments and often neglected in simulations, become extremely
important under particular geometric conditions. For instance, the tip
of a scanning tunneling microscope is attracted towards a surface tens of
Angstroms underneath by van der Waals attraction.
In these cases, there is no cuto distance and, more importantly, each parti-
cle interacts with all the images of the other particles in the nearby MD boxes:
the \minimum image rule" introduced in 2.2.1 breaks down. A particle would
also interact with its own images.
A known method to deal with long-range forces is that of Ewald sums.
Brie y, and assuming we are dealing with charged ions, one assumes that each
ion is surrounded by a spherically symmetric charge distribution of opposite sign
which neutralizes the ion. This distribution is localized, that is, it is contained
within a cuto distance. This e ectively screens ion-ion interactions, which are
then treated with conventional short-range techniques.
To restore the original system, one then considers the Coulomb interaction
of similar charge distributions centered on the ions, but now with the same
sign as the original ions. These distributions exactly cancel those that we
have considered before, so that the total interactions that we shall obtain will
be those of the original system made of point charges only. The interaction
of the cancelling distributions is computed in reciprocal space. The required
Fourier transforms are particularly simple when the (arbitrary) shape of the
distributions is chosen to be a gaussian.
The interested reader is referred to ref. [3], x5.5, for more details. An exten-
sion of Finnis-Sinclair potentials for metals to incorporate van der Waals forces
has been made by Sutton and Chen [33].

4.7 But do we really need potentials?

About ten years ago, Car and Parrinello [34] developed a powerful method al-
lowing to perform MD simulations where the ions are still moved classically, but
under the action of forces obtained by solving the electronic structure problem,
thus eliminating the need for empirical potentials at the expense of much larger

39
computational requirements. This important development spawned one of the
most important research lines in today's condensed matter physics, known as
rst-principles (or ab initio) molecular dynamics (FPMD), or simply as Car-
Parrinello method.
Will the existence of the certainly more accurate FPMD technique imply
the death of classical molecular dynamics based on empirical potentials? It is
safe to say that the answer is no. While at present classical MD is being run
on systems containing million of particles on the most powerful computers, or
for simulation times of the order of nanoseconds (10;9 s), the current limit for
FPMD is of the order of a thousand atoms, and its typical simulation times are
measured in picoseconds (10;12 s). This means that several problems requiring
large sizes and/or large times|such as many of those mentioned in x1.4|can
be attacked only by classical methods. Even if the speed of computers keeps
increasing at a breathtaking pace, so does the size of problems of interest, and
it is very likely that all the MD techniques that we know today will remain in
use for quite some time from now, and new ones may appear.
In fact, classical MD and Car-Parrinello MD already represent the extremes
of a spectrum of MD methods, di ering in the degree of approximation used
to solve the electronic problem in order to obtain atomic forces. In particular,
a promising area is that of tight-binding molecular dynamics (TBMD), where
systems containing thousand of particles are now being simulated thanks to the
recent development of O(N )-scaling computer codes.
The reader that wants to be introduced to the fast-growing areas of simu-
lations from rst-principles could start from the review articles present in the
literature [35, 36, 37].

4.8 Fitting to ab initio data by \force matching"

We have seen in x4.5.1 how analytical forms can become pretty complicated
when we strive for realism in modeling. Parametrizing the function appearing in
a complex Hamiltonian can then become a tiresome task, and people often resort
to very simple parametrizations. For instance, Terso potentials have only 12
parameters: an average of 2 parameters to describe each of the 6 functions
constituting the model.
A rst-principles molecular dynamics program can easily produce huge quan-
tities of useful and reliable data, such as forces acting on atoms in a variety of
di erent geometries and coordinations. These data can complement the ex-
perimental quantities normally used to t potentials, helping to construct more
realistic potentials. Recently, a scheme called \force matching method" [38] has
been developed. It consists of a numerical optimization procedure which tries
to match as closely as possible ab initio-derived forces with a classical potential,
using a rich parametrization (10-15 parameters for each function in the poten-
tial). By explicitly including di erent geometries and di erent temperatures in
the data set, one can attack the transferability problem at its very heart. One
can says that the potential is constructed by learning from rst-principles.
Several groups are now working on this front, so that potentials for classes

40
of materials so far unexploited by classical molecular dynamics are likely to
appear in the next years.

41
Appendix A
Choosing hardware and
software
A.1 What kind of computer?
Compared with other applications in today's computational condensed matter
physics, MD simulations using classical potentials are not particularly demand-
ing, unless one's goal is to treat systems of the order of 106 atoms. For a
typical system containing 10000 atoms|a size suitable to tackle a large vari-
ety of problems|and involving short-range potentials, one can estimate a total
memory occupation of the running program between 8 and 16 megabytes (MB).
This could t well in a typical current workstation with a RISC processor1 , in-
cluding entry-level machines.
Using code which scales linearly with the number of particles (as described
later), a RISC CPU can yield a total simulation time of hundred of picoseconds
or perhaps a few nanoseconds in a week of production, which is a reasonable
rate for productive work.
In fact, a popular choice at many sites running di erent applications is to
dedicate the smallest machines to classical MD, leaving the biggest machines
with larger amounts of memory to the more demanding electronic structure
codes.
A.1.1 Selecting a computer for molecular dynamics
Technology advances at an astounding rate. Every few months new, more
powerful machines appear on the market. The market situation is therefore
in continuous evolution, and one should always try to decide what to buy as
closely as possible to the day at which the purchase is actually made.
Trying to select the \best computer for MD" is never easy. The choice is
always conditioned by various factor, such as
1
The name RISC (Reduced Instruction Set Computer) denotes an architecture for com-
puter processor units which ourished in the last decade, consisting of reducing the number
of di erent machine instructions, and optimizing their execution time|possibly at the level
of one clock cycle per executed instruction.

42
  Size of physical problems to study
 Will the machine be used to run other kinds of programs?
 Is fast graphics important? (The site may have another machine speci -
cally dedicated to visualization)
 \Friendliness" of operating systems, compilers and development tools
 Good customer and service support given by the company in your area
 Existence of other similar machines at your site, which would decrease
the system management duties
Obviously, one has to carefully weight these factors keeping into account the
peculiarities of the local situation.
Nevertheless, performance remains perhaps the most important point! To
evaluate performance, it is common to use benchmarks: programs which accom-
plish a certain task (sometimes very simple), are never (or seldom) modi ed,
and are run in controlled conditions. When running a benchmark on a machine
for evaluation, the execution time is measured, and the output is compared
with a reference output to make sure that the machine produced the correct
result.
There are many scienti c benchmarks around. The most well known are
the Linpack2 benchmark, and the SPEC3 program suite. A problem is that
di erent machines behave di erently on di erent programs. In other words, the
ratio between the execution time of two di erent benchmarks can vary widely
between one machine architecture and another. Linpack, for instance, is a linear
algebra benchmark. Machines which are able to perform matrix multiplications
very eciently|perhaps due to a special organization of their oating point
unit or to compilers highly optimized for this problem|will obtain a very good
Linpack score. While extremely representative of an important class of physical
problems (for instance those in which diagonalizations are involved), this is of
limited usefulness to assess MD performance, since there is no linear algebra in
classical MD!
The other extreme, and by far the best choice, is to run yourself a copy of
the program you are going to use on the machines you are interested in. This
will not only give you an idea of a machine's performance, but also of how easy
is using it. However, this is often impractical to do, due to the need to nd
access to these computers, and the time to run the program on each of them
and verify the results.
To partially relieve this problem, the author developed in 1988 a molecular
dynamics benchmark called MDBNCH4 . A database of execution times on sev-
eral machines has been collected over the years, and it keeps growing thanks to
the contribution of many people. Keeping in mind that this benchmark may
not exactly represent the applications which you are interested to run, it is
2
http://netlib.bell-labs.com/
3
http://www.specbench.org/
4
http://www.sissa.it/furio/mdbnch.html

43
probably an indicator of molecular dynamics performance more accurate than
the more general benchmarks available around.
A.1.2 Storage space
On the other hand, the need for disk space should not be underestimated.
Storing positions, velocities and accelerations of 10000 atoms in 64-bit preci-
sion takes 720 KB. In a typical project, a large number of such con gurations
are saved for later reference, and additional statistical quantities are often also
computed, making it easy to ll hundreds of MB of mass storage. Projects
involving dynamical analysis over ne-grained trajectory les, such as calcula-
tions of vibrational spectra, are even more demanding. While the use of tapes
is always possible, availability of a comfortable amount of online storage is an
important ingredient to work smoothly and eciently.

A.2 What language?

A.2.1 High level languages
A computer language is an interface between man and machine.
Machines execute sequences of instructions belonging to their instruction
set. Such instructions are rather simple and tightly bound to features of the
machine hardware. For instance, \load into register A the contents of mem-
ory location 00900752", or \multiply the numbers contained into registers A
and B, leaving the result in register C" could be two examples of instructions
at the machine level. A program written in machine language (or assembler)
would not be easily transferable to a computer produced by a di erent manu-
facturer. Moreover, such a program would be very hard to decipher for a human
examining it.
So called high level computer languages are the solution. A language is a
machine-independent way to specify the sequence of operations necessary to
accomplish a task. A language can be designed to express in a concise way a
\common" sequence of operations. A line in a high level language can embody
powerful operations, and correspond to tens, or hundreds, of instructions at the
machine level. Of course, what is \common" depends on the kind of problem
at hand. For this reason, many di erent languages exist in the computer world,
aimed at solving problems in di erent areas.
The compiler is the tool to convert a program written in a high level language
into the sequence of machine instructions required by a speci c computer to
accomplish the task. Users typically control details of the operation of compilers
by means of options supplied on the command line, or directives embedded in
the program source, but they seldom need to examine the resulting machine
language code.
However, one should not forget the fact that programs are also an interface
with humans, and a way to communicate between humans. A program carries
with it the information of how a problem is being solved, and it is extremely
important for this information to be presented in a clear, complete way to

44
people examining the code. In the scienti c world, quite often the source code
is also the main documentation for the program. For this reason, it is very
important to have in mind people needs|and not only machines needs|when
writing a program. In practice, this means writing informative comments in
the code, using meaningful names for symbols, avoiding \dirty tricks" unless
absolutely necessary. This documenting activity should always be done at the
time of code writing, not deferred to a later time that may never come.
A.2.2 Fortran or C?
Today, the two more popular languages for numeric computation are Fortran
and C.
Fortran has a long tradition in computational physics, chemistry and en-
gineering. Born at the end of the 50s, during the 60s and the 70s it became
the most popular choice for the development of numerically intensive code. An
enormous body of knowledge has been embedded in Fortran codes.
C is a more recent language. It was developed in the 70s and provides
the foundation for the Unix operating system. The popularity of Unix is a
direct consequence of the fact that it is largely written in C, and therefore
easily portable between di erent brands of computers. C is a simple still very
powerful language, and its scope is rather broad: the solution to problems in
very di erent areas has been coded in C with success.
Languages evolve much more slowly than computer architectures. The main
reason for that is that we can easily upgrade to new computers: what we have
to do is to recompile our programs with the compiler supplied with the new
machine. Changing language is a much more expensive and painful business:
we must rewrite our codes, and make sure that the new version does what it
is supposed to do. However, advances in computer systems often create a need
for new features in a language.
For this reason, even a stable and well-proven language such as Fortran
undergoes periodic revisions in order to keep it up to date. Fortran was revised
in 1966 (Fortran 66), then 1977 (Fortran 77), and more recently in 1990 (Fortran
90). At present, there is an immense amount of scienti c codes written in
Fortran 77. Many of them are of public domain, and among them there are
libraries of subroutines for solving a large variety of problems.
Both Fortran 90 and C are perfectly valid choices for a computationally
intensive application. Fortran 77 lacks many features which are crucial today,
and should be regarded as an obsolescent language.
In favor of C:
 ecient compilers are universally available for all the main architectures
in use, and a good public compiler also exists (gcc). C compilers often
come free with machines, while Fortran 90 compilers must be purchased,
and are often expensive
 C is very broad in scope, and is more powerful of Fortran 90 in some
areas, such as pointers and manipulation of strings of characters

45
  acquired coding experience can be directly used outside the scienti c
world: C is very common in the commercial world.
In favor of Fortran 90:
 language designed with easiness of coding numerical problems as one of
the main goals (example: array operations)
 language designed with high performance in numerical applications as one
of the main goals
 interfacing with previously developed Fortran code or directly reusing
portions of existing Fortran code is easier
 No need to master simultaneously two languages when developing new
codes while maintaining old Fortran 77 programs at the same time
 better portability between di erent architectures (for instance, between
32- and 64-bit machines)
The nal choice can be in favor of either one of them, depending to the weight
assigned to the various factors.
For this writer, the advantages of Fortran 90 outweigh those of C, and the
code examples in this notes are written in Fortran 90. Some will disagree with
this choice.

46
Bibliography
[1] D. W. Heermann, Computer simulation methods, Springer, 1986.
[2] Molecular dynamics simulation of statistical-mechanical systems, G. Cic-
cotti and W. G. Hoover (Eds.), North-Holland, 1986. Contains chapters
written by experts of the eld on many aspects of simulation. A place to
dig in into speci c topics, with plenty of technical details.
[3] M. P. Allen and D. J. Tildesley, Computer simulation of liquids, Oxford,
1987. This is a classic how to book, containing a large number of technical
details on a vast array of techniques.
[4] R. W. Hockney and J. W. Eastwood, Computer simulation using particles,
IOP Publishing, 1988. This book mostly addresses simulation methods
where particles are immersed in an external eld as well as with each
other. Such methods are commonly used for instance in plasma physics. A
section on MD simulations on liquids and solids is however present, and
the discussion of the cell method in section 8-4 is particularly useful.
[5] W. G. Hoover, Computational statistical mechanics, Elsevier, 1991. This
book focuses on the foundations of simulation from the statistical mechan-
ics point of view, and contains interesting chapters on ensembles, hydro-
dynamics, non-equilibrium molecular dynamics, with several examples.
[6] J. M. Haile, Molecular dynamics simulation, Wiley, 1992. An excellent and
recent general introduction to molecular dynamics. Not as thorough as
[3] in explaining techniques, but very deep in explaining the fundamentals,
including \philosophical" aspects. Probably the best starting point. A good
amount of Fortran 77 code is included.
[7] D. C. Rapaport, The art of molecular dynamics simulation, Cambridge
Univ. Press, 1995. Excellent recent \recipe book", describing in detail a
large number of techniques, and containing a large amount of C code. The
book has a home page5 .
[8] Simulation of liquids and solids, G. Ciccotti, D. Frenkel and I. R. McDonald
(Eds.), North-Holland, 1987. A collection of key papers on simulation from
the origins to 1986.
5
http://www.cup.cam.ac.uk/onlinepubs/ArtMolecular/ArtMoleculartop.html

47
 [9] W. W. Wood, in [2], p. 3.
[10] B. J. Alder and T. E. Wainwright, J. Chem. Phys. 27, 1208 (1957).
[11] J. B. Gibson, A. N. Goland, M. Milgram, and G. H. Vineyard, Phys. Rev.
120, 1229 (1960).
[12] A. Rahman, Phys. Rev. 136, A405 (1964).
[13] L. Verlet, Phys. Rev. 159, 98 (1967); Phys. Rev. 165, 201 (1967).
[14] J.-P. Hansen and L. Verlet, Phys. Rev. 184, 151 (1969).
[15] S. Yip, in [2], p. 523.
[16] H. J. C. Berendsen, in [2], p. 496; Comp. Phys. Commun. 44, 233 (1987).
[17] J. P. Hansen and I. R. McDonald, Theory of simple liquids, 2nd Ed., Aca-
demic, 1986. A classic book on liquids, with a chapter devoted to computer
simulation.
[18] See for instance C. Kittel, Introduction to solid state physics, Wiley.
[19] H. J. C. Berendsen and W. F. van Gunsteren, in [2], p. 43.
[20] Interestingly, one can nd a description of the leap-frog algorithm in R. P.
Feynman, R. B. Leighton and M. Sands, The Feynman Lectures on Physics,
Vol. 1, Addison-Wesley, 1963, Chapter 9 (\Newton's Laws of Dynamics").
[21] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, Science 220, 671 (1983).
[22] H. C. Andersen, J. Chem. Phys. 72, 2384 (1980).
[23] M. Parrinello and A. Rahman, Phys. Rev. Lett. 45, 1196 (1980); J. Appl.
Phys. 52, 7158 (1981); M. Parrinello, in [2], p. 204.
[24] S. Nose, Molec. Phys. 52, 255 (1984); W. G. Hoover, Phys. Rev. A 31,
1695 (1985).
[25] F. Ercolessi, M. Parrinello, and E. Tosatti, Philos. Mag. A 58, 213 (1988).
[26] A. Carlsson, Solid State Phys. 43, 1 (1990).
[27] M. S. Daw and M. I. Baskes, Phys. Rev. B 29, 6443 (1984); S. M. Foiles,
M. I. Baskes, and M. S. Daw, Phys. Rev. B 33, 1986.
[28] M. W. Finnis and J. E. Sinclair, Philos. Mag. A 50, 45 (1984).
[29] K. W. Jacobsen, J. K. Nrskov and M. J. Puska, Phys. Rev. B 35, 7423
(1987).
[30] F. Stillinger and T. A. Weber, Phys. Rev. B 31, 5262 (1985).
[31] J. Terso , Phys. Rev. B 37, 6991 (1988).

48
[32] D. W. Brenner, Phys. Rev. B 42, 9458 (1990).
[33] A. P. Sutton and J. Chen, Philos. Mag. Lett. 61, 139 (1990).
[34] R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985).
[35] D. K. Remler and P. A. Madden, Molec. Phys. 70, 921 (1990).
[36] M. C. Payne, M. P. Teter, C. C. Allan, T. A. Arias, and J. D. Joannopoulos,
Rev. Mod. Phys. 64, 1045 (1992).
[37] G. Galli and A. Pasquarello, in Computer simulation in chemical physics,
M. P. Allen and D. J. Tildesley (Eds.), Kluwer, 1993, p. 261.
[38] F. Ercolessi and J. B. Adams, Europhys. Lett. 26, 583 (1994).

49