M

ETODOS COMPUTACIONALES EN
FSICA DE LA MATERIA CONDENSADA
(COMPUTATIONAL METHODS IN
CONDENSEDMATTER PHYSICS)
http://www.uam.es/enrique.velasco/master/cmcmp
Master en Fsica de la Materia Condensada
y Nanotecnologa
PART I
E. Velasco
Departamento de Fsica Te orica de la Materia
Condensada
Facultad de Ciencias
Universidad Aut onoma de Madrid
1
CONTENTS
Bibliography 3
Tools used in this course 3
A brief introduction to computation in Condensed Matter Physics 4
The Monte Carlo method 6
The method of Molecular Dynamics 54
List of computer codes 69
2
BIBLIOGRAPHY
Most of the topics discussed in these notes are reviewed in more detail in one of the
following textbooks:
D. Frenkel and B. Smit, Understanding molecular simulation: from algorithms to
applications (Academic Press Inc., San Diego, 1996).
M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Oxford University,
1987)
The following research paper may also be useful:
L. Verlet, Computer Experiments on Classical Fluids. I. Thermodynamical Prop-
erties of Lennard-Jones Molecules, Phys. Rev. 159, 98 (1967).
Francis H. Ree and William G. Hoover, Fifth and Sixth Virial Coecients for Hard
Spheres and Hard Disks, J. Chem. Phys. 40, 939 (1964).
Copies of these papers can be downloaded from the course web page:
http://www.uam.es/enrique.velasco/master/cmcmp
In this web page a more complete bibliography can be found.
TOOLS USED IN THIS COURSE
These notes contain a number of computer codes written in the FORTRAN language, which
can be downloaded from the course web page. Therefore access to a FORTRAN compiler
is required. Also, use of a graphics package, such as gnuplot, will be necessary. The
practical problems of the course, to be done in the PC lab, will require use of these tools.
3
1. A BRIEF INTRODUCTION TO COMPUTATION IN CONDENSED
MATTER PHYSICS
Computational methods are nowadays essential to explore the behaviour of condensed
phases of matter. By computational methods we do not mean methods that simply rely
heavily on the use of the computer to solve a particular problem; for instance, a meso-
scopic Hamiltonian of the GinzburgLandau type, in terms of elds representing local
coarsegrained order parameters, could be formulated to obtain the structure of a com-
plex uid limited in space by some complicated boundaries. For sure, this problem would
demand heavy computational resources since the functional representing the Hamiltonian
would have to be minimised with respect to a huge number of variables.
Rather, the methods we are set ourselves to describe (albeit in a sketchy way) are methods
that are formulated in terms of microscopic variables dening the state of assemblages of
particles, either quantum or classical, and interactions between these particles. Particle
variables appear explicitely in the problem, which is then microscopic in nature, and we
would like to follow the evolution of these particles in time, and/or simply obtain average
properties of the system in order to make contact with the macroscopic properties. In
essence, we are looking here at statistical mechanics from a purely microscopic point of
view, using the basic methods of classical and quantum mechanics, and discuss methods
to solve the microscopic equations describing the microscopic variables. These methods
can be properly called computer simulation.
In the last 60 years computers have brought about a complete upheaval of the eld.
Today it is possible to follow the dynamical evolution of systems of millions of particles
interacting in more or less complicated ways on a computer. This is still far from the
thermodynamic limit, or may even be far from typical numbers in surface science, but
nite-sized systems of this kind already contain many of the essential properties of real
systems. The aims and methods of statistical mechanics can now be explicitely put to
work on a computer. The two basic methods of computational statistical mechanics, the
Monte Carlo (MC) method and the Molecular Dynamics (MD) method will be introduced
in the rst part of this course. The two basic approaches of statistical mechanics, based
on time averages and ensemble averages, are played by the MD and MC methods, respec-
tively.
The role played by computer simulation is two-fold. Given a physical system, one de-
vises a model to describe the essential features of the system. Now this model can be
analysed using theoretical approaches or computer simulation. Since the rst necessarily
entail approximations, whereas the second are essentially exact, a comparison of these two
allows to draw conclusions as to the accuracy of the theoretical approximations involved
in the theoretical approach. Experiments, on the other hand, are made directly on the
real system, and a comparison with simulation results can be made to obtain information
on the adequacy of the model used to describe the real system. Computer simulation thus
plays a twofold role that has helped improve our understanding of matter very signi-
cantly. Computer simulation does not lack limitations and problems, and has become a
specialised eld which requires good training, in much the same way as an experimenter
or theoretician are experts in their own elds. The division between simulator and the-
4
oretician or experimenter, however, is not clearcut.
Limitations of computer simulations
There are spatial limitations in computer simulation. Systems of interest are usually
very large (macroscopic) in size, and properties may sometimes be sensitive to this. It
is then necessary to apply nite-size scaling theory to extract the thermodynamic limit.
Also, close to secondorder phase transitions or in uid interfaces, large uctuations must
be treated with care and may require larger system sizes than amenable to computational
resources available.
Also, there are time limitations. For example, in aggregation dynamics, phenomena such
as critical slowing down near secondorder phase transitions, or dynamics in inhomo-
geneous systems, accessible computational times may sometimes be too short. When
dealing with time scales close to typical hydrodynamical time scales (seconds), i.e. when
there is macroscopic ow, long simulations are required. Relaxation times associated to
covalent bonds, molecular conformations or translations are usually dealt with by present
computers without too much eort.
5
2. THE MONTE CARLO METHOD
In physics and mathematics, any method that makes use of random numbers is called
Monte Carlo method. We will start with a short introduction to dierent methods to
generate random deviates, or random numbers, in computational physics, and then will
illustrate the importance-sampling technique in thermal physics using the Ising model.
2.1. How to generate uniformly distributed random numbers
The general problem that we set ourselves to solve is how to generate random devi-
ates x, distributed according to some particular probability function f(x) in the interval
[a, b]. Remember that
f(x)dx = probability that x lies within the interval [x, x + dx]
All FORTRAN compilers come with intrinsic functions that generate uniformly distributed
random numbers, but very often we lack information as to which methodology is used, or
what the performance of these generators is. Sometimes these routines are not suciently
accurate for computations in statistical physics.
Virtually all methods that are available to solve this problem rely on methods to generate
uniformly distributed random numbers, i.e. numbers distributed according to a constant
probability function in the interval [0, 1]. This is the problem on which we now focus our
attention.
f(x)
x
f(x)
x
Figure 1: Left: uniform distribution function. Right: histogram obtained from a nite sequence
of uniformly distributed random numbers (i.e. numbers distributed according to the function
on the left).
First we illustrate the method of congruences, using integer arithmetics. We rst note
that the digits in the product of these two integer numbers,
12345 65539 = 809078955
look random. If we truncate the last ve digits, we get
08090[78955 78955 65539 = 5174631745
51746[31745 31745 65539 = 2080535555
20805[35555 35555 65539 = ...
6
The number 12345 is called the seed, while the number 65539, which is used as a constant
factor to multiply in all steps is called the base. If we divide the results of these successive
truncation operations by 10
5
, we obtain the following sequence of numbers
0.78955, 0.31745, 0.35555, ...
which are apparently random; actually, since we have obtained the numbers through a
well-dened sequence of steps (an algorithm), they are called pseudo-random numbers. A
histogram of a nite sequence of these numbers would give a plot similar to that shown
on the right panel in Fig. 1. On taking N , the histogram would tend to atten out
to unity in the whole interval.
0 2147483647 -2147483648
Figure 2: Range of integer numbers that can be stored in a register of a 32bit computer.
This procedure may pose a problem when implemented on some computers. For example,
in a 32-bit computer (4 bytes; each byte contains 8 bits) can only store in their internal
CPU registers signed numbers between 2
31
and 2
31
1 (see Fig. 2). Any number larger
than 2
31
1 is truncated, and can even be turned negative! For example, the product
78955 65539 = 5174631745
is represented, in binary format, as
5174631745
10
= 100110100011011101001110101000001
2
which is 33bit number. On a 32bit computer, the 32th bit (counting from the right),
which is a 0 in our case, indicates negative sign (a 1 would be for a positive number),
but the 33th bit (the rst from the left), which is a 1, is lost forever. The number would
be stored in the register as
00110100011011101001110101000001
2
which, in decimal base, is 879664449 (far from the correct result, 5174631745). An even
more dramatic example is
51326 65539 = 931112582
In order to avoid this problem, if the number is negative we add the period, 2147483648 =
2147483647 + 1 (it is necessary to use the second form, 2147483647 + 1, since again
2
31
= 2147483648 would not be accepted by the computer).
Things are a little easier in real arithmetics. In the following a possible FORTRAN code,
based on the above method (but with a base of 16807 and using real numbers), is pro-
vided. These lines of code generate a random number, uniformly distributed: note that
7
the singleprecision variable r is a random number in [0, 1]. The variable dseed is a
doubleprecision number (e.g. dseed=173211d0) provided by the user.
001. !
002. ! P1: Random number generator
003. !
004. subroutine ggub(dseed,r)
005. real*8 z,d2p31m,d2pn31,dseed
006. data d2p31m/2147483647./,d2pn31/2147483648./
007. z = dseed
008. z = dmod(16807.*z,d2p31m)
009. r = z / d2pn31
010. dseed = z
011. end
The crucial lines are 008 and 009. In line 008, the function dmod multiplies by the base,
divides by the number stored in d2pn31 (which plays the role of the 10
5
of the previous
example) and takes the remainder. Then division by d2pn31 puts the number in the
interval (0, 1). The values of the parameters are optimised to give a stable and reliable
random number generator.
Once one knows how to generate uniformly distributed numbers, i.e. numbers distributed
according to the probability function f(x) = 1 in the interval [0, 1], it becomes an easy
task to generate numbers distributed uniformly but according to the constant (normalised)
probability function f(x) = 1/(b a) in [a, b]: we just have to generate y [0, 1] and
make x = a + (b a)y.
2.2. How to generate random numbers distributed according to a general
f(x)
There are dierent methods to obtain random numbers that are distributed following
a general probability function f(x). We will review three methods.
Acceptancerejection method. This is a not very ecient method, since it
involves generating two random numbers to obtain a single number distributed
according to f(x). We limit ourselves to stating the procedure involved, without
mentioning the basis of the algorithm. The method is based on the fact that if
we generate two sequences of uniform random numbers, x and y, then the set of
numbers
x/f(x) y
is distributed according to f(x). The steps are:
1. we generate uniform numbers x and y in [a, b]
2. if f(x) y, we accept x; otherwise x is rejected
3. we go to step 1
8
Inversion method. This method uses the cumulative function:
P(x) =

x
a
dx

f(x

), P(a) = 0, P(b) = 1
Since P(x) is a monotonically increasing function, it can be inverted to give x =
P
1
(y). Then, if y is uniformly distributed in the interval [0, 1], then x is distributed
according to the probability function f(x) in the interval [a, b]. The feasibility of the
method relies on the possibility of explicitely obtaining the inverse function of P(x)
(we know it can be done, but it may be necessary to do it numerically, which would
slow down the process considerably and be a disadvantage from a computational
point of view).
Method of Metropolis et al. This is a very general and elegant method, although
it has some disadvantages also. The method is of such importance in statistical and
condensedmatter physics that it deserves to be given more time than the others.
Again, the theoretical basis of the method, the theory of Markov chains, will not be
mentioned at all, and we will focus on the algorithm and on a few tricks that have
to be followed in order to implement eciently the algorithm.
The method involves generating a sequence of numbers,
x
[0]
, x
[1]
, x
[2]
, ...
that, asymptotically, are distributed according to the probability function f(x). To
obtain this, we begin with a seed, x
[0]
, and obtain from it a test number, which
is accepted or not as a new number in the sequence with some probability. The
process is repeated as many times as necessary. In each step n, the value x
[n+1]
is
then obtained from x
[n]
, using the following algorithm:
we choose a test number as
x
[t]
= x
[n]
+ x
where [1, 1] is a uniform random number, and x a constant number.
This value is accepted as a new member of the sequence of random numbers
with probability
r =
f(x
[t]
)
f(x
[n]
)
which means: we get a random number [0, 1], so that

si r > the test value is accepted, i.e. x

[n+1]
= x
[t]
si r < the test value is not accpeted, i.e. x
[n+1]
= x
[n]
A few important points to bear in mind about the algorithm:
The method generates numbers distributed according to f(x), but only asymp-
totically. Therefore, we need a previous warmup period, long enough to reach
the asymptotic regime as much as possible, and it is after this initial period
that the random numbers can be taken as acceptable
9
The value of x is usually adjusted so that approximately 50% of the tested
numbers are accepted. This insures a quick convergence to the asymptotic
regime (note that rejections are as equally needed as acceptances!)
The workings of the method only rely on the ratio r = f(x
[t]
)/f(x
[n]
). In
particular, r does not depend on the normalisation of the distribution function
(which is obviously constant). It is this feature that makes the Metropolis et al.
algorithm so useful in applications to statistical and condensedmatter physics.
The Metropolis et al. algorithm may be easily extended to more than one
random variable, x = (x
1
, x
2
, ...). In this case we obtain a sequence of vectors
x
[0]
, x
[1]
, ..., and generate test vectors as x
[t]
= x
[n]
+x, following exactly the
same steps as in the case of a single variable
As an example, we apply the method to obtain a sequence of random numbers with
distribution function f(x) exp(x
4
) in the innite interval (, ). This is
clearly a dicult or utterly impossible task for the methods mentioned above: the
inversion method does not apply since the cumulative function cannot be explicitely
calculated, let alone inverted; and the acceptancerejection method is inecient).
Consider the function
f(x) = ce
x
4
,
where c is a constant such that

dxf(x) = 1.
The value of c does not need to be known to implement the Metropolis et al.
algorithm; however, we will calculate c to compare the results with the normalised
distribution function. The normalisation integral can be expressed in terms of the
function by making the substitution x
4
= t:

dxe
x
4
= 2

0
dxe
x
4
=
1
2

0
dtt
3/4
e
t
=
1
2

1
4

= 1.8128...
A FORTRAN code that implements the algorithm for this problem follows.
001. !
002. ! P2: Metropolis et al. algorithm
003. ! for f(x)=c exp(-x**4)
004. !
005. parameter (n=1000)
006. real*8 dseed,his(-n:n),h,x,rr,xt,amp,darg,anorma
007. real*4 r
008. integer i,m,j,k,n
009. h=0.1d0
010. dseed=53211d0
011. amp=0.2d0
012. m=10**6
013. nacc=0
10
014. do i=1, n
015. his(i)=0
016. end do
017. x=0.2
018. do j=1, m
019. call ggub(dseed,r)
020. xt=x+(2*r-1)*amp
021. darg=xt**4-x**4
022. rr=dexp(-darg)
023. call ggub(dseed,r)
024. if(dble(r).lt.rr) then
025. x=xt
026. nacc=nacc+1
027. end if
028. k=dnint(x/h)+1
029. if(k.gt.n.or.k.lt.-n) stop k fuera de intervalo
030. his(k)=his(k)+1
031. end do
032. anorma=0d0
033. do k=-n, n
034. anorma=anorma+his(k)
035. end do
036. do k=-n, n
037. write(*,*) (k-1)*h,his(k)/(anorma*h)
038. end do
039. end
The histogram is stored in vector his. The step interval is h, set to 0.1 in line
009, and the amplitude x is set to 0.2 in line 011 (in the variable amp. The num-
ber of Monte Carlo steps is set to 10
6
in line 012, and in lines 014-016 the histogram
is reset (since it is an accumulator, it is convenient to set it to zero at the beginning).
Monte Carlo steps are done in lines 018031. First, at test value xt (line 020) is
chosen, and the probability ratio r (here rr) is calculated in line 022. The accep-
tance criterion is applied in line 024, and lines 025026 deal with the case where
acceptance results (updating the value of the random variable and adding unity to
the number of accepted moves (nacc). Note that, after line 027, x contains the new
value, if it was accepted, or the old one, and that in either case the value of x is
accumulated in the histogram (in line 030). After all the m Monte Carlo steps have
been done, lines 032035 calculate the norm of the histogram (in order to normalise
it), and in lines 036038 results are dumped to the screen.
In Fig. 3 the resulting histogram is shown. The value of x (amp in the code)
can be adjusted so that the acceptance percentage oscillates about 50% (this can be
done in the warmup period, by computing the acceptance rate along the simulation
run and changing x according to how this ratio varies: if less than 50%, x is
decreased, thus increasing the acceptance ratio, and the other way round); this is
not implemented in the present code. In the gure calculations are presented for
11
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
-2 -1 0 1 2
x
x
x
(a) N=10
3
(b) N=10
4
(c) N=10
7
f
(
x
)
f
(
x
)
f
(
x
)
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
-2 -1 0 1 2
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
-2 -1 0 1 2
Figure 3: Histogram for the distribution function f(x) = ce
x
4
, x (, ), as obtained from
the Metropolis et al. algorithm. (a) N = 10
3
, (b) N = 10
4
and (c) N = 10
7
random numbers.
three values of the number of steps, N = 10
3
, 10
4
and 10
7
. In the rst case the
distribution is still far from the real distribution; we would still be in the warm
up period. At least an order of magnitude more steps are necessary to enter the
asymptotic regime. Then we can conclude that at least 10
4
steps are necessary in
the warmup period; these initial steps are not to be included in the nal period,
where we can use the sequence of random numbers to do whatever we intend to do
with them (either construct a histogram, as in the present case, or evaluate an inte-
gral, or something else). Of course, if we do not know the real distribution function
(which is normally the case) other criteria must be used to ascertain whether the
asymptotic regime has been reached or not.
2.3. Applications in statistical physics
Problems in statistical and condensedmatter physics are formulated in terms of mod-
els that generally involve a huge number of independent variables. There are two basic
problems, thermal and nonthermal, which normally require dierent treatments. Non
thermal problems are generally easier to deal with, since the Boltzmann distribution
associated with a temperature is not involved. We will begin with a few examples of
12
nonthermal problems, such as the random walk model in a number of versions and the
percolation model, with a view to showing some typical techniques involving the genera-
tion of uniform random numbers.
2.3.a. Random walk (RW)
This model is used to describe a number of problems in statistical physics, for exam-
ple, it can be used to characterise Brownian motion, ot the global properties of a polymer
macromolecule in solution in the limit of high exibility. Although the model can be
formulated in many dierent ways, we will use a discrete two-dimensional square lattice.
0
3
4
2
1
Figure 4: Central point (origin), marked with 0, from which a random walk is generated as a
sequence of N steps. In the rst step a random direction (either 1, 2, 3 or 4) is chosen; the new
point is used to again generate another move in the same way. After N steps we get a chain of
N segments: a random walk of N steps.
Consider the lattice represented in Fig. 4. We start from an origin, set at some particular
node of the lattice and labelled with 0 in the gure. We generate a random walk (RW)
of N steps by rst dening the four unit vectors
v
1
= (+1, 0), v
2
= (0, +1), v
3
= (1, 0), v
4
= (0, 1)
and then using the following algorithm:
1. set r
0
= (0, 0) and n = 0
2. choose an integer random number m
n
from the set 1, 2, 3, 4
3. make r
n+1
= r
n
+v
mn
4. if n = N, set R = r
n
; otherwise we go back to step 2
The vector R connects the origin with the last point generated in the random walk.
Clearly, if we generate two dierent random walks, the vector R will be dierent. Fig 5
shows three dierent realisations of a random walk with N = 500 steps.
13
For a lattice with a coordination (number of neighbours) z the number of dierent Nstep
random walks is
Z
RW
N
= z
N
If the random walk were a representation of a polymer chain, Z
RW
N
would be the partition
function of the polymer (assuming the congurations of the polymer had zero energy).
-30
-20
-10
0
10
20
30
-30 -20 -10 0 10 20 30
-30
-20
-10
0
10
20
30
-30 -20 -10 0 10 20 30
-30
-20
-10
0
10
20
30
-30 -20 -10 0 10 20 30
Figure 5: Three realisations of a random walk with N = 500 steps.
Let us now look at the statistical behaviour of the distance between the origin and the
last point generated, i.e. the end-to-end distance of the random walk [R[. We would like
to calculate the average '[R[
2
` over a large enough number of RW chains. We have:

R
2

[R[
2

n=1
v
mn

n=1
[v
mn
[
2

n=1
N

p = 1
n = p
v
mn
v
mp

By denition [v
mn
[
2
= 1 and, since moves are independent and therefore uncorrelated,

n=1
N

p = 1
n = p
v
mn
v
mp

= 0
We then have

R
2

RW
= N, r
N

'R
2
`
RW
=

N
where r
N
is a way to dene an end-to-end distance. This result holds in the limit N .
The following code generates a number of random walks of dierent lengths.
001. !
002. ! P3: Random walk
003. !
14
004. implicit real*8(a-h,o-z)
005. real*8 dseed
006. real*4 r
007. dimension i(4),j(4)
008.
009. data i/1,0,-1,0/
010. data j/0,1,0,-1/
011.
012. dseed=51323d0
013. Ntot=20
014. M=10**6
015. ar2=0
016.
017. do l=1,M
018. x=0
019. y=0
020. do n=1,Ntot
021. call ran (dseed,r)
022. k=4*r+1
023. x=x+i(k)
024. y=y+j(k)
025. end do
026. r2=x**2+y**2
027. ar2=ar2+r2
028. end do
029. write(*,(N, M, <r2>=,2i8,f12.5)) Ntot,M,ar2/M
030.
031. end
The vectors i(4), j(4) give the xy coordinates of the four possible directions the RW
can take from a given lattice site. The length of the RW chain is stored in variable Ntot.
Variable M is the number of congurations for the chain generated. ar2 is an accumulator
for the end-to-end distance of the RW squared, and is set to zero on line 015. The main
loop, over dierent realisations of the RW, is started on line 017. In lines 020025 a
particular chain is generated. The square of the end-to-end distance of that chain is then
computed on line 021, and on the following line accumulated. Results are then presented
on the screen.
With this code, or with slight modications thereof, the following results were obtained.
Fig. 6(a) represents the value of 'R
2
` /N with respect to the number of chains generated
for a value N = 20. It can be seen that the results tend to unity as M, the number of
chains, increases. In Fig. 6(b) the value of 'R
2
` is plotted with respect to N, the length
of the RW. We can also see how the results tend to conrm the theoretical prediction that
'R
2
` increases linearly with N (the continuous line has slope one).
15
0.6
0.7
0.8
0.9
1.0
1.1
1.2
1.3
1.4
10
0
10
1
10
2
10
3
10
4
10
5
10
6
10
7
10
8
no. RW chains
<
R

>
/
N
2
0
50
100
150
200
250
0 50 100 150 200 250
N
<
R

>
2
(b) (a)
Figure 6: (a) Ratio of

R
2

to RW length with respect to number of chains generated. The

ratio tends to unity as the latter number increases. (b)

R
2

as a function of chain length; the

number of chains generated was 10
6
for each value of N. The continuous line is a straight line
with slope equal to one.
2.3.b. Modied random walks
The random walk is a convenient model when we would like to describe step conduc-
tion in solids, diusion processes in lattices, etc. But in applications to polymer physics
the model presents some shortcomings, namely:
a polymer chain cannot turn back on itself
a polymer chain cannot intersect itself
Both these problems stem from the fact that the random walk model neglects the excluded
volume between the units that make up the polymeric chain. The non-reversal random
walk (NRRW) model corrects for the rst problem, by rst dening
v
n4
= v
n
and then modifying step 2 of the RW for n > 1:
1. make r
0
= (0, 0) and n = 0
2. choose an integer random number m
n
from the set m
n1
1, m
n1
, m
n1
+ 1
3. make r
n+1
= r
n
+v
mn
4. if n = N we set R = r
n
; otherwise go to step 2
The partition function is now
Z
NRRW
N
= (z 1)
N
16
To avoid the second problem, we can introduce the socalled self-avoiding random walk
(SARW) model, which does not generate congurations where the chain can intersect
itself. The following additional condition in step 3 of the NRRW algorithm is added:
3

if the node r
n+1
has been already visited, the process stops and we start it over
again.
The partition function of the SARW model is much more complicated; actually an analytic
expression is known only in the limit N :
Z
SARW
N
N
1
z
N
e
N
, z
e
z 1
where is a critical exponent, and z
e
is the eective coordination number, neither of
which can be calculated exactly except in some particular cases; for instance, in one di-
mension, z
e
= z 1, with z = 2. Once again, this is demonstration of the importance of
numerical methods in problems of condensedmatter physics.
As an illustration of the dierence between the three models, RW, NRRW and SARW,
Fig. 7 gives some examples of possible chains for each model.
1
2
3
4
5
6
7 8 9
10
11
1
2
3
4
5
6
20
7 8 9 14
15
16 11
12
13
10
17 18 19
1
2 3
4 5
6,8 7
9,12
10 11
13 14 15 16 19,20
18 17
RW NRRW SARW
Figure 7: Three possible chains for the RW, NRRW and SARW models. All cases were generated
using a chain length N = 20. The labels indicate the sequence of nodes visited. In the SARW
case the same sequence of random numbers as in the NRRW case was used, but the process
stopped at the 12th step because a cross path would have followed.
When N is large, the SARW algorithm becomes very inecient since most of the attempts
to make the chain grow are rejected. In this limit, the ratio of accepted to attempted
steps can be calculated as:
P
N
=
Z
SARW
N
Z
NRRW
N
N
1

z
e
z 1

N
= e
N log
z1
z
e
+(1) log N
N
Thus, the probability of having a rejected move decreases exponentially. This exponential
ineciency makes the algorithm impractical for N 100 or larger (see Fig. 8(a) for
17
actual results from Monte Carlo simulation).
There are no analytical expressions for the NRRW and SARW models, due to the fact
that moves are highly correlated and the mathematics of the problem becomes intractable.
However, it is very easy to perform computer simulations. For instance, from Monte Carlo
simulations it is known that

R
2

SARW
= N
2
, 0.59 (1)
where is another critical exponent. Fig. 8(b) show results from Monte Carlo simulations
using the code below; a power law with the correct exponent is obtained.
We now present a possible FORTRAN code for the SARW algorithm.
001. !
002. ! P4: Algorithm for the SARW model
003. !
004. implicit real*8(a-h,o-z)
005. real*4 r
006. dimension i(4),j(4)
007. dimension ic(-200:200,-200:200)
008. dimension xx(0:200),yy(0:200)
009. data i/1,0,-1,0/
010. data j/0,1,0,-1/
011.
012. dseed=263d0
013. Nchains=2000
014. do Ntot=5,60,5
015. mm=0
016. ar2=0
017. 3 n=0
018. x=0
019. y=0
020. ix=0
021. iy=0
022. do i1=-Ntot,Ntot
023. do j1=-Ntot,Ntot
024. ic(i1,j1)=0
025. end do
026. end do
027. xx(n)=0
028. yy(n)=0
029. 1 call ggub(dseed,r)
030. if(n.gt.0) then
031. kk=3*r-1
032. k1=k+kk
033. if(k1.gt.4) k1=k1-4
034. if(k1.lt.1) k1=k1+4
18
035. k=k1
036. else
037. k=4*r+1
038. end if
039. ix1=ix+i(k)
040. iy1=iy+j(k)
041. if(ic(ix1,iy1).eq.1) then
042. go to 3
043. end if
044. ic(ix1,iy1)=1
045. ix=ix1
046. iy=iy1
047. x=x+i(k)
048. y=y+j(k)
049. n=n+1
050. xx(n)=x
051. yy(n)=y
052. if(n.gt.Ntot) then
053. ar2=ar2+x**2+y**2
054. mm=mm+1
055. if(mm.gt.Nchains) then
056. ar2=ar2/Nchains
057. write(*,*) dlog(dble(Ntot)),dlog(ar2)
058. go to 5
059. end if
060. go to 3
061. end if
062. go to 1
063. 5 end do
064.
065. end
The strategy is to set up a square lattice with NtotNtot nodes, and introduce an occu-
pation matrix called ic, dened on the lattice nodes, which equals one if a given has been
visited, and equals zero if the node has not. In this particular code the number of chains
generated was Nchains=2000, and the chain lengths Ntot varied from 5 to 60 (line 014
opens a do loop over chain lengths). mm is a variable that contains the number of chains
generated for some particular length Ntot, while ar2 accumulates the end-to-end distance
of a successfully generated chain (this is set to zero on line 016). Vectors xx, yy contain
the lattice coordinates of all the points that belong to the chain (this is actually not used
in the present code, but stored for other purposes in expanded versions of the code). On
lines 022026 all the occupation variables are set to zero. The strategy is to start from
the origin (lines 020021, where ix, iy are the position of the current node visited, and
lines 027028) and move so as to avoid visiting the previous point (this is done using the
algorithm on lines 030038) and then check whether the node has been visited already
(lines 041043). If not, the occupation of the node is set to unity, and the coordinates
of the node stored. n contains the number of nodes contained in the given chain that is
19
being generated (line 049). Then it is checked whether the chain has already grown to
the required length (line 052), and if so its end-to-end distance is accumulated. Line 052
checks if the number of chains of length Ntot has been already reached. If this is the case
the average is computed and shown on the screen (lines 056057).
0
20
40
60
80
100
0 10 20 30 40 50
N
%

a
t
t
e
m
p
t
s

r
e
j
e
c
t
e
d
1
10
100
1000
1 10 100
N
N
1.18
<
R


>
2
(a) (b)
Figure 8: Monte Carlo generation of SARW chains. (a) Percentage of rejected attempts to
build a chain as a function of chain length. The continuous line is a t to a Gaussian function.
(b) Average end-to-end distance squared for chains as a function of chain length; each point is
an average over 2 10
3
chains. The theoretical power law N
2
, with 2 = 1.18, is represented
as a continuous line.
Results are shown in Figs. 8(a)(b). The percentage of rejected attempts to build a chain
as a function of chain length N is shown in Fig. 8(a). We can see that as the chains
get longer, the eciency of the algorithm decreases quite dramatically (actually, the data
points in the gure have been tted to a Gaussian function, although for large N the
function decreases exponentially). In Fig. 8(v) the value of 'R
2
` is represented with re-
spect to the number of steps N in logarithmic scale, to better check that the exponential
law, with the correct exponent 2 1.18, is obtained (we have to bear in mind, however,
that the chains generated for each value of N are only 2000 in number, so that averages
over chains are short, while the power law in Eqn. (1) is only correct in the limit N ).
2.3.c The percolation problem
Percolation is a simple geometrical problem, but the solution is far from trivial. The
simplest form of percolation is the following: consider a square lattice with N sites, each
of which can be lled with probability p and left empty with probability 1 p. We dene
a cluster of size l as a group of l connected nearest neighbours. Fig. 9 shows an example.
Let n
l
(p) be the number of clusters of size l when the occupancy probability of each site
is p, and P(p) the fraction of occupied sites (number of occupied sites with respect to the
20
total number of sites N) that belong to the largest cluster.
Figure 9: Example of a possible conguration in a percolation problem. The circles represent
occupied sites, while the bonds between nearest neighbours have been indicated. In this case
there is a percolating or spanning cluster that goes throughout the lattice, along with cluster of
dierent sizes.
The single most important problem associated with percolation is the existence of a per-
colation transition with respect to p, which occurs when there appears a cluster in the
system that percolates, i.e. spans the whole system. The transition takes place when
p = p
c
, where p
c
is a critical probability. A spanning or percolating cluster is a cluster
with a size similar to that of the system, and therefore one can go from one side of the
system to the other, in at least one direction, by just jumping from one neighbour to the
other without leaving the percolating cluster. In Fig. 10 the fraction of occupied sites
that belong to the largest cluster, P(p), is shown as a function of p, as results from a
simulation on a (two-dimensional) square lattice. If the system is nite (N < ) P(p) is
a monotonically increasing function of p, but exhibits an abrupt change about p
c
, which
can be associated with the percolation transition of the innite system, N . In the
innite system, by contrast, P(p) increases abruptly from zero with innite derivative at
p = p
c
, and then slowly increases up to unity (in the latter case the percolating cluster
spans the whole system).
There is a simple relation between P(p) and n
l
(p):

l
(nite clusters)
l
n
l
(p)
N
+ p P(p) = p
21
The percolation problems can be analysed easily by means of the Monte Carlo method:
we ll each site of the lattice with probability p. To do that, we visit the sites sequentially
(one after the other), and for each site we generate a uniform random number in [0, 1]. If
< p, we occupy the site with an atom, otherwise we leave it empty. After visiting the N
sites we will have generated one conguration. By repeating the process we can generate as
many congurations as we wish, say M congurations, all of them resulting from the same
value of p. The complicated problem is how to identify the clusters of each conguration.
There are a number of dierent algorithms to perform this operation. Assuming we have
been able to do this, the nal step is to average over all the M congurations, and we will
obtain P(p). It is known from Monte Carlo simulation that, for the (twodimensional)
square lattice, p
c
0.59.
0
1
0 p 1
c
p
P
(
p
)
Figure 10: Qualitative behaviour of P(p) for an innite system (continuous line) and a nite
system (dotted line).
The percolation problem is rather similar to a secondorder transition. The following
identications are necessary:
P(p) order parameter
G(p) =

l
n
l
(p) free energy
Then the problem has critical and scaling properties. For instance, there appear singu-
larities at p = p
c
, which take the form of power laws:
P(p) (p p
c
)

, G(p) (p p
c
)
2
,
where and are critical exponents. It is clear that any quantitative analysis of the
percolation problem rests on the use of the computer. The problem is interesting not just
from a theoretical point of view, but also in practical applications; for example, in diluted
magnetic alloys, in conduction problems in disordered and inhomogeneous media and, in
22
general, in all nonthermal problems that depend on the geometry of random clusters
(res in forests, petroleum extraction, etc.)
A crucial step in the percolation problem is the analysis of clusters. There are many
algorithms to identify clusters on a lattice, some more ecient than others. Code P5 in
the course web page is an example.
2.3.d Numerical integration with many variables
Randomnumber generation can be easily used to estimate the value of integrals in mul-
tidimensional volumes. It is precisely this problem, that of averaging in multidimensional
volumes, that statistical mechanics is primarily concerned with. Traditional numerical
integration methods (based on the exact integration of families of orthogonal or non
orthogonal polynomials) cannot be used in multidimensional spaces. In order to under-
stand the reason, let us dene a square lattice in D dimensions. The number of sides of
this hypercube is 2D. Let n be the number of nodes along any of the D axes; then we
will have 2Dn
D1
nodes on the surface of the hypercube. But there are n
D
nodes in the
bulk (the remaining space of the hypercube, i.e. not counting the surface). The ratio of
the two is:
2Dn
D1
n
D
=
2D
n

D
Therefore the method becomes highly inecient since, as D increases, nodes on the sur-
face of the hypercube are sampled more and more often than nodes in the bulk.
The Monte Carlo does not present this diculty. Suppose we would like to numerically
estimate the value of the Ddimensional integral

V
D
g(x)d
D
x
where x (x
1
, x
2
, ..., x
D
) is a vector in D dimensions, d
D
x the corresponding volume
element, and V
D
is some particular hypervolume in D dimensions. Then the Monte
Carlo method writes

V
D
g(x)d
D
x
W
D
M
M

i=1
f(x
i
)
i
where W
D
V
D
is a volume that contains the volume V
D
(in particular they may be
the same volume), M is the number of uniform random numbers x
i
that we generate in
W
D
, and
i
is a number such that

i
=

1, if x
i
V
D
0, otherwise
If W
D
= V
D
(i.e. we are able to sample in V
D
) then
i
= 1 always. This method may
be regarded as a kind of rectanglebased method, where the rectangles have random sizes.
23
This method enjoys great advantages in case the integrand includes a probability dis-
tribution function, i.e. a positive-denite, normalised function f(x). To see this for
easily, we consider a one-dimensional problem. Consider the integral
I =

b
a
dxf(x)g(x) (2)
where f(x) is a probability distribution function. Then, we would evaluate the integral
by simply writing
I =

b
a
dxf(x)g(x)
b a
M
M

i=1
f(x
i
)g(x
i
)
where x
i
are uniformly distributed numbers in the interval [a, b]. But let us make use
of the cumulative function P(x):
y = P(x) =

x
a
dx

f(x

) dy = f(x)dx
If we make the change in variables x y, we have:
I =

b
a
dxf(x)g(x) =

1
0
dyg

P
1
(y)

But now, if we remind ourselves that, according to the acceptance-rejection method,

P
1
(y) is a random variable distributed according to f(x) provided y is uniformly dis-
tributed in [0, 1], we can write
I =

b
a
dxf(x)g(x)
1
M
M

j=1
g(x
j
) (3)
where x
j
are now a set of numbers distributed according to the function f(x) in the
interval [a, b]. We will make use of this result shortly. This latter result may be directly
extended to multidimensional integrals with an integrand weighted by a distribution func-
tion.
2.3.e Thermal problem: Ising model in 2D
Thermal problems are usually more complicated technically. The reason behind this
is related to the central problem of statistical physics. In a system with many degrees of
freedom (for N particles say 3N degrees of freedom, save contraints) in contact with a
thermal bath at temperature T, the probability of each state is given by the normalised
Boltzmann factor:
P e
H(q,p)
where
q q
1
, q
2
, ..., q
3N
, p p
1
, p
2
, ..., p
3N

and = 1/kT. We know from statistical physics that this distribution is highly peaked
about a few typical congurations; in fact, one way to obtain the Boltzmann factor
theoretically is to assume that only one conguration contributes, i.e. that conguration
24
Figure 11: A possible conguration of a 2DIsing model of 1/2 spins on a square lattice.
that maximises the degeneracy. In a nite system the distribution is wider, but even
then its width is narrow compared to the mean. Now uniform sampling will produce
congurations that, in most cases, will have a negligible statistical weight (i.e. small
Boltzmann factor), which will give rise to very inecient (or even helpless!) simulations.
To illustrate this important point, we will use an Ising model in 2D. Remember that the
Ising model is formulated on a lattice, in our case chosen as a square lattice, on every
node of which we place a 1/2spin, i.e. pointing either up or down (Fig. 11). The variable
associated to the ith node is s
i
= 1, with +1 meaning spin up and 1 spin down. The
Hamiltonian (energy) of the system is
H(s
i
) = J

nn
s
i
s
j
B

i
s
i
where J (here assumed positive) is a (ferromagnetic) coupling constant, and B an exter-
nal magnetic eld. nn stands for nearest neighbours, i.e. each spin interacts with their
four neighbours (up, down, right, left). s
i
, hereafter labelled simply as s, represents
a conguration of N spins, i.e., s = s
1
, s
2
, ..., s
N
. At zero magnetic eld, B = 0, this
model presents a phase transition from an ordered phase, with average magnetisation per
spin m = 0 for T < T
c
, to a disordered phase where m = 0 above T
c
. It is a second-order
phase transition, with a critical point at (B
c
, T
c
) = (0, 2....).
Our aim is the computation of thermal averages of quantities that can be dened for
each conguration s of the N spins; two obvious examples are the energy and the mag-
netisation. Let A(s) be one of such quantities. Then the statistico-mechanical expression
for the thermal average is
'A` =

s
A(s)P(s) =

s
A(s)e
H(s)

s
e
H(s)
(4)
In principle, the sums

s
are extended over all possible congurations of the N spins,
25
which are 2
N
in number. Even for very small N, the number of congurations is as-
tronomically large. For example, for N = 10
4
(a rather modest system on macroscopic
grounds),
2
N
= 2
10
4
10
3010.3
It is therefore out of question to sum over all congurations, and we should choose a
reduced set; the obvious solution is to randomly select the congurations to be included
in the sum, but in view of the presence of the weighting factor P(s), an importance
sampling scheme is in order. Using the same steps that led to (3) from (2), the average
(4) would be
'A` =
1
M

s
A(s), (5)
where M is the number of congurations generated. If, instead, we were to choose a set
of random uniform congurations, we would nd that most of them would have negligible
probability P(s) and would contribute negligibly to the sum (4). Since this is an inter-
esting and illustrative point, let us digress a little in this direction.
We will perform the above average, Eqn. (4), by two methods:
Uniform sampling. Here we obtain congurations by assigning spins to the nodes
of the lattice using the following algorithm:
1. We choose a uniform random number [0, 1]
2. If < 0.5, we set s
i
= +1; otherwise, we set s
i
= 1
3. We repeat for all spins i = 1, 2, ..., N
In this way we generate completely random congurations.
Importance sampling. The algorithm to be used here will be the Metropolis et al.
version, already introduced previously as a method to generate random deviates
distributed according to a particular probability distribution. The particular im-
plementation for the problem at hand will be explained in detail later. For the
moment we limit ourselves to presenting the results in comparison with those from
the previous uniform-sampling scheme.
Fig. 12 shows two energy histograms, f(e), with e = E/NJ the reduced energy
per spin, each obtained according to the two above sampling methodologies. The
histograms are constructed by computing the energy of each sampled conguration,
and assigning the energy to a given energy subinterval along the energy axis (the
distributions are normalised to unity).
The reduced temperatures were set to kT/J = 2.1 (upper panel) and 2.5 (lower
panel), which happen to be respectively below and above the critical temperature
T
c
of the phase transition (kT
c
/J 2.269). In the rst case the equilibrium phase
is therefore an ordered phase with nonzero magnetisation, whereas the second cor-
responds to a phase with zero magnetisation. The histogram based on the uniform
sampling method has zero mean and a small variance (as compared to the mean),
26
0
1
2
3
4
5
6
7
0
1
2
3
4
5
6
-2 -1.5 -1 -0.5 0 0.5 1
e
f
(
e
)
kT/J=2.1
kT/J=2.5
Figure 12: Energy distribution of a two-dimensional Ising model at two dierent reduced tem-
peratures (indicated as labels), below (upper panel) and above (lower panel) the critical tem-
perature (kT
c
/J 2.269). Energy e = E/NJ is expressed per spin and in units of the coupling
constant J. The shaded histogram was obtained by the uniform sampling method, whereas the
other corresponds to the true thermal distribution as obtained from Monte Carlo simulation
using importance sampling.
whereas that based on importance sampling has a mean correctly located at a neg-
ative value (we have to bear in mind that, in the perfectly ordered state T = 0,
the reduced energy per spin would be E/NJ = 2). The variance (i.e. width)
of both distributions are a little dierent, although both should be 1/

N (here
N = 20 20, which is a small number). What is remarkable is the negligible over-
lap between both distributions, the thermal and the uniform distributions. Uniform
sampling is not generating any single signicant conguration at all!
In the case where the temperature is above the phase transition temperature (T >
T
c
), the same situation holds, although this time the thermal distribution extends
over a much wider range of magnetisations, the mean one being zero, as is the case
with the uniform distribution. However, the latter misses a lot of congurations
with importante statistiscal weights.
As far as the distribution in magnetisation is concerned, Fig. 13 shows the distribu-
tions as obtained from the two sampling techniques. Note that, for the temperature
kT/J = 2.1, the temperature is below the critical temperature, and therefore two
(ideally symmetric) peaks appear, located at values of magnetisation diering in
sign (due to the short length of the Monte Carlo simulation with importance sam-
pling, the peaks are not completely symmetric). Again, the histogram obtained from
the uniform-sampling method does not overlap with the real, thermal distributions.
The average magnetisation obtained from this method would be zero, whereas we
27
0
1
2
3
4
5
6
7
8
0
1
2
3
4
5
6
7
-1 -0.5 0 0.5 1
m
f
(
m
)
kT/J=2.1
kT/J=2.5
Figure 13: Magnetisation distribution of a two-dimensional Ising model at two dierent re-
duced temperatures (indicated as labels), below (upper panel) and above (lower panel) the
critical temperature (kT
c
/J 2.269). Magnetisation m = M/N is given per spin. The shaded
histogram was obtained by the uniform sampling method, whereas the other corresponds to the
true thermal distribution as obtained from Monte Carlo simulation using importance sampling.
know that, since we are below the critical temperature, the equilibrium magnetisa-
tion is non-zero. Note that, strictly speaking, the importance-sampling distribution
also has zero mean, because it is bimodal and symmetric, but the correct average
magnetisation is the mean of one of the peaks
1
! This example demonstrates that it
is critical to use importance sampling in statistical simulations of condensedmatter
systems.
2.4. Monte Carlo code for the 2D Ising model
Now that we have introduced the Ising model in 2D and discussed the signicance of
using an importance-sampling technique to generate congurations for performing the
averages, we will discuss how to implement the Metropolis et al. algorithm for the Ising
model. Again, we repeat that this method is the only one feasible in statistical mechanics,
since it does not require knowledge of the partition function (norm of Boltzmann weight),
which cannot be calculated.
Before discussing the code, we give a few general comments on the basis of the method,
which can also be extended to other applications.
Initialisation. Ideally one should start from an initial conguration which is repre-
1
For temperatures T well below the critical temperature this eect does not appear since any reason-
ably long simulation, even long enough to correctly explore congurations with, say positive magnetisation
this will depend on the initial conguration, will not be that long as to explore states with opposite
magnetisation, and correct means will result.
28
-1
-0.5
0
0.5
1
-1
-0.5
0
0.5
1
5 10 15 20 25 30 35 40 45 50
MC time
M
/
N
M
/
N
kT/J=2.5
kT/J=1.8
Figure 14: MC time evolution of magnetisation per spin, M/N, at two dierent tempera-
tures: kT/J = 1.8 (upper panel), and kT/J = 2.5 (lower panel). In both cases the starting
congurations were chosen as allup, random and alldown.
sentative of the thermodynamic conditions. This is not always possible. In the Ising
model, we may start from allup, alldown, or random congurations. Clearly, if
the temperature is below the critical temperature, then the rst two congurations
are to be preferred. Fig. 14 (upper panel) shows how the magnetisation per spin
evolves in MC time or MC steps (see later) for kT/J = 1.8: the initial random
conguration, with very small magnetisation per spin, evolves towards congura-
tions with high magnetisation, which are more representative of the equilibrium
state, but it takes some time for the Monte Carlo algorithm to reach this situation
(typically we have to wait for a relaxation time , which will become longer as the
temperature is decreased). Allup or alldown congurations, on the other hand,
are more representative of the equilibrium state since they have a higher Boltzmann
factor, and they seem to converge rapidly towards equilibrium. For kT/J = 2.5
(Fig. 14, lower panel), which corresponds to the disordered phase, it is now the all
up and alldown congurations the ones that are not representative, so that they
show a slow relaxation towards equilibrium (the initial random conguration also
shows slow relaxation; this is due to the close proximity of the phase transition).
Warming up. As we know, the Metropolis et al. algorithm tends to generate
correctly distributed congurations, but only asymptotically. A warmup period is
needed. Depending on the quality of the initial conguration, this period may be
short or long (for example, we may have a look at how the magnetisation or the
energy are distributed, and wait until Gaussian functions are obtained).
29
Trial moves. In the usual Monte Carlo method for the Ising model, spins are
moved one at a time. Normally one visits the spins sequentially, and moves the
visited spin (from up to down or the other way round). Say we are on the ith spin,
with variable s
i
. Then the change of energy involved in turning the spin is
E = 2s
i

n.n.
s
j
where n.n. are the four neighbour spins (at right, left, on top and bottom). The
probability with which this singlespin move is accepted is taken as exp (E).
After we have attempted to turn the N spins, we have one MC step, so that M
steps involve M attempts per spin. Trial moves may involve more than one spin
(for instance, blocks of spins). However, this methodology may not work, as energy
changes involved are proportional to the number of spins involved in the move, and
the corresponding acceptance probability may be too small.
Acceptance criterion. In the Ising model there is only one possibility to move
the spin. In other models (e.g. models with continuous variables) we have a range
of possible moves to generate a test value for the random variables. In this latter
case we have to adjust the possible moves so as to have an acceptance ratio of
approximately 50%, since this ratio optimises the rate at which equilibrium (i.e.
the asymptotic regime) is reached.
001. !
002. ! P6: 2D-Ising model Monte Carlo with
003. ! importance sampling
004. !
005. parameter (n=20)
006. integer*1 is(n,n)
007. real*8 dseed
008. dimension histm(101),histe(101)
009.
010. data dseed /186761d0/
011. nsteps = 100000000
012. nblock=10
013.
014. n2 = n*n
015. T=2.1
016.
017. ama = 0
018. amaa = 0
019. ene = 0
020. enea = 0
021.
022. do i=1,101
023. histm(i)=0
024. histe(i)=0
025. end do
026.
30
027. do i = 1, n
028. do j = 1, n
029. is(i,j)=1
030. call ran (dseed,random)
031. if(random.lt.0.5) is(i,j)=-1 ! comment to set ordered state
032. ama = ama + is(i,j)
033. end do
034. end do
035.
036. do i = 1, n
037. ip = i + 1
038. if (ip .gt. n) ip = 1
039. do j = 1, n
040. jp = j + 1
041. if (jp .gt. n) jp = 1
042. ene = ene - is(i,j)*(is(ip,j)+is(i,jp))
043. end do
044. end do
045.
046. do k = 1, nsteps/nblock
047. do l = 1, nblock
048.
049. do i = 1, n
050. do j = 1, n
051. sij = is(i,j)
052. ip = i + 1
053. if (ip.gt.n) ip=1
054. im = i - 1
055. if (im.lt.1) im=n
056. jp = j + 1
057. if (jp.gt.n) jp=1
058. jm = j - 1
059. if (jm.lt.1) jm=n
060. de = 2*sij*(is(ip,j)+is(im,j)+is(i,jp)+is(i,jm))
061. call ran (dseed,random)
062. if (exp(-de/T) .gt. random) then
063. is(i,j) = -sij
064. ama = ama - 2*sij
065. ene = ene + de
066. end if
067. end do
068. end do
069.
070. end do
071.
072. enea = enea + ene
073. amaa = amaa + ama
31
071.
071. im=(ama/n2+1)*0.5e2+1
071. histm(im)=histm(im)+1
074.
075. ie=(ene/n2+2)*50.0+1
076. histe(ie)=histe(ie)+1
077. end do
078.
079. write(*,(T= ,f5.3, M/N= ,f5.3, E/N= ,f6.3))
080. >T, abs(amaa/(nsteps/nblock*n2)),enea/(nsteps/nblock*n2)
081.
082. open(3,file=histe.dat)
083. sum=0
084. do i=1,101
085. sum=sum+histe(i)
086. end do
087. do i=1,101
088. emm=2*(i-1)/100.0-2
089. write(3,(2f9.3)) emm,histe(i)/sum*50.0
090. end do
091. close(3)
092.
093. open(3,file=histm.dat)
094. sum=0
095. do i=1,101
096. sum=sum+histm(i)
097. end do
098. do i=1,101
099. amm=2*(i-1)/100.0-1
100. write(3,(2f9.3)) amm,histm(i)/sum*50.0
101. end do
102. close(3)
103.
104. end
Values of the spins are stored in the two-index matrix is(n,n), where n is the num-
ber or rows and columns of the lattice (so that n2, dened on line 014 is the number of
spins). The temperature T is set on line 015. Accumulators for magnetisation and energy,
amaa and enea, are set to zero, as well as variables for the magnetization and energy
histograms (lines 022025). On lines 027034 the spins are assigned, either in all-up con-
guration or randomly (this is changed simply by commenting or uncommenting line 031),
and the magnetisation of this initial conguration is computed and stored in variable ama.
Lines 036044 compute the energy of this conguration. On lines 046 and 047 the loop
over MC steps is initiated. The loop is divided into blocks of nblock steps, so that there
are nsteps/nblock blocks (nsteps is the number of total MC steps in the simulation).
This is made so that quantities can be accumulated not every single MC step but every
nblock steps to avoid correlation eects between congurations. On lines 049-068 spins
32
are chosen sequentially, one at a time, and switched from up to down or the other way
round. The change of energy involved is computed in de, on line 060, and on line 062
the move is either accepted or rejected depending on the Metropolis et al. criterion. If
accepted, the magnetisation and energy of the new conguration are computed on lines
064 and 065, respectively. Lines 072073 accumulate these quantities every block, and
also the histograms are built on lines 071076. The end do instruction on line 077 ends
the main loop over MC blocks. After that, results are written; note that, whereas the
temperature, magnetisation per spin and energy per spin are dumped on the screen (lines
079 and 080), the histograms are written on les histe.dat and histm.dat for energy
and magnetisation, respectively. Also note that the warm-up stage is not included in
this version of the code.
The results obtained from running this code are shown in Figs. 15. In Fig. 15(a) the
evolution of the magnetisation per spin and the energy for some particular experiment
(with an allup initial conguration) is displayed. The temperatures are kT/J = 1.80 and
2.19. The second one is particularly interesting, as it is close to the critical temperature
(but still below T
c
), and it shows that the system visits congurations with positive and
negative values of M/N. The system has been caught in the act: in a matter of 30
MC steps all the spins are ipped so that M/N changes sign. This behaviour is typical
of states close to critical points.
-2.0
-1.5
-1.0
-0.5
0.0
0.5
1.0
0 100 200 300 400 500 600 700 800
M/N
E/N
MC steps
kT/J=1.80
kT/J=2.19
E
/
N
,

M
/
N
-2.0
-1.5
-1.0
-0.5
0.0
0.5
1.0
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5
M/N
E/N
kT/J
E
/
N
,

M
/
N
(a) (b)
Figure 15: (a) MC time evolution of magnetisation per spin, M/N, and energy per spin, E/N,
at two dierent temperatures. (b) Phase diagram (M/N vs. T) and dependence of energy per
spin, E/N, on T.
Fig. 15(b) is the phase diagram (M/N vs. T). The number of MC steps used to obtain
the magnetisation for each value of the temperature was 10
7
. We obtain a smooth function
of the temperature: due to the nite size of our system (N = 20 20 spins) the disconti-
nuity in the derivative of M/N at the critical temperature is smoothed out. The critical
temperature can be obtained using some practical criterion (e.g. looking for the temper-
ature for which M/N = 0.5), but there are rigorous ways to extract information from
the MC data (for example, the method based on cumulants) and obtain more meaningful
information, such as the value of T
c
. The energy per spin, E/N, is also plotted in the
gure. Note that it is a smooth function of temperature, and always negative, even in the
33
disordered phase; this means that there is some degree of local ordering even in this phase.
2.5. Model with continuous variables: hard spheres
Now we would like to illustrate the Monte Carlo technique in the context of a thermal
problem with continuous degrees of freedom. We have chosen the hardsphere model,
since it is both simple and suciently important. The model has played a crucial role in
the development of the statistical mechanics of uids and solids, and as a consequence, it
is a central model in condensedmatter physics.
The model approximates the interaction potential between two spherical molecules, (r),
by:
(r) =

, r <
0, r >
is the diameter of the hard spheres; these are termed hard because they cannot penetrate
each other. r is the distance between the centres of mass of the two spheres. The model
is an attempt to approximate the repulsive part of the interaction potential in simple
systems (spherically symmetric molecules, either in reality or eectively), for instance
noble gases, some metals, CH
4
, NH
3
, ... It is an athermal model (i.e. thermodynamic
properties depend trivially on temperature) since the corresponding Boltzmann factor
does not depend on temperature T:
e
(r)
=

0, r > (spheres overlap)

1, r > (spheres do not overlap)
Therefore, the only allowed congurations of N spheres are congurations where there is
no overlap between the spheres; all these congurations have the same statistical weight,
and possess zero energy.
Despite its trivial appearance, the hardsphere model contains highly nontrivial physics.
The only relevant parameters are density = N/V and pressure p. The former is more
conveniently represented by the socalled packing fraction, which is the ratio of volume
occupied by the N spheres to total volume V . If v = (4/3) (/2)
3
=
3
/6 is the
molecular volume, then the packing fraction, , is:
=
Nv
V
= v =

6

3
In terms of , the model has (Fig. 16):
A uid phase in the interval [0, 0.494], with longrange positional disorder (there
is, however, shortrange order; this will be discussed later on).
A crystalline phase, in the interval [0.545, 0.740], with longrange positional
order. The crystal is a fcc (face-centred cubic) lattice. Note that the maximum
possible value of the packing fraction is 0.740.
34
Figure 16: Phase diagram of hard spheres in terms of the packing fraction.
An amorphous (metastable) phase, in the interval [0.59, 0.64], with longrange
positional disorder but a negligible viscosity coecient and virtually zero diusion.
There is strongly rstorder phase transition between the uid and the crystal, at a re-
duced pressure p

p
3
/kT 8.9
Before presenting the Monte Carlo algorithm as implemented for hard spheres, we in-
troduce a method to calculate the pressure, along with two simple theories, applicable to
the uid and the crystal, respectively, that describe in simple (albeit approximate) terms
the equation of state of the system, i.e. the function p(). In the theory for the crystal
phase calculations are most easily done using the (unweighted, i.e. uniform sampling)
Monte Carlo integration technique; in the former, it may also be necessary to use it in
case renements be required.
Pressure and radial distribution function. The radial distribution function
g(r) is a function of the distance between a pair of particles. It is a measure of
the local order of the uid. In practice it is computed as the mean number of
particles around a given one at distance r, divided by the expected number should
interactions be turned o (the latter is just the bulk density ). For a dense uid,
the function shows a strong rst peak at distances , and a damped oscillatory
structure tending to one as r increases. As the density is increased, this structure
becomes more and more sharpened. Knowledge of g(r) for a uid allows dierent
thermodynamic properties to be computed. For example, the reduced pressure is
given by
p
kT
= 1
2
3kT

0
drr
3

(r)g(r) (6)
where

(r) is a derivative with respect to r. This is valid for any spherically

symmetric interaction potential. In particular, for the hardsphere potential (where

(r) exhibits deltatype behaviour at r = ), the above expression can be written

in terms of the value of g(r) at r = (the socalled contact value):
p
kT
= 1 +
2
3

3
g(
+
)
35
In the simulation, g(r) is computed as a histogram; we will explain this later.
Clausius equation of state for the uid. We start from the idealgas equation
of state, which we known is exact in the lowdensity limit = N/V 0:
pV = NkT
Here p is the pressure, V the volume, N the number of particles, k Boltzmann
constant, and T the temperature. At higher density the volume occupied by the
spheres has to be taken into account, and the equation must be modied as
p(V Nb) = NkT p =
NkT
V Nb
where b is the excluded volume per particle. This volume arises from the fact that
each sphere has an associated volume within which the centres of mass of other
spheres are not allowed. Every pair of spheres contributes to this volume with
b =
1
2

4
3

3

=
2
3

3
= 4v
Then:
pV
NkT
=
1
1 Nb/V
=
1
1 4
which is Clausius equation. This equation predicts the following virial expansion:
pV
NkT
=
1
1 4
= 1 + 4 + 16
2
+ ...
The coecient of the term
2
is incorrect; actually its correct value is 10. The
predicted higherorder coecients are also in error. In fact, the virial coecients,
B

n
= B
n
/v
n1
, dened by the virial expansion
pV
NkT
= 1 +

n=1
B
n+1

n
=

n=1
B

n+1

n
can be obtained from Clausius equation as B

n
= 4
n
; scaling with the second virial
coecient, we get
B
n
B
n1
2
= 1, n > 1
which is a very gross approximation. In particular, Clausius equation predicts a
divergence of the pressure at a packing fraction = 1/4 = 0.250, which is clearly
wrong since we know that the uid is stable at least up to = 0.494. Clausius
equation is approximate, and can be used only at low densities.
The free energy can be obtained from the pressure by integration; the nal result is
F
NkT
= log

3
1 4

1 (7)
Clausius equation can be improved in various ways. Much eort has been devoted
over the last decades to this goal. A typical approach incorporates knowledge of
36
higher-order virial coecients to build up a more accurate equation of state. One of
these approaches leads to the socalled CarnahanStarling equation of state which,
when properly integrated in density, gives the following free energy per particle:
F
NkT
= log

1 +
(4 3)
(1 )
2
(8)
A more general approach involves constructing Pade approximants. For example, if
the rst two virial coecients, B
2
and B
3
, were known, the coecients a
1
and a
2
in
the following Pade approximant,
pV
NkT
=
1 + a
1
+ a
2

2
1 v
= 1 + B
2
+ B
3

2
+ ...
can be easily obtained. The approximant is only adjusted for low densities, but its
functional form is expected to be valid for a much larger density range.
Expressions for the virial coecients, based on statistical mechanics, can be written.
They include the socalled Mayer function, f(r) = exp [(r)] 1. The rst two
coecients read:
B
2
=
1
2

drf(r)
B
3
=
1
3

dr

dr

f(r)f(r

)f ([r r

[) (9)
Expressions for the higher-order coecients are considerably more complicated. For
hard spheres, the rst three coecients are known exactly:
B
2
=
1
2

dr( r) = 2

0
drr
2
=
2
3

3
= 4v
B
3
=
5
2
18

6
2.74156
B
4
=

89
280
+
219

2
2240
+
4131
2240
arccos

B
3
2
(10)
The higher-order coecients, B
n
with n > 4, must be calculated numerically. Note,
in passing, that Clausius equation predicts B
3
= 16v
2
= 16 (
3
/6)
2
= 4.3865,
to be compared with the exact value 2.7416; not a fair comparison really!
As an example of how the higher-order virial coecients can be evaluated numer-
ically, let us evaluate the B
3
coecient (which, as mentioned, is known exactly)
using a technique based on the multidimensional Monte Carlo integration. We use
Eqn. (9), which extends over all space, but only when three spheres overlap at the
same time does the integrand contribute. In spherical coordinates:
B
3
=
1
3

0
drr
2

1
1
d (cos )

2
0
d

0
dr

r
2

1
1
d (cos

2
0
d

f(r)f(r

)f ([r r

[)
37
M B
MC
3
B
MC
3
/B
exact
3
10
2
2.51933 0.91894
10
3
2.65935 0.97002
10
4
2.77402 1.01184
10
5
2.74633 1.00174
10
6
2.74411 1.00093
10
7
2.74242 1.00032
10
8
2.74170 1.00005
Table 1: Third virial coecient as obtained from the Monte Carlo integration method for
dierent number of congurations generated, M, and ratio to the exact value, given by Eqn.
(10).
Note that the jacobians in coordinates have been absorbed into dierentials, so
that it is cos that is to be sampled uniformly. Let us place a sphere at the origin,
and choose a location for a second sphere randomly so that it overlaps with the rst.
This is easily done if we randomly and uniformly generate numbers (r, cos , ) in
the intervals [0, ] [1, 1] [0, 2]. Then, one of the Mayer functions, say f(r),
is always equal to minus one. If we do the same for a third sphere, i.e. randomly
generate a position so that it always overlaps with the rst, we have f(r

) = 1.
Now there would remain to check whether the second and third spheres overlap, i.e.
whether the value of the Mayer function f ([r r

[) is 0 or 1. The approximation
based on Monte Carlo integration would, in this case, be
B
3
=
1
3

1
M
[() (2) (2)]
2
M

i=1
r
2
i
r
2
i

i
=
(4)
2
3M
M

i=1
r
2
i
r
2
i

i
where i = 1, ..., M is a Monte Carlo step (which includes sampling two spheres within
the excluded volumen of the central one) and = 1 or 0 depending on whether the
two spheres overlap. Table 1 contains results obtained with dierent number of MC
steps generated. As can be seen, the results for B
3
converge nicely to the exact
value. This technique is obviously useful beyond the fourth virial coecient, for
which there are no exact results.
Freevolume equation for the crystalline phase. For the other stable phase,
the crystal, a simple approximation, called freevolume approximation, can be made.
Let us consider a crystalline lattice, with the spheres vibrating about their equilib-
rium positions, which coincide with the sites of the lattice. The freevolume approx-
imation assumes the neighbours of a given sphere to be xed in their equilibrium
positions, and that our central sphere moves in the cage left by the neighbours. The
problem is then that of an ideal particle in an external potential; the latter is innite
whenever the particle overlaps with the neighbours, and zero otherwise. This gives
rise to a free volume, v
f
, and the partition function of the particle will be z = v
f
.
The partition function for the N spheres (which we assume to be distinguishable,
38
since a particular sphere has to be chosen to calculate v
f
) will be
Z =
1
h
3N
N

i=1

v
f
dr
i

dp
i
=
v
N
f

3N
where is the thermal wavelength (this comes from the integration over momenta).
The free energy will be:
F = kT log Z = NkT log

v
f

,
F
NkT
= log

3
v
3
f

and the pressure:

p =

F
V

The dependence of F on V (or density ) comes through v

f
, which depends on how
close or how far the neighbours of the central sphere are, which in turn depends on
density.
To calculate v
f
we can make approximations, or we can use the Monte Carlo in-
tegration method, in the same manner as we used it to compute the third virial
coecient B
3
. In Fig. 17 the situation is explained on a triangular lattice. Six
neighbours of a central sphere (the latter is not shown) are depicted in light grey.
The region where the centre of mass of the central particle could be located without
ever overlapping with its neighbours is also depicted (the region is bounded by cir-
cular sectors, each centred on each of the six neighbours and with radius 2). On a
two-dimensional triangular lattice the free volume v
f
can be obtained analytically.
The three-dimensional case requires considerably more eort, so that Monte Carlo
integration is well suited to the problem: we would generate positions for the central
particle within a volume W that contained the volume v
f
we want to calculate; any
point generated such that there is overlap with any of the neighbouring spheres is
rejected (or given zero weight). The free volume would be:
v
f

W
M
M

m=1

i
,
i
=

1, no overlap
0 overlap
We have to make sure that the volume W is suciently large as to include the free
volume. The quantity
i
is easily calculated by simply generating the positions of
the n nearest neighbours of the lattice in question. For example, for the fcc lattice
n = 12, and their positions are given by:

a
2
,
a
2
, 0

a
2
, 0,
a
2

0,
a
2
,
a
2

where a is the lattice constant (related to density by = 4/a

3
). This model, among
other things, explains why the stable crystalline phase is a fcc lattice and not, for
example, the bcc body-centred cubic lattice.
39
Figure 17: The six neighbours of a central sphere (the latter is not shown) depicted in light
grey. The region where the centre of mass of the central particle could be located without ever
overlapping with its neighbours is also depicted.
Phase transition between a uid and a crystal in hard spheres. Fig. 18
shows the free energies of the uid and two crystal structures for the hard-sphere sys-
tem. The rst one is obtained from the Carnahan-Starling approximation, whereas
the latter two were obtained using the free-volume approach. Note that the fcc phase
is always more stable than the bcc phase, which is metastable (actually computer
simulations have proved it to be unstable against shear modes, not considered in the
simple freevolume approach). Also note that the uid branch has been computed
even for high densities, for which the uid no longer exists (this could be considered
as a highly metastable uid). The uid and fcc curves intersect at some density.
This demonstrates the uidtosolid transition is of rst order. In order to obtain
the densities (or packing fractions) of the coexisting uid and crystal phases, we
would have to apply a commontangent construction, i.e. obtain
f
and
c
for uid
and crystal by solving the equations p
f
= p
c
and
f
=
c
simultaneously ( is the
chemical potential). Results from Monte Carlo estimations of free energies give the
following values for the coexistence densities:
f

3
= 0.943 and
c

3
= 1.041, which
correspond to packing fractions
f
= 0.494 and
c
= 0.545.
FORTRAN codes and results
In this section we present complete FORTRAN codes to implement various calculations
presented before, and discuss some of the results.
Calculation of B
3
coecient
The code used is the following:
001. !
002. ! P7: Evaluation of B3 by Monte Carlo integration
003. ! hard-sphere diameter sigma=1
004. !
40
0
2
4
6
8
10
12
14
16
18
20
0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70
F
/
V
k
T

fcc
bcc
fluid
Figure 18: Free energy density per unit thermal energy F/V kT for the uid phase (using the
Carnahan-Starling approximation), and for two crystalline phases (fcc and bcc). For the latter
two the freevolume approximation was used.
005. implicit real*8(a-h,o-z)
006. real*4 r,t,f
007. integer*8 M,i
008.
009. dseed=173211d0
010. pi=4.0d0*datan(1.0d0)
011. pi2=2*pi
012.
013. do mpot=2,8
014. M=10**mpot
015. sum=0
016. do i=1,M
017. call ggub (dseed, r)
018. call ggub (dseed, t)
019. call ggub (dseed, f)
020. cp=dcos(pi2*f)
021. sp=dsin(pi2*f)
022. ct=2*t-1
023. st=dsqrt(1-ct**2)
024. x=r*st*cp
025. y=r*st*sp
026. z=r*ct
027. r1a=x**2+y**2+z**2
028. call ggub (dseed, r)
41
029. call ggub (dseed, t)
030. call ggub (dseed, f)
031. cp=dcos(pi2*f)
032. sp=dsin(pi2*f)
033. ct=2*t-1
034. st=dsqrt(1-ct**2)
035. x1=r*st*cp
036. y1=r*st*sp
037. z1=r*ct
038. r1b=x1**2+y1**2+z1**2
039. r2=(x-x1)**2+(y-y1)**2+(z-z1)**2
040. if(r2.lt.1.0d0) sum=sum+r1a*r1b
041. end do
042. vol=4*pi
043. B3=sum*vol**2/(3*M)
044. B3ex=5*pi**2/18
045. write(*,(M=,i9, B3=,f8.5, B3/B3ex=,f8.5))
046. > M,B3,B3/B3ex
047. end do
048. end
The integration is done with dierent numbers of Monte Carlo steps, in order to
visualise the convergence. Therefore, there is a loop in the variable mpot, which
gives the number of steps in powers of 10. x, y, z are the random coordinates of
the second sphere, which is chosen within the excluded volume of the rst (placed
at the origin). x1, y1, z1 are the random coordinates of the third sphere, also
chosen to always overlap with the rst. In lines 039 and 040 the distance between
the second and third spheres is checked, and if less than one ( is set to unity)
the conguration is accepted. vol in line 042 is the volume sampled per sphere (so
in total a volume vol**2 is sampled). In line 043 the third virial coecient B3 is
obtained and then compared with the exact result B3ex.
Calculation of free volume v
f
This code is prepared for the fcc lattice.
001. !
002. ! P8: Free volume for the fcc lattice
003. !
004. implicit real*8(a-h,o-z)
005. dimension xx(12),yy(12),zz(12)
006. data xx/+0.0d0,+0.0d0,+0.0d0,+0.0d0,+0.5d0,-0.5d0,+0.5d0,
007. >-0.5d0,+0.5d0,-0.5d0,+0.5d0,-0.5d0/
008. data yy/+0.5d0,-0.5d0,+0.5d0,-0.5d0,+0.0d0,+0.0d0,+0.0d0,
009. >+0.0d0,+0.5d0,+0.5d0,-0.5d0,-0.5d0/
010. data zz/+0.5d0,+0.5d0,-0.5d0,-0.5d0,+0.5d0,+0.5d0,-0.5d0,
011. >-0.5d0,+0.0d0,+0.0d0,+0.0d0,+0.0d0/
012. real*4 r
42
013. pi=4d0*datan(1d0)
014. sq2=dsqrt(2d0)
015.
016. rho=1.3d0
017. a=(dsqrt(2d0)/rho)**(1d0/3d0)
018. eta=rho*pi/6
019. dseed=173211d0
020. M=2e7
021. sum=0
022.
023. do i=1,M
024. call ggub (dseed, r)
025. x=a/2*(2*r-1)
026. call ggub (dseed, r)
027. y=a/2*(2*r-1)
028. call ggub (dseed, r)
029. z=a/2*(2*r-1)
030.
031. do k=1,12
032. xk=a*xx(k)*sq2
033. yk=a*yy(k)*sq2
034. zk=a*zz(k)*sq2
035. r2=(xk-x)**2+(yk-y)**2+(zk-z)**2
036. if(r2.lt.1d0) go to 1
037. end do
038. sum=sum+1
039. 1 end do
040.
041. vf=sum/M*a**3
042. write(*,(eta=,f10.5, vf=,f10.5, ff=,f12.7))
043. >eta,vf,-rho*dlog(vf)
044.
045. end
The coordinates of the twelve neighbours of the central sphere, in the fcc lattice,
are dened in vectors xx, yy and zz (lines 005011). The coordinates are in units
of the cubic lattice parameter. The density rho is set in line 016, and the nearest
neighbour distance a is obtained in line 017. eta is the packing fraction, and M (line
020) the number of MC steps (i.e. the number of random points generated). Now
a loop is opened in line 023, and coordinates for the random point x, y and z are
chosen within a cube of side a centred on the central site. Then in lines 031037
it is checked whether the generated position for the central sphere is such that an
overlap occurs. If this is the case, the point is rejected (line 036). If not, an ac-
cumulator is advanced by one (line 038); this is the variable sum, which after the
M steps will contain the number of points that lie within the free volume. This is
calculated in line 041, and in lines 042043 the packing fraction, the free volume
and the associated free energy per unit volume and unit thermal energy are writ-
43
ten on the screen (the thermal wavelength is set to one for the sake of convenience).
Monte Carlo code for thermal simulation of hard spheres
This FORTRAN code implements the MC simulation of thermal averages for the hard-
sphere system using the importancesampling technique. The code performs the
MC dynamics, and calculates the radial distribution function from which, in a (sep-
arate) nal calculation, the pressure is evaluated. The results of the simulations
obtained with this code can be compared with the approximate equations of state
for uid and crystal phases obtained in the previous sections. Before presenting the
code we will digress a little on the peculiarities of the Monte Carlo algorithm for
the system of hard spheres.
In Section 2 we presented the Metropolis et al. algorithm, and extended it to
more than one random variable. In the case of N hard spheres contained in a rect-
angular box, a move consists of randomly moving a sphere and checking whether
this displacement gives an overlap or not (Fig. 19). If there is no overlap, the move
is accepted; otherwise it is rejected. This is so because the acceptance probability
in this case is simply
e
E
=

0, there is at least one overlap

1, there are no overlaps
Figure 19: The chosen sphere (in gray) is given a random displacement to a test position (in
gray line). The neighbouring spheres are represented as circles in black line. In this case the
displacement results in no overlap, and the move will be accepted.
The algorithm is then as follows:
1. We prepare an initial (nonoverlapping) conguration for the N spheres within
the rectangular box, at the density we wish to investigate (a highdensity
44
disordered conguration is dicult to obtain; that is why normally we will
start from a perfectly ordered conguration, even if the density is too low for
the crystal phase to be stable).
2. We choose a sphere, at random or sequentially, and displace it by a random
amount (within some distance from the original position)
3. We check for overlaps with neighbouring spheres. If there is at least one overlap,
the move is rejected; if there are no overlaps, the move is accepted.
4. We repeat the process until the N spheres have been visited (or N moves have
been performed). Then we will have completed one MC step.
5. We iterate the procedure as many times as needed
As in any implementation of the Metropolis et al. algorithm, an initial warmup
period is required before taken averages. Also, the ratio of accepted moves to total
moves should be chosen to be about 0.5, and this is realised by adjusting the value
of the maximum displacement (clearly, if is too large moves will mostly result
in overlaps, and most moves will be rejected, while if is very small the moves will
almost always be accepted).
i
j j
j j
j j j
j
j
Figure 20: Central simulation box surrounded by replicas (in the case of twodimensional sys-
tems and square boxes there are 8 replicas; for threedimensional systems there are 27 replicas).
Particles in the replicas have the same positions and velocities as in the central box; then, when
exiting one box they supersede their clon that, in turn, is exiting its box (periodic boundary
conditions). Distance between the ith and jth particles is calculated as the minimum distance
between the ith and all the clones of the jth particles, named j

in the gure (minimum image

convention). In the example, the minimum distance corresponds to one of the clones of j (the
one right above the central box), not to j itself.
There are some technicalities of the method that we discuss in the following.
a. Periodic boundary conditions.
When simulating extended (bulk) systems it is convenient to use periodic
45
boundary conditions (Fig. 20). These boundary conditions minimise the ef-
fects of surfaces on the results, and consist of placing replicas or images of the
system next to it and along all directions. Let L
x
, L
y
, L
z
be the sides of the
(rectangular) simulation box. In practice, the method redenes the coordinates
of any particle that gets outside of the box; for example, if at some instant of
time the x coordinate of a particle is x > L
x
, then we make x L
x
x. If
x < 0, then we redene x +L
x
x. This is done on every coordinate of every
particle along the simulation.
b. Minimum image convention.
For pairwise interaction potentials, we have to calculate the distance between
the two particles of a given pair. For rectangular simulation boxes the sim-
ulation box has 26 replicas. This means that, when calculating the distance
between the ith and the jth particles, one has decide which of the 27 possible
distances between the ith particle and the jth particles, including itself and
its 26 replicas, has to be taken into account, since all or some of them can be
within the interaction range. The minimum image convention chooses the one
which gives the least distance [r
j
r
i
[ (Fig. 20). See the code for a possible
algorithm that implements this.
c. Calculation of radial distribution function.
We focus on a particular sphere. Then g(r)d
3
r is the number of particles to be
found in a dierential volume element d
3
r a distance r from that sphere, with
respect to the number that would obtain should interactions between particles
be turned o (or equivalently, with respect to the corresponding number that
would obtain very far from the sphere, where interactions and for that matter
correlations are zero or negligible). In practice, given the expected spherical
symmetry of the problem, we consider a spherical dierential volume element
of width r around a particular particle i (see Fig. 21). Put in a mathematical
expression:
g(r) =
1
N

i=1
n
i
(r, r)

4r
2
r
The average '...` is over dierent congurations of the spheres. Note that
an average over particles is also explicit in the expression (sum over particles,
i = 1, 2, ..., N) and division by N; this can (and should) be done, as all particles
are equivalent and it helps improve the accuracy. n
i
(r, r) is the number of
particles within a spherical shell of width r, centred on the ith particle, and
at a distance r from this particle. This can be easily computed for any given
congurations of sphere. Finally, 4r
2
r is the volume of such a spherical
shell, which can also be approximated by
4
3

r +
r
2

r
r
2

When multiplied by , we get the number of particles within the volume at

very large distance (rigorously speaking, we have to multiply by (N 1)/V ,
not by = N/V , since the central particle cannot be counted). Then, for
46
large r, the ratio should tend to unity. In the simulation a number of discrete
distances r
i
is sampled, and g(r) is constructed as a normalized histogram.
i
r
r
Figure 21: The radial distribution function at distance r is calculated by counting the number
of particles in a spherical shell of radius r (here represented by a ring) and thickness r when
the origin is at the ith particle, giving n
i
(r, r). This is done for all the N particles and then
averaged (dividing by N). The resulting quantity is normalised to the number of particles
expected in the shell if there were no interactions, i.e. 4r
2
r.
We now present the Monte Carlo code for hard spheres. The code has a main body
of code and a few subroutines. Subroutine ggub was explained previously so it is
not included here.
001. !
002. ! P9: MC simulation of hard spheres
003. !
004. implicit real*8(a-h,o-z)
005. parameter (num=4000, ngrx=1000)
006.
007. real*4 rr
008. dimension r1x(num), r1y(num), r1z(num)
009.
010. common gr(ngrx)
011.
012. data dseed /1467383.d0/
013. data hr /0.05d0/
014. data pi /3.141592654d0/
015.
016. open (1, file = hs.par)
017. read (1,*) ini, lmn
018. read (1,*) nstepeq, nstep0, nblock0
019. read (1,*) dens, amp
47
020. close (1)
021.
022. if (lmn .gt. num) stop lmn .gt. num
023.
024. open (2, file = hs.inp)
025. open (3, file = hs.out)
026. open (4, file = hs.gr)
027.
028. if (ini .eq. 0) then
029. call fcc (lmn, dens, r1x, r1y, r1z, xlx, yly, zlz)
030. else
031. do i = 1, lmn
032. read(2,*) r1x(i), r1y(i), r1z(i)
033. end do
034. read(2,*) xlx,yly,zlz
035. end if
036.
037. xll = xlx/2
038. Yll = yly/2
039. zll = zlz/2
040.
041. volume = xlx*yly*zlz
042. ro = lmn/volume
043.
044. nox = xll/hr
045. if (nox .gt. ngrx) stop nox .gt. ngrx
046.
047. do istage = 1, 2
048.
049. if (istage .eq. 1) then
050. nblock = nstepeq
051. nstep = nblock
052. else
053. nblock = nblock0
054. nstep = nstep0
055. end if
056.
057. do j = 1, lmn
058. call pbc (r1x(j),r1y(j),r1z(j),xlx,yly,zlz)
059. end do
060.
061. do j = 1, ngrx
062. gr(j) = 0
063. end do
064.
065. naccept = 0
066. numgr = 0
48
067.
068. do ii = 1, nstep/nblock
069.
070. numgr = numgr + 1
071. call gdr (lmn,r1x,r1y,r1z,xlx,yly,zlz)
072.
073. nacceptb=0
074.
075. do jj=1,nblock
076.
077. do j = 1, lmn
078. call pbc
079. > (r1x(j),r1y(j),r1z(j),xlx,yly,zlz)
080. end do
081.
082. do k = 1, lmn
083. call ggub(dseed, rr)
084. desx = (rr-0.5)*amp
085. call ggub(dseed, rr)
086. desy = (rr-0.5)*amp
087. call ggub(dseed, rr)
088. desz = (rr-0.5)*amp
089. r1xk = r1x(k)
090. r1yk = r1y(k)
091. r1zk = r1z(k)
092. r3x = r1xk + desx
093. r3y = r1yk + desy
094. r3z = r1zk + desz
095.
096. do j = 1, lmn
097. if (j .ne. k) then
098. xi = r3x
099. yi = r3y
100. zi = r3z
101. xx = xi - r1x(j)
102. yy = yi - r1y(j)
103. zz = zi - r1z(j)
104. call mic (xx,yy,zz,xlx,yly,zlz)
105. r = xx*xx + yy*yy + zz*zz
106. if (r .lt. 1d0) go to 1
107. end if
108. end do
109.
110. naccept = naccept + 1
111. nacceptb = nacceptb + 1
112. r1x(k) = r3x
113. r1y(k) = r3y
49
114. r1z(k) = r3z
115.
116. 1 end do
117. end do
118.
119. acceptratio = real(nacceptb)/(nblock*lmn)*100
120. if (acceptratio .gt. 50.0) amp = amp * 1.05
121. if (acceptratio .lt. 50.0) amp = amp * 0.95
122. end do
123.
124. do i=1,lmn
125. write(3,*) r1x(i),r1y(i),r1z(i)
126. end do
127. write(3,*) xlx,yly,zlz
128.
129. write(8,(% Mov acept =,f10.6)) naccept*100./(nstep*lmn)
130.
131. if (istage .eq. 2) then
132.
133. do i = 1, ngrx
134. r = (i-1+0.5d0) * hr
135. dvol = 4*pi/3d0 * ((r+hr/2)**3-(r-hr/2)**3)
136. grr = gr(i)/(dvol*lmn*numgr*(lmn-1)/volume)
137. write(4,*) r,grr
138. end do
139.
140. end if
141.
142. end do
143.
144. end
145. !
146. ! Periodic boundary conditions
147. !
148. subroutine pbc (r1x,r1y,r1z,xlx,yly,zlz)
149.
150. implicit real*8 (a-h,o-z)
151.
152. r1x=r1x-xlx*dnint(r1x/xlx)
153. r1y=r1y-yly*dnint(r1y/yly)
154. r1z=r1z-zlz*dnint(r1z/zlz)
155.
156. end
157. !
158. ! Minimum image convention
159. !
160. subroutine mic (xx,yy,zz,xlx,yly,zlz)
50
161.
162. implicit real*8(a-h,o-z)
163.
164. xx=xx-xlx*dnint(xx/xlx)
165. yy=yy-yly*dnint(yy/yly)
166. zz=zz-zlz*dnint(zz/zlz)
167. end
168. !
169. ! Generation of fcc lattice
170. !
171. subroutine fcc (lmn, dens, r1x, r1y, r1z, xlx, yly, zlz)
172.
173. implicit real*8 (a-h, o-z)
174. parameter (num = 4000)
175.
176. dimension r1x(num), r1y(num), r1z(num)
177. dimension sx(4), sy(4), sz(4)
178.
179. data sx /0d0, 0.5d0, 0.5d0, 0d0/
180. data sy /0d0, 0.5d0, 0d0, 0.5d0/
181. data sz /0d0, 0d0, 0.5d0, 0.5d0/
182.
183. data sh /0.01d0/
184.
185. a = (4d0/dens)**(1d0/3d0)
186. n = (lmn/4)**(1d0/3d0) + 0.001
187.
188. xlx = n*a
189. yly = n*a
190. zlz = n*a
191.
192. m = 1
193. do i = 1, n
194. do j = 1, n
195. do k = 1, n
196. do l = 1, 4
197. r1x(m) = (i - 1 + sx(l) + sh) * a - xlx/2
198. r1y(m) = (j - 1 + sy(l) + sh) * a - yly/2
199. r1z(m) = (k - 1 + sz(l) + sh) * a - zlz/2
200. m=m+1
201. end do
202. end do
203. end do
204. end do
205.
206. end
207. !
51
208. ! Histogram for g(r)
209. !
210. subroutine gdr(lmn,r1x,r1y,r1z,xlx,yly,zlz)
211.
212. implicit real*8(A-H,O-Z)
213. parameter (num=4000, ngrx=1000)
214.
215. dimension r1x(num), r1y(num), r1z(num)
216. common gr(ngrx)
217. data hr /0.05d0/
218.
219. do i = 1, lmn-1
220. do j = i + 1,lmn
221. xx = r1x(i) - r1x(j)
222. yy = r1y(i) - r1y(j)
223. zz = r1z(i) - r1z(j)
224. call mic (xx,yy,zz,xlx,yly,zlz)
225. r = xx*xx + yy*yy + zz*zz
226. rr = dsqrt(r)
227. if (rr .lt. xlx/2) then
228. k = rr/hr + 1
229. gr(k) = gr(k) + 2
230. end if
231. end do
232. end do
233.
234. end
The necessary input parameters are read in from le hs.par. Here, ini is a variable
such that a fcc lattice is built up when ini=0; otherwise the initial congurations
are read in from le hs.inp. The number of particles is lmn (this must be of the
form 4n
3
, with n an integer, to conform with fcc symmetry). The number of steps
in the equilibration stage is nstepeq, and the number of averaging steps are in
nstep0. These steps are divided into blocks of size nblock0 (averages are accu-
mulated every nblock0 steps). In lines 028035 the initial conguration is set up.
The two stages, equilibration and averaging, are controlled by the doloop in line
047. Periodic boundary conditions are applied at the beginning (lines 057059),
and accumulators set to zero (lines 061066). Blocks are opened in line 068, and
subblocks in line 075. Trial moves are applied sequentially (line 082), by randomly
choosing a displacement vector with amplitude amp along all three directions, in
lines 083094 (this variable is set in le hs.par). Lines 096108 perform the overlap
test. If positive, a new particle move is tested (lines 106 and 116); otherwise the
coordinates of the particle tested are updated (lines 112114). In lines 119121 the
amplitude amp for the maximum displacement is changed according to whether the
acceptance ratio in the last block of steps was dierent from 50%. In line 122 the
loop on steps is closed, and particle positions are dumped to le hs.out. In case
we are in the averaging stage, the histogram for the radial distribution function is
52
normalised and written on le (lines 133138). After this the necessary subroutines
are dened: pbc (periodic boundary conditions), mic (minimum image convention),
fcc (generation of fcc lattice), and gdr (calculation of histogram for g(r)).
The main result from this code is the calculation of the radial distribution func-
tion. Fig. 22 collects six radial distribution functions obtained, with this code,
for six dierent densities. As the density is increased the oscillatory structure of
the function becomes more apparent, and between
3
= 0.9 and 1.0 the function
develops peaks located below r = , which is a clear signature of angular order and
hence of crystalline order.
53
0
1
2
3
4
5
6
1 1.5 2 2.5 3
0
1
2
3
4
5
6
1 1.5 2 2.5 3
0
1
2
3
4
5
6
1 1.5 2 2.5 3
0
1
2
3
4
5
6
1 1.5 2 2.5 3
0
1
2
3
4
5
6
1 1.5 2 2.5 3
0
1
2
3
4
5
6
1 1.5 2 2.5 3
r r

3
=1.1

3
=1.0

3
=0.9
3
=0.6

3
=0.7

3
=0.8
Figure 22: Radial distribution functions g(r) of the hardsphere system for six dierent reduced
densities
3
(indicated in each panel).
54
3. THE MOLECULAR DYNAMICS METHOD
3.a Introduction
In this method the equations of motion of a system of N interacting particles are solved
numerically. Depending on the imposed constrainst, 3N or less equations will have to be
solved, one for each degree of freedom. Suppose there are 3N of such degrees, taken here
to be the positions r
i
, i = 1, ..., N. The central quantity is the potential energy,
U = U(r
i
) = U(r
1
, r
2
, ..., r
N
).
The equations of motion are:

dr
i
dt
= v
i
m
dv
i
dt
=
i
U(r
j
)
where v
i
is the velocity of the ith particle, and m is the mass of the particles. Since the
dynamics is energyconserving, the evolution of the system proceeds as if the state point
moved on the hypersurface dened by H(r
i
, v
i
) =const. in the 3Ndimensional phase
space.
3.b Methods of integration
The numerical solution of the equations of motion is performed with an algorithm based
on nite dierences. Time is discretised using a time interval h. Knowledge of the posi-
tions and velocities of all particles at time t allows calculation of the forces at time t, F
i
,
on all particles, and the problem is how to obtain the positions and velocities at a later
time t + h. The simplest algorithm is the Verlet algorithm, which we now derive.
Verlet algorithm. We write the following Taylor expansions:
r
i
(t + h) = r
i
(t) + hv
i
(t) +
h
2
2m
F
i
(t) +
h
3
6
v
i
(t) + ...
r
i
(t h) = r
i
(t) hv
i
(t) +
h
2
2m
F
i
(t)
h
3
6
v
i
(t) + ...
Adding, neglecting terms of order O(h
4
), and reshuing,
r
i
(t + h) = 2r
i
(t) r
i
(t h) +
h
2
m
F
i
(t)
This recurrence formula allows to calculate the new position if one knows the po-
sitions r
i
(t) and r
i
(t h). It is an O(h
3
) algorithm [i.e. the new positions contain
errors of O(h
4
)].
This version of the algorithm is the Newtonian version. The velocities are not
required to calculate the new positions, but of course they are needed in case we
55
wanted to calculate the kinetic energy, related to an important quantity such as the
temperature. Velocities can be approximated with the expansion
r
i
(t + h) = r
i
(t h) + 2hv
i
(t) + O(h
3
)
and from here
v
i
(t) =
r
i
(t + h) r
i
(t h)
2h
which contains errors of O(h
3
). The kinetic energy, at time t, can be estimated from
E
c
(t) =
N

i=1
1
2
m[v
i
(t)[
2
and the temperature, using the equipartition theorem, is

1
2
mv
2

=
kT
2
T =

i=1
m[v
i
(t)[
2
kN
f

where denotes an arbitrary degree of freedom, and N

f
is the number of degrees
of freedom of the system (N
f
= 3N if no contraints are imposed; in simulations
it is quite common to x the position of the centre of mass, so we would have
N
f
= 3N 3).
leap-frog version. This version is numerically more stable than Verlets. One rst
denes:
v
i

t
h
2

=
r
i
(t) r
i
(t h)
h
, v
i

t +
h
2

=
r
i
(t + h) r
i
(t)
h
.
Positions are then updated with
r
i
(t + h) = r
i
(t) + hv
i

t +
h
2

Using the Verlet algorithm,

r
i
(t + h) r
i
(t) = r
i
(t) r
i
(t h) +
h
2
m
F
i
(t)
so that
v
i

t +
h
2

= v
i

t
h
2

+
h
m
F
i
(t).
From here the Hamiltonian version of the Verlet algorithm follows:

r
i
(t + h) = r
i
(t) + hv
i

t +
h
2

,
v
i

t +
h
2

= v
i

t
h
2

+
h
m
F
i
(t).
56
Note that the positions and velocities are outofphase by a time interval h/2.
Velocities at t can be calculated from the expression
v
i
(t) =
v
i

t +
h
2

+v
i

t
h
2

2
. (11)
Another version. The problem with the latter version is that r
i
(t) and v
i
(t) are
not obtained at the same order in h at the same time t. A possible modication
that avoids this is

r
i
(t + h) = r
i
(t) + hv
i

t +
h
2

,
v
i
(t + h) = v
i
(t) +
h
2m
[F
i
(t + h) +F
i
(t)] .
(12)
This numerical scheme is completely equivalent to Verlets. To check it out, note
that
r
i
(t + 2h) = r
i
(t + h) +v
i
(t + h)h +
h
2
2m
F
i
(t + h),
r
i
(t) = r
i
(t + h) v
i
(t + h)h
h
2
2m
F
i
(t).
Adding,
r
i
(t + 2h) +r
i
(t) = 2r
i
(t + h) + h [v
i
(t + h) v
i
(t)] +
h
2
2m
[F
i
(t + h) F
i
(t)] .
Using the second of Eqns. (12), we obtain
r
i
(t + h) = 2r
i
(t) r
i
(t h) +
h
2
m
F
i
(t),
which is precisely Verlets algorithm.
3.c Stability of trajectories
Systems with many degrees of freedom are prone to being unstable in the following sense.
Take two possible trajectories in phase phase that start from very close initial conditions.
Let us calculate the distance (t) between these two trajectories as if phase space were
Euclidean. Expressing the distance in terms of an exponential function (Fig. 23) as
(t) e
t
, we say that the system is
STABLE if < 0 or decreases more slowly than exponentially
UNSTABLE if > 0; now the system is said to be chaotic
The coecient is known as Lyapunov exponent. The criterion for the separation be-
tween trajectories based on an exponential comes naturally from a linear analysis of the
eect a perturbation of amplitude would have on a given trajectory; the linear equation
57
~ e
t
P
H
A
S
E
S
P
A
C
E
Figure 23: Exponential increase of the separation between two trajectories that started very
close in phase space.
would be

, which gives an exponential as a solution.
The Lyapunov instability occurs when > 0; there are two reasons why it is impor-
tant:
gives a limiting upper time beyond which an accurate and trustable numerical so-
lution can be found
in order to reach a high precision when computing a trajectory after a time t, we
would need to know the initial condition with an unreasonably accuracy, since the
number of digits grows linearly with time: if is the number of digits,
e
t
10

t
log 10
It turns out that systems with many degrees of freedom, such as those encountered
in condensedmatter physics, are intrinsically unstable, and hence intrinsically chaotic.
Which are the practical consequences of all this? Obviously use of highprecision numer-
ical algorithms to integrate the equations of motion with high accuracy is not only costly
from the computational point of view but also useless in practical terms.
The basic properties that a numerical integration scheme must meet in order to be useful
in condensedmatter physics are:
1. It should be timereversible
Reversibility is a basic property of the microscopic equations of motion:

dr
i
dt
= v
i
m
dv
i
dt
=
i
U
58
If we make the transformation
t t
v
i
v
i
(i.e. the arrow of time is reversed and, at the same time, the sign of the velocities
is reversed, the equations remain unchanged. In Newtonian language,
m
d
2
r
i
dt
2
=
i
U
reversibility with respect to the transformation t t is also obvious. It is easy
to see that the Verlet algorithm respects the property of invariance with respect to
the operation h h:
r
i
(t + h) = 2r
i
(t) r
i
(t h) +
h
2
m
F
i
(t)
r
i
(t h) = 2r
i
(t) r
i
(t + h) +
h
2
m
F
i
(t)
If terms are rearranged in the second equation we obtain the rst. The built-in
irreversibility of some integration algorithms induces an intrinsic energy dissipation
and nonenergyconserving dynamics. This problem may give rise to unwanted and
sometimes even spurious results in the simulations. But there is still the problem
of the irreversibility introduced by the roundo errors in the computer due to use
of oatingpoint arithmetics, which necessarily involves a nite number of decimal
digits. These eects are usually negligible.
2. It should be symplectic
This means the following. Let f(q, p, t) be the probability distribution of the system
so that
f(q, p, t)dqdp
is the number of points (congurations) in phase space contained in a dierential
volume dqdp centred at (q, p) at time t. Each point propagates in time according
to the dynamical equations (e.g. Hamiltons equations). The Liouville theorem
says that the probability distribution function f(q, p, t) evolves in time like an in-
compressible uid, which in mathematical terms means that

f = 0, where the dot
implies total time derivative. Fig. 24 represents pictorially this ow. This property
can be shown to be contained in the dynamical equations, and of course has to be
respected by the numerical integration scheme used.
A symplectic numerical scheme respects Liouvilles theorem, i.e. it conservs phase
space volume. If the numerical algorithm is written in matrix form, this means that
the Jacobian transformation that takes the system from time t to time t +h has to
be equal to unity:

r(t + h)
v(t + h)

= M

r(t)
v(t)

, Jac(M) = 1
59
P
H
A
S
E
S
P
A
C
E
Figure 24: Liouville theorem: volume of phase space is a constant in time and therefore behaves
as an incompressible uid.
where M is the dynamical matrix associated with the particular numerical algo-
rithm. For instance, the Verlet algorithm, written in Hamiltonian form,

r
i
(t + h) = r
i
(t) + hv
i

t +
h
2

v
i

t +
h
2

= v
i

t
h
2

+
h
m
F
i
(t)
can be written, with the help of the vectors
(x
1
, y
1
, z
1
, x
2
, y
2
, z
2
, ..., x
N
, y
N
, z
N
, v
x
1
, v
y
1
, v
z
1
, v
x
2
, v
y
2
, v
z
2
, ..., v
x
N
, v
y
N
, v
z
N
)
in terms of the matrix
M =

1 h.1
0 1

,
the determinant of which is obviously equal to unity.
3.d. The LennardJones potential: practical implementation
The LennardJones potential,
LJ
(r), is a pair potential that depends on the distance
r between two particles. Its functional form is

LJ
(r) = 4

12

and it contains two parameters: , an energy parameter which is equal to the depth of the
potential well (see Fig. 25), and a length parameter, , which is the root of the potential
and can be taken roughly as the diameter of the particles. The potential has a repulsive
part at short distances, which accounts for the repulsion felt by two spherically-symmetric
60
atoms at close distance due to the overlap of electronic clouds and to the Pauli exclusion
principle. At long distances the potential is attractive, due to van der Waals interaction.
At intermediate distances, there is a potential well.
-1
0
1
2
0 1 2
3
r

(
r
)
Figure 25: The LennardJones pair potential.
We will illustrate the Molecular Dynamics technique using this model potential. The
interest in the LennardJones potential is due to the fact that a system of particles inter-
acting via this pair potential presents a realistic phase diagram, containing liquidvapour
and uidsolid transitions. The rst is a due to the presence of attractions in the poten-
tial. Note that the hardsphere model does not contain this feature.
The central problem of Molecular Dynamics is the calculation of forces. In our case
the forces can be calculated analytically by dierentiation of the potential. Let N be the
number of particles in the system. The potential energy of the system is:
U(r
1
, r
2
, ..., r
N
) =
1
2
N

i=1
N

j=1
j = i

LJ
([r
j
r
i
[) =
N

i=1
N

j<i

LJ
([r
j
r
i
[)
The rst version includes a 1/2 prefactor in order not to count the same pair of particles
twice. The force on particle i is:
F
i
=
r
i
U =

j=i

r
i

LJ
([r
j
r
i
[)

j=i
F
ij
Note that F
ij
=
r
i

LJ
([r
j
r
i
[) is the contribution of particle j to the force on particle
i, and that, by Newtons third law, F
ij
= F
ji
; this property helps save valuable computer
time. Now:

r
i

LJ
([r
j
r
i
[) =

LJ
([r
j
r
i
[)
r
i
[r
j
r
i
[ =

LJ
([r
j
r
i
[)
r
j
r
i
[r
j
r
i
[
61
0
0.5
1
1.5
2
2.5
3
-7
-6.5
-6
-5.5
-5
-4.5
-5
-4.5
-4
-3.5
-3
-2.5
-3.37
-3.36
-3.35
0 50 100 150 200 250 300
MD steps
E
tot
E
tot
E
pot
E
kin
Figure 26: Time evolution of (from top to bottom): kinetic, potential, total energy (all shown
with the same vertical scale) and total energy (at a much lower scale) from a MD simulation of a
LennardJones liquid at reduced density
3
= 0.8 and reduced (initial) temperature kT/ = 2.
Since

LJ
(r) =
48

13

1
2

we nally have
F
i
=
48

j=i

r
ij

14

1
2

r
ij

(r
j
r
i
) .
Once the forces are known, the dynamics can be calculated by means of the numerical
integration algorithm.
Quantities can be calculated as time averages (instead as conguration averages, like
in the MC method). For example, the pressure can be calculated from the virial theorem
which, for pairwise additive forces, reads:
p = kT
1
3V

j>i
r
ij
F
ij

t
,
where the brackets indicate time average. The dot product is:

j>i
r
ij
F
ij
=

LJ
([r
j
r
i
[)
r
ij
r
ij
[r
j
r
i
[
=
48r
ij

r
ij

13

r
ij

62
The temperature is computed, as discussed before, with
T =

i=1
m[v
i
(t)[
2
kN
f

t
,
and the potential energy is
U =

i=1
N

j<i

LJ
([r
j
r
i
[)

t
.
The total energy,
E =

i=1
1
2
m[v
i
(t)[
2

t
+

i=1
N

j<i

LJ
([r
j
(t) r
i
(t)[)

t
,
is useful to check that the MD code is free from errors, since it has to be constant during
the simulation. If the leapfrog algorithm is used, on has to be careful to use the velocities
at time t, i.e. at the same time as the positions, using Eqn. (11).
Let us now discuss the question of units. Natural units of energy and length are and ,
respectively, and from these a LJ time scale is = (m
2
/)
1/2
. Therefore dimensionless
positions, velocities and forces are:
r

i
=
r
i

, v

i
=
v
i
/
= v
i

1/2
, F

i
=
F
i
/
=
F
i

.
The equation that updates positions in the leapfrog algorithm becomes:
r

i
(t + h) = r

i
(t) + h

1/2
v

t +
h
2

i
(t + h) = r

i
(t) + h

t +
h
2

,
where h

= h/. The equation for velocities is:

v

t +
h
2

1/2
= v

t
h
2

1/2
+
h

m
2

1/2

i
(t)
so that
v

t +
h
2

= v

t
h
2

+ h

i
(t).
The dimensionless temperature is T

= kT/:
T

=
k

1
kN
f

i=1
m

1/2
v

i
(t)

t
=
1
N
f

i=1
[v

i
(t)[
2

t
.
Note that, since all particles are equivalent,
T

=
1
3N
N

[v

(t)[
2

[v

(t)[
2

= 3T

(13)
63
0
0.5
1
1.5
2
2.5
3
0 50 100 150 200 250 300
MD steps
k
T
/

(
r
e
d
u
c
e
d

t
e
m
p
e
r
a
t
u
r
e
)
Figure 27: Time evolution of temperature for the same experiment as in Fig. 26.
for a generic particle (here we took N
f
= 3N). Dimensionless kinetic and potential
energies are E

c
= E
c
/,
E

c
=
1

1
2
m

i=1

1/2
v

i
(t)

t
=
1
2

i=1
[v

i
(t)[
2

t
,
and U

= U/. The dimensionless density is

=
3
, and the dimensionless pressure
p

= p
3
/ is calculated from the virial theorem:
p

3
=

1
3V

j>i

ij

ij

t
,
so that
p

1
3V

j>i
r

ij
F

ij

t
.
FORTRAN code
Only the main body of code is given, since all the necessary subroutines (implementation
of periodic boundary conditions, minimum image convention, calculation of histogram
for g(r), and initial conguration for positions) are the same as those used in the Monte
Carlo code for hard spheres.
001. !
002. ! P10: MD for Lennard-Jones particles
003. !
64
004. implicit real*8(a-h,o-z)
005. parameter (n=256, ngrx=1000)
006. dimension x(n),y(n),z(n)
007. dimension vx(n),vy(n),vz(n)
008. dimension fx(n),fy(n),fz(n)
009. dimension ig(1000)
010. real*4 r
011. common gr(ngrx)
012.
013. rho=0.8
014. T=2
015. npasos=2000
016.
017. dum = 17367d0
018. pi = 4d0 * datan(1d0)
019. call fcc (n, rho, x, y, z, aL, aL, aL)
020. do i=1,n
021. vr = dsqrt(3*T)
022. call ggub(dum,r)
023. cost = 2*r-1
024. sint = dsqrt(1-cost**2)
025. call ggub(dum,r)
026. fi = r*2*pi
027. vx(i) = vr*sint*dcos(fi)
028. vy(i) = vr*sint*dsin(fi)
029. vz(i) = vr*cost
030. end do
031.
032. aL2 = aL/2d0
033. cut2 = (2.5d0)**2
034. cutr2 = (aL/2)**2
035. ngr = 100
036. hgr = aL2/ngr
037. dt = 0.01d0
038.
039. ec = 0
040. u = 0
041. ap = 0
042.
043. do i = 1,1000
044. ig(i) = 0
045. end do
046.
047. do k = 1,npasos
048.
049. do i = 1, n
050. call pbc (x,y,z,aL,aL,aL)
65
051. end do
052. do i = 1, n
053. fx(i) = 0
054. fy(i) = 0
055. fz(i) = 0
056. end do
057. epot=0
058. do i = 1, n-1
059. xi = x(i)
060. yi = y(i)
061. zi = z(i)
062. do j = i+1, n
063. xx = xi-x(j)
064. yy = yi-y(j)
065. zz = zi-z(j)
066. call mic (xx,yy,zz,aL,aL,aL)
067. r2 = xx**2+yy**2+zz**2
068. if (r2 .lt. cut2) then
069. r1 = 1/r2
070. r6 = r1**3
071. pot=4*r6*(r6-1)
072. u = u+pot
073. epot=epot+pot
074. rr = 48*r6*r1*(r6-0.5d0)
075. fxx = rr*xx
076. fyy = rr*yy
077. fzz = rr*zz
078. ap = ap+rr*r2
079. fx(i) = fx(i)+fxx
080. fy(i) = fy(i)+fyy
081. fz(i) = fz(i)+fzz
082. fx(j) = fx(j)-fxx
083. fy(j) = fy(j)-fyy
084. fz(j) = fz(j)-fzz
085. end if
086. end do
087. end do
088.
089. ekin=0
090. do i=1,n
091. vxi = vx(i)+dt*fx(i)
092. vyi = vy(i)+dt*fy(i)
093. vzi = vz(i)+dt*fz(i)
094. vxx = 0.5d0*(vxi+vx(i))
095. vyy = 0.5d0*(vyi+vy(i))
096. vzz = 0.5d0*(vzi+vz(i))
097. en = vxx**2+vyy**2+vzz**2
66
098. ekin = ekin+en
099. ec = ec+en
100. vx(i) = vxi
101. vy(i) = vyi
102. vz(i) = vzi
103. x(i) = x(i)+dt*vx(i)
104. y(i) = y(i)+dt*vy(i)
105. z(i) = z(i)+dt*vz(i)
106. end do
107.
108. call gdr(x,y,z,aL,aL,aL)
109.
110. end do
111.
112. temp =ec/(3*n*npasos)
113. u = u/(n*npasos)
114. ap = rho*temp+ap/(3*aL**3*npasos)
115.
116. write(*,*) temperature=,temp
117. write(*,*) energy=,u
118. write(*,*) pressure=,ap
119.
120. do k=1,ngr
121. r=(k-1)*hgr+hgr/2d0
122. vol=4*pi/3*((r+hgr/2d0)**3-(r-hgr/2d0)**3)
123. gdr=gr(k)/(n*(n-1)/aL**3*npasos*vol)
124. write(1,*) r,gdr
125. end do
126.
127. end
The code is quite similar to that of a Monte Carlo simulation, except that one has to
take account of the velocities and update coordinates and velocities according to the
Verlet algorithm (in the leap-frog version). Coordinates are initialised in line 019, using
the same subroutine as in the Monte Carlo code for hard spheres (which generates a
fcc lattice). Here we have to generate values for the initial velocities. This is done by
assigning random directions for the velocity vectors with a constant modulus obtained
from the initial temperature [vr, line 021; this relation is obtained from Eqn. (13)], in
lines 020-030. The time interval h, called dt in the code, is set to 0.01 in line 037. The
histogram for the radial distribution function is stored in variable ig. In line 047 the loop
over MD steps is opened. The forces on the particles, the potential energy and the virial
are computed in lines 058-087. Particle coordinates and velocities are updated in lines
090-106, and then the histogram for g(r) is calculated for the current MD step. Finally,
the average temperature temp, total energy u and pressure ap (this is obtained from the
virial) are calculated and printed out, and the radial distribution function is normalised
and dumped to le fort.1.
67
Let us discuss some results generated by this code. In Fig. 26 kinetic, potential and
total energies are represented as a function of MD steps, i.e. simulation time (with a
time step of h = 0.01 in reduced units ). This corresponds to a simulation at uid
conditions, but starting with a crystalline conguration with random velocity vectors of
the same moduli. We can see how both kinetic and potential contributions uctuate quite
substantially, while the total energy remains constant (except for some initial time where
the system readjusts itself due to the quite atypical conguration that we started from).
Note that the total energy does uctuate, but with a much lower amplitude than the
kinetic and potential contributions separately (see bottom panel). We can conclude that
the Verlet algorithm is energyconserving and works well.
Fig. 27 show the evolution of the temperature; this is the same simulation as before,
and the time behaviour of the temperature follows that of the kinetic energy (since they
are proportional). Note that the initial temperature of kT/ = 2 decreases down to an
average of 1.12 along the run, as a result of the internal readjustments of the dierent
contributions to the energy.
In Fig. 28 the uctuating part of the pressure (without the idealgas contribution kT),
i.e. the virial, is plotted against MD time. The horizontal line is the average value over
2000 MD steps. We can see that the virial is a highly uctuating quantity; therefore, in
order to have a reliable value for the pressure, suciently long simulations have to be run.
0
0.5
1
1.5
2
2.5
0 50 100 150 200 250 300
MD steps
r
e
d
u
c
e
d

v
i
r
i
a
l
Figure 28: Time evolution of the virial part of the pressure for the same experiment as in Fig.
26. The horizontal line is the average value over 2000 MD steps.
Finally, in Fig. 29, the radial distribution function is represented for the same conditions;
in this case an initial period of 1000 MD steps for equilibration has been run, to give the
system time to reach typical (uid, i.e. globally disordered) congurations. The average
68
time interval was 2000 MD steps. It has the usual features of the radial distribution
function of a uid phase: a rst pronounced peak and oscillatory damped behaviour at
larger distances.
0
0.5
1
1.5
2
2.5
0 0.5 1 1.5 2 2.5 3 3.5
g
(
r
)
r/
Figure 29: Radial distribution function for a LennardJones uid at conditions
3
= 0.8 and
kT/ = 1.12.
69
LIST OF COMPUTER CODES
P1. Random number generator (page ???)
P2. Metropolis et al. algorithm for f(x) = c exp (x

4)
P3. Random walk
P4. Algorithm for the SARW model
P5. Clustering algorithm (in the course web page)
P6. 2D-Ising model Monte Carlo with importance sampling
P7. Evaluation of B3 by Monte Carlo integration
P8. Free volume for the fcc lattice
P9. MC simulation of hard spheres
P10. MD for LennardJones particles
70