4 Monte Carlo Simulation
The first simulation technique we will study is the Monte Carlo method. The primary “flavor” of MC most appropriate for this introductory level course is the original Metropolis method, explained in Sec. 4.1.
We will also use MC to explain basic technical aspects of molecular simulation code in Sec. 4.3. These aspects are not all restricted to MC, but following the text, we will introduce and discuss these technical aspects here. These include (1) periodic boundary conditions, (2) energy evaluation, (3) representation of data, among others.
Finally, four case studies will be presented and investigated. The first is the Ising magnet (Sec. 4.2). The second considers hard-disks confined in a circle (Sec. 4.4), and the third is hard-disk-dumbbells (Sec. 4.5). The fourth is MC to explore the equation of state of a model liquid known as the Lennard-Jones fluid (Sec. 4.6).
4.1 The Metropolis Monte Carlo Method
Monte Carlo is a type of numerical integration. Consider the following integral: \begin{equation} I = \int _a^b f\left (x\right ) dx \end{equation} Now imagine that we have a second function, \(\rho (x)\), which is positive in the interval \([a,b]\). We can also express \(I\) as \begin{equation} I = \int _a^b \frac {f\left (x\right )}{\rho \left (x\right )}\rho \left (x\right ) dx \end{equation} If we think of \(\rho \left (x\right )\) as a probability density, then what we have just expressed is the expectation value of the quantity \(f/\rho \) on \(\rho \) in the interval \([a,b]\): \begin{equation} I = \left <\frac {f\left (x\right )}{\rho \left (x\right )}\right >_\rho \end{equation} This implies that we can approximate \(I\) by picking \(M\) values \(\left \{x_i\right \}\) randomly out of the probability distribution \(\rho (x)\) and computing the following sum: \begin{equation} I \approx \frac {1}{M} \sum _{i=1}^{M} \frac {f\left (x_i\right )}{\rho \left (x_i\right )} \end{equation} Note that this approximates the mean of \(f/\rho \) as long as we pick a large enough number of random numbers (\(M\) is large enough) such that we “densely” cover the interval \([a,b]\). If \(\rho \left (x\right )\) is uniform on \([a,b]\), \begin{equation} \rho \left (x\right ) = \frac {1}{b-a},\ \ \ a < x < b \end{equation} and therefore, \begin{equation} \label {eq:unif_avg} I \approx \frac {a-b}{M} \sum _{i=1}^{M} f\left (x_i\right ) \end{equation}
The next question is, how good an approximation is this, compared with other one-dimensional numerical integration techniques, such as Simpson’s rule and quadrature? A better phrasing of this question is, how expensive is this technique for a given level of accuracy, compared to traditional techniques? Consider this means to compute \(\pi \): \begin{equation} \label {eq:get_pi} I = \int _0^1 \left (1-x^2\right )^{1/2}dx = \frac {\pi }{4} \end{equation} Allen and Tildesley [2] mention that, in order use Eq. 47 to compute \(\pi \) to an accuracy of one part in 10\(^6\) requires \(M\) = 10\(^7\) random values of \(x_i\), whereas Simpson’s rule requires three orders of magnitude fewer points to discretize the interval to obtain an accuracy of one part in 10\(^7\) (Fig. 2). So the answer is, integral estimation using Monte Carlo estimation with uniform random variates is expensive.
But, the situation changes radically when the dimensionality of the integral is large, as is the case for an ensemble average. For example, for a system of 100 particles comprising 300 coordinates, the configurational average \(\left <G\right >\) (Eq. 39) could be discretized using Simpson’s rule. If we did that, requesting only a modest 10 discrete points per axis in configurational space, we would need to evaluate the integrand \(Ge^{-\beta {U}}\) 10\(^{300}\) times. This is an almost unimaginably large number. Using a direct numerical technique to compute statistical mechanical averages is simply out of the question.
We therefore return to the idea of evaluating the integrand at a discrete set of points selected randomly from a distribution. Here we call upon the idea of importance sampling. Let us try to use whatever we know ahead of time about the integrand in picking our random distribution, \(\rho \), such that we minimize the number of points (i.e., the expense) necessary to give an estimate of \(\left <G\right >\) to a given level of accuracy.
Now, clearly the states that contribute the most to the integrals we wish to evaluate by configurational averaging are those states with large Boltzmann factors; that is, those states for which \(\rho _{NVT}\) is large. It stands to reason that if we randomly select points from \(\rho _{NVT}\), we will do a pretty good job approximating the integral. So what we end up computing is the “average of \(G\rho _{NVT}\) over \(\rho _{NVT}\)”: \begin{equation} \left <G\rho _{NVT}/\rho _{NVT}\right > \approx \left <G\right >, \end{equation} which should give an excellent approximation for \(\left <G\right >\). The idea of using \(\rho _{NVT}\) as the sampling distribution is due to Metropolis et al. [3]. This makes the real work in computing \(\left <G\right >\) generating states that randomly sample \(\rho _{NVT}\).
Metropolis et al. [3] showed that an efficient way to do this involves generating a Markov chain of states which is constructed such that its limiting distribution is \(\rho _{NVT}\). A Markov chain is just a sequence of trials, where (i) each trial outcome is a member of a finite set and (ii) every trial outcome depends only on the outcome that immediately precedes it. By “limiting distribution,” we mean that the trial acceptance probabilities are tuned such that the probability of observing the Markov chain atop a particular state is defined by some equilibrium probability distribution, \(\rho \). For the following discussion, it will be convenient to denote a particular state \(n\) using \(\Gamma _n\), instead of \(\nu \).
A trial is some perturbation (usually small) of the coordinates specifying a state. For example, in an Ising system, this might mean flipping a randomly selected spin. In a system of particles in continuous space, it might mean displacing a randomly selected particle by a small amount \(\delta r\) in a randomly chosen direction \((\theta ,\phi )\). There can be a large variety of such “trial moves” for any particular system.
The probability that a trial move results in a successful transition from state \(n\) to \(m\) is denoted \(\pi _{nm}\) and \(\bf \pi \) is called the “transition matrix.” It must be specified ahead of time to execute a traditional Markov chain. Since the probability that a trial results in a successful transition to any state, the rows of \(\bf \pi \) add to unity: \begin{equation} \sum _i \pi _{ni} = 1 \end{equation} With this specification, we term \(\bf \pi \) a “stochastic” matrix. Furthermore, for an equilibrium ensemble of states in state space, we require that transitions from state to state do not alter state weights as determined by the limiting distribution. So the weight of state \(n\): \begin{equation} \rho _n \equiv \rho _{NVT}\left (\Gamma _n\right ) \end{equation} must be the result of transitions from all other states to state \(n\): \begin{equation} \label {eq:flow} \rho _n = \sum _m \rho _m \pi _{mn}. \end{equation} For all states \(n\), we can write Eq. 52 as a post-op matrix equation: \begin{equation} \label {eq:balance} {\bf \rho \pi } = {\bf \rho } \end{equation} where \(\bf \rho \) is the row vector of all state weights. Eq. 53 constrains our choice of \(\bf \pi \). This means there is still more than one way to specify \(\bf \rho \). Metropolis et al. [3] suggested: \begin{equation} \label {eq:detailed_balance} \rho _m\pi _{mn} = \rho _n\pi _{nm} \end{equation} That is, the probability of transitioning from state \(m\) to \(n\) is exactly equal to the probability of transitioning from state \(n\) to \(m\). This is called the “detailed balance” condition, and it guarantees that the state weights remain static. Observe: \begin{equation} \sum _m \rho _m\pi _{mn} = \sum _m\left (\rho _n\pi _{nm}\right ) = \rho _n\left (\sum _m\pi _{nm}\right ) = \rho _n \end{equation} Detailed balance is, however, overly restrictive; is the only the conceptually simplest way to guarantee that a limiting distribution is obtained. This fact is of little importance in this course, but you may encounter other balance-enforcing conditions in the literature.
Metropolis et al. [3] chose to construct \(\bf \pi \) as \begin{equation} \pi _{nm} = \alpha _{nm} {\rm acc}\left (n\rightarrow m\right ) \end{equation} where \(\alpha \) is the probability that a trial move is attempted, and \(\rm acc\) is the probability that a move is accepted. If the probability of proposing a move from \(n\) to \(m\) is equal to that of proposing a move from \(m\) to \(n\), then \(\alpha _{nm} = \alpha _{mn}\), and the detailed balance condition is written: \begin{equation} \rho _n{\rm acc}\left (n\rightarrow m\right ) = \rho _m{\rm acc}\left (m\rightarrow n\right ) \end{equation} from which follows \begin{equation} \frac {{\rm acc}\left (n\rightarrow m\right )}{{\rm acc}\left (m\rightarrow n\right )} = \frac {\rho _m}{\rho _n} = \frac {e^{-\beta {U}\left (\Gamma _m\right )}}{e^{-\beta {U}\left (\Gamma _n\right )}} \end{equation} giving \begin{equation} \label {eq:acc} \frac {{\rm acc}\left (n\rightarrow m\right )}{{\rm acc}\left (m\rightarrow n\right )} = \exp \left \{-\beta \left [{U}\left (\Gamma _m\right )-{U}\left (\Gamma _n\right )\right ]\right \} \equiv \exp \left (-\beta \Delta {U}_{nm}\right ) \end{equation} where we have defined the change in potential energy as \begin{equation} \Delta {U}_{nm} = {U}\left (\Gamma _m\right )-{U}\left (\Gamma _n\right ) \end{equation}
There are many choices for \({\rm acc}\left (n \rightarrow m\right )\) that satisfy Eq. 59. The original choice of Metropolis is used most frequently: \begin{equation} {\rm acc}\left (n \rightarrow m\right ) = \left \{\begin {array}{ll} \exp \left (-\beta \Delta {U}_{nm}\right ) & \Delta {U}_{nm} > 0\\ 1 & \Delta {U}_{nm} < 0 \end {array} \right . \end{equation} So, suppose we have some initial configuration \(n\) with potential energy \(U_n\). We make a trial move, temporarily generating a new configuration \(m\). Now we calculate a new energy, \(U_m\). If this energy is lower than the original, (\(U_m < U_n\)) we unconditionally accept the move, and configuration \(m\) becomes the current configuration. If it is greater than the original, (\(U_m > U_n\)) we accept it with a probability consistent with the fact that the states both belong to a canonical ensemble. How does one in practice decide whether to accept the move? One first picks a uniform random variate \(x\) on the interval \([0,1]\). If \(x \le {\rm acc}\left (n \rightarrow m\right )\), the move is accepted.
The next section is devoted to an implementation of the Metropolis Monte Carlo method for a 2D Ising magnet.
4.2 Case Study 1: The 2D Ising Magnet
4.2.1 Introduction
In the introductory lecture, I introduced state-counting using a simple, one-dimensional Ising system. To be precise, this was a particular case known as an Ising magnet of noninteracting spins, because each spin made its own independent contribution to the Hamiltonian. The state of the magnet is specified by specifying each spin variable as either +1 or -1:
Now, we will consider a more interesting Ising system; namely, that of interacting spins on a 2d square lattice. What do we mean by interacting? We mean that the Hamiltonian depends on pairs of spins: \begin{equation} \label {eq:2dising} \mathscr {H} = -J\sum _{\left <ij\right >}s_is_j \end{equation} where the notation \(\left <ij\right >\) denotes that we are summing over unique pairs of nearest neighbors, and, as before a spin \(s_i\) is either +1 or -1. \(J\) is the “coupling constant” and represents our default unit of energy. How many unique pairs of nearest neighbors are there on a lattice of \(N\) spins (assuming periodic boundary conditions)? To answer this question, we must know the coordination, \(z\), of the lattice; that is, how many nearest neighbors does one lattice position have? For a square lattice, \(z\) = 4. Each spin therefore contributes two unique spin pairs to the system, so there are \(Nz/2\) unique nearest neighbor pairs.
Think about this Hamiltonian. It says that energy is minimal when all \(N\) spins have the same alignment, either all up or all down. Imagine a microscopic observable called the magnetization, or average spin orientation, \(M\): \begin{equation} M \equiv \frac {1}{N} \sum _i^N s_i \end{equation} Then, \(M\) = +1 and -1 are two energetically equivalent microstates. We can then expect that if there is any thermal energy in the system which randomly flips spins, the amount of this thermal energy (that is, the temerature) will somehow control the observable magnetization. We “observe” magnetization using an ensemble average: \begin{equation} \label {eq:magnetization} \left <M\right > = \frac {\displaystyle \sum _{s_1=-1}^{+1} \cdots \sum _{s_N=-1}^{+1} \left [\sum _i^N s_i\right ]\exp \left [-\beta \mathscr {H}\left (s_1,\dots ,s_N\right )\right ]} {\displaystyle \sum _{s_1=-1}^{+1} \cdots \sum _{s_N=-1}^{+1} \exp \left [-\beta \mathscr {H}\left (s_1,\dots ,s_N\right )\right ]} \end{equation} Our physical intuition tells us that as \(T\rightarrow \) 0, \(\left |\left <M\right >\right |\rightarrow \) 1, and as \(T\rightarrow \infty \), \(\left <M\right >\rightarrow \) 0. The fascinating thing about an Ising magnet is that there is a finite temperature called the critical temperature, \(T_c\). If we start out with a “hot” system, and cool it to just below \(T_c\), the absolute value of the magnetization spontaneously jumps from 0 to some finite value. In other words, the system undergoes a phase transition from a disordered phase to a partially ordered phase. A Metropolis Monte Carlo simulation can allow us to probe the behavior of an Ising system and learn how the system behaves near criticality. The rest of this case study will be devoted to showing you the inner workings of a C program which simulates the Ising magnet using Metropolis MC, as a first implementation of this technique. In the suggested exercises appearing at the end of this case study, you will modify this code slightly to compute averages values of certain observables.
But first, I recommend you check out this Java implementation of a Monte Carlo simulation of a 2-dimensional Ising magnet (one of many on the web; google “ising simulation” and you’ll get a nice sample.) This is a fun little Java applet that lets you play with an Ising magnet. You can change the temperature of the simulation: making it cold will “freeze” the system, and making it hot “melts” it. Near the critical temperature, \(T_c\), relatively large regions of mostly-up spins compete with regions of mostly down spin. In one of the exercises, we’ll learn how to measure an observable called the correlation length, which characterizes the size of these domains and is a useful signature of criticality.
4.2.2 A C Code for the 2D Ising Magnet
In this section, we will dissect piece-by-piece a small C program that implements an NVT Metropolis Monte Carlo simulation of a 2D Ising magnet.
If you have not already done so, clone the Abrams-Teaching/instructional-codes repository from Github.
This repository will be updated throughout the term with source code in the originals subdirectory. You
should create a subdirectory called my_work inside this repository and do all editing, compiling, and running in
there. This directory is specifically excluded from git revision control by its inclusion in the .gitignore
file. This way, when I put new codes in the repository, you only have to git pull to download
them.
cd cd cheT580 git clone git@github.com:Abrams-Teaching/instructional-codes.git cd instructional-codes mkdir my_work cd my_work cp ../originals/ising_mc.c . code .
From a terminal command-line inside VSCode, or outside, you can compile ising_mc.c via
gcc -O3 -o ising ising_mc.c -lm -lgsl
and you can then run it at the command line as ./ising. It has a lot of options for controlling the size
of the magnet (the number of spins), the temperature, and other parameters, which I’ll go over
now.
ising_mc.c conducts a canoncial Metropolis Monte Carlo simulation of an Ising magnet of size \(L\times L\) at
temperature \(T\) (both specified by the user at run time on the command line), and it computes both the average
energy per spin \(\left <\epsilon \right >\) and the average spin value, \(\left <s_1\right >\).
Periodic boundaries are employed in calculating the nearest-neighbor interactions. Consider an \(L\times L\) magnet; each row \(i\) indexed from 0 to \(N-1\) has \(L\) columns also indexed from 0 to \(N-1\). If the cell at (5,5) queries its neighbors, they are at (4,5), (6,5), (5,4), and (5,6). However, the cell at (0,5) would have a southern neighbor at (\(N-1\),5) instead of (-1,5), since there is no row indexed -1! Fig. 3 demonstrates this.
An abbreviated listing of the code follows. Some comments in the full, downloadable code have been omitted for space, and I have instead explained each code fragment in accompanying text.
Headers
These are some standard headers we include in most C programs, along with the GNU scientific library.
#include <stdio.h> #include <stdlib.h> #include <math.h> #include <gsl/gsl_rng.h>
The Change in Energy upon a Proposed Spin Flip
dE() is function that accepts as arguments the 2D array of spins, M; the side-length, L, and a position \((i,j)\), and
returns the change in energy (in units of \(J\)) upon flipping spin \((i,j)\), without actually flipping it. All variables are of
integer type, int. The syntax ** M means that M is a pointer to a pointer to an int, signifying that M is a 2D
array. To access element in row i and column j in M, the syntax is M[i][j]. In C, all array indices start at 0, so
M[0][0] refers to the “upper-left” corner of the magnet (assuming row numbers increase going down the
magnet). The inline conditional syntax i?(i-1):(L-1) expands as, “If i is non-zero, return i-1; otherwise,
return L-1.” This implements periodic boundaries when looking to the west of spin M[i][j]. The syntax
(i+1)\%L returns computes \((i+1)\mod L\), which also implements periodic boundaries when looking to the east of spin
M[i][j].
int dE ( int ** M, int L, int i, int j ) {
return -2*(M[i][j])*(M[i?(i-1):(L-1)][j]+M[(i+1)%L][j]+
M[i][j?(j-1):(L-1)]+M[i][(j+1)%L]);
}
Initialize all Spin Values
The function init() visits every spin and assigns it either +1 or -1 randomly.
void init ( int ** M, int L, gsl_rng * r ) {
int i,j;
for (i=0;i<L;i++) {
for (j=0;j<L;j++) {
M[i][j]=2*(int)gsl_rng_uniform_int(r,2)-1;
}
}
}
Main Program: Variable Declarations
In C functions, including main(), all variables must be declared by type before they are used. Often, it makes
sense to initialize them to some default values upon declaration. Here, we make the default magnet 20\(\times \)20 in
size, set the temperature to 1 (in dimensions of \(J/k_B\)), and provide some defaults fo the number of cycles (10\(^6\)) and
sampling interval (10\(^3\)) in cycles (a cycle is one round of \(N\) attempted flips). Variables for holding a
couple of observables are initialized. Finally, the pseudorandom number generator from the GSL is
declared.
int main (int argc, char * argv[]) {
/* System parameters */
int ** M; /* The 2D array of spins */
int L = 20; /* The sidelength of the array */
int N; /* The total number of spins = L*L */
double T = 1.0; /* Dimensionless temperature = (T*k)/J */
/* Run parameters */
int nCycles = 1000000; /* number of MC cycles to run;
one cycle is N consecutive attempted
spin flips */
int fSamp = 1000; /* Interval between taking samples */
/* Computational variables */
int nSamp; /* Number of samples taken */
int de; /* energy change due to flipping a spin */
double b; /* Boltzman factor */
double x; /* random number */
int i,j,a,c; /* loop counters */
/* Observables */
double s=0.0, ssum=0.0; /* average magnetization */
double e=0.0, esum=0.0; /* average energy per spin */
/* Pseudorandom Number Generator */
gsl_rng * r = gsl_rng_alloc(gsl_rng_mt19937);
unsigned long int Seed = 23410981;
Main Program: Argument Parsing
Here, we parse the command line arguments. The user running the program can specify the magnet
side-length \(L\), the temperature (in units of \(J/k_B\)), the number of cycles, the sampling interval, and the seed value for
the pseudorandom number generator. The array argv[] and its count of elements argc appear as parameters
in the definition of main, as is customary in C. The built-in function strcmp() returns 0 if the two
arguments match, so the logical expression strcmp(a,b)! evaluates to TRUE if strings a and b
match. The built-in functions atoi() and atof() convert their string arguments to integers and
floating-points, respectively. The notation argv[++i] means that first the current value of i is
incremented by 1 and then it is used to access the value in argv[]. So, for example, if the string "-L"
is detected at the ith position in the array of command-line arguments, the code immediately
jumps to the next position in the array and tries to intepret that argument as an integer to assign to
L.
for (i=1;i<argc;i++) {
if (!strcmp(argv[i],"-L")) L=atoi(argv[++i]);
else if (!strcmp(argv[i],"-T")) T=atof(argv[++i]);
else if (!strcmp(argv[i],"-nc")) nCycles = atoi(argv[++i]);
else if (!strcmp(argv[i],"-fs")) fSamp = atoi(argv[++i]);
else if (!strcmp(argv[i],"-s")) {
Seed = (unsigned long)atoi(argv[++i]);
}
}
Main Program: Initial Outputs
Next, we echo the command entered back to the terminal, and then output any parsed variables values.
printf("# command: ");
for (i=0;i<argc;i++) printf("%s ",argv[i]);
printf("\n");
printf("# ISING simulation, NVT Metropolis Monte Carlo\n");
printf("# L = %i, T = %.3lf, nCycles %i, fSamp %i, Seed %lu\n",
L,T,nCycles,fSamp,Seed);
Main Program: Initialization
Next, we initialize the pseudorandom number generator r by setting the seed value using Seed.
Then the number of spins \(N\) is computed assuming a square magnet. We then allocate memory
needed to hold the magnet; this is a standard method of allocating a 2D array of integers. Next
we call our init() function, and finally we convert \(T\) to \(1/T\) since multiplication is more efficient than
division.
/* Seed the random number generator */ gsl_rng_set(r,Seed); /* Compute the number of spins */ N=L*L; /* Allocate memory for the system */ M=(int**)malloc(L*sizeof(int*)); for (i=0;i<L;i++) M[i]=(int*)malloc(L*sizeof(int)); /* Generate an initial state */ init(M,L,r); /* For computational efficiency, convert T to reciprocal T */ T=1.0/T;
Main Program: Monte Carlo Loop
And now the loop. Note the outer loop counts cycles, and the inner loops counts flip attempts in one cycle. The
variables i and j are randomly assigned between 0 and \(L\)-1, identifying a random spin. This random
spin is tagged and sent to our dE() function to calculate the change in energy we would expect if
that spin were flipped (+1\(\rightarrow \)-1 or -1\(\rightarrow \)+1). We then calculate the Boltzmann factor in b, and then use
the Metropolis criterion to decide whether to accept or reject this spin flip: choosing a random
variable on [0,1] and accepting the move if the Boltzmann factor is greater than this number. If it is
accepted, we actually peform the flip by multiplying the spin value by -1. Finally, if the current
cycle is one in which we collect a sample, we do so by calling the samp() function. The logical
expression (c%fSamp)! evaluates to TRUE if the cycle counter c is divisible by fSamp. In the call
to samp(), the arguments s and e are passed by reference, signified by the & qualifier. This is
necessary because inside samp() we modify both variables. Each of s and e are added to their
respective tallies, ssum and esum, and the number of samples taken nSamp is incremented by
1.
s = 0.0;
e = 0.0;
nSamp = 0;
for (c=0;c<nCycles;c++) {
/* Make N flip attempts */
for (a=0;a<N;a++) {
/* randomly select a spin */
i=(int)gsl_rng_uniform_int(r,L);
j=(int)gsl_rng_uniform_int(r,L);
/* get the "new" energy as the incremental change due
to flipping spin (i,j) */
de = dE(M,L,i,j);
/* compute the Boltzmann factor; recall T is now
reciprocal temperature */
b = exp(de*T);
/* pick a random number between 0 and 1 */
x = gsl_rng_uniform(r);
/* accept or reject this flip */
if (x<b) { /* accept */
/* flip it */
M[i][j]*=-1;
}
}
/* Sample and accumulate averages */
if (!(c%fSamp)) {
samp(M,L,&s,&e);
fprintf(stdout,"%i %.5le %.5le\n",c,s,e);
fflush(stdout);
ssum+=s;
esum+=e;
nSamp++;
}
}
Main Program: Final Outputs and Program Termination
After the outer loop finishes, we can finish up by reporting the averages \(\left <s_1\right >\) from the tallies of s and
e.
fprintf(stdout,"# The average magnetization is %.5lf\n", ssum/nSamp); fprintf(stdout,"# The average energy per spin is %.5lf\n", esum/nSamp); fprintf(stdout,"# Program ends.\n"); }
4.2.3 Example: Average Energy and Magnetization vs. Temperature
Let’s run the code for an \(L\)=40 lattice for the following values of temperature: 5.0, 4.0, 3.0, 2.0, 1.0. At each temperature, we’ll run six separate simulations with a unique random number generator seed. Fig. 4 shows the average magnetization \(\langle s_1\rangle \) and the average energy per spin \(\langle e\rangle \) vs. temperature, with each simulation contributing a unique point.
What is happening here? Clearly, as the temperature decreases, the average energy per spin approaches -2; this makes sense because the spins will tend to align with their neighbors. However, when we look at the average magnetization, we see that \(\langle s_1\rangle \) is zero at high temperatures but then can seemingly approach either +1 or -1 as the temperature goes down. This is because of symmetry breaking: although either all-up or all-down is equally likely, once it falls toward one it won’t ever climb back out and go to the other. So we see some magnets go to all -1 and others to all +1.
Symmetry-breaking is an important phenomenon in molecular simulations. The impact is has is to prevent exploration of state space because of barriers that are only extremely rarely overcome when resolving state-to-state transitions microscopically. In the low-\(T\) Ising magnet, the all-up and all-down states are equally likely, but once one is committed to, the other will never practically be explored. Why is this significant? Many systems have state spaces like this, where there are two or more high-probability regions separated by large regions of low probability. Sampling based on local moves in state space can almost never overcome such barriers, meaning such simulations are likely never truly ergodic. None of the MC simulations below about \(T\) = 2.2 for an Ising simulation are ergodic!
4.2.4 Suggested Exercises
- 1.
- Modify the code so that when samples are taken in accumulating statistics for \(\left <s_1\right >\) and \(\left <\epsilon \right >\), the current sample values are output to the terminal. You’ll want to find the right place to add the following line: fprintf(stdout,"%i %.5lf %.5lf\(\backslash \)n",c,s,e);
- 2.
- The current version of the code initializes the Ising lattice with random spins. What temperature does this correspond to? Modify the code so that the initial lattice has two well-defined domains, all spin-up for \(i < L/2\) and all spin-down for \(i > L/2\). Re-run at the various temperatures. Do you see any differences?
- 3.
- (Advanced) Modify the code ising.c to compute the quantity \(\left <s_is_j\right > - \left <s_i\right >\left <s_j\right >\) as a function of various distances between spins \(i\) and \(j\).
4.3 Elements of a Continuous-Space MC program
The Ising system serves us well as a simple introduction to the technique of Monte Carlo simulation. Now, we move on to the more advanced case of continuous-space Monte Carlo; that is, Monte Carlo on a system composed of particles whose position and velocity vector components are real numbers. We will first consider the simple case of a “hard-sphere” liquid, and then the more realistic Lennard-Jones liquid. These distinctions have to do with the potential energy function used to compute the potential energy \(U\) of a configuration \({\bf r}^N\).
Regardless of the potential used, all continuous-space MC programs have common elements:
- 1.
- data representation and input/output;
- 2.
- energy calculation;
- 3.
- trial move generation.
4.3.1 Data Representation and Input/Output
The information that must be stored and handled in an MC program are the particle positions. In 3D space, each particle has three components of its position. The simplest way to represent this information in a program is by using parallel arrays; that is, three arrays, dimensioned from 1 to \(N\), one for each dimension. In C, we might declare these arrays for 1,000 particles as
int N = 1000; double rx[N], ry[N], rz[N];
This data can be stored in standard ASCII text files, or in unformatted binary files. In fact, a large part of most molecular simulation is producing and storing configurations which can be processed “offline” (away from the simulation program) to perform analyses. A simple and widely used ASCII configuration file format is called “XYZ”. A code fragment to write a configuration of \(N\) hydrogen atoms in XYZ format appears below:
FILE * fp;
...
fp = fopen("my_config.xyz","w");
fprintf(fp,"%i\n",N);
fprintf(fp,"This is a comment or just a blank line\n");
for (i=0;i<N;i++)
fprintf(fp,"%s %.8le %.8le %.8le\n","H",rx[i],ry[i],rz[i]);
The main feature of XYZ files is that one can specify the element type of each atom; for simple visualization purposes of our toy systems, this can be anything, so we’ll just use H for now.
It’s quite easy to read ASCII data in that was written in XYZ format, e.g., using the fragment above, ignoring the element types (for now):
FILE * fp;
int N;
double * rx, * ry, * rz;
char dummy[4];
...
fp = fopen("my_config.xyz","r");
fscanf(fp,"%i\n",&N);
/* allocate if not done already */
rx=(int*)malloc(N*sizeof(int));
ry=(int*)malloc(N*sizeof(int));
rz=(int*)malloc(N*sizeof(int));
fscanf(fp,"\n"); /* reads the comment line but throws it away */
for (i=0;i<N;i++)
fscanf(fp,"%s %.8le %.8le %.8le\n",dummy,&rx[i],&ry[i],&rz[i]);
4.3.2 Analytical Potentials
The total system potential energy is most often computed by means of analytical potentials. Such a potential energy can be written as a continuous function of the positions of particle centers-of-mass. (Such potentials are often termed “empirical” when applied to atomic systems because they formally neglect quantum mechanics.) In most molecular simulations, the total system potential can be decomposed in the following way: \begin{equation} {U}\left ({\bf r}^N\right ) = U_0 + \sum _i U_1( {\bf r}_i ) + \sum _{i<j} U_2( {\bf r}_i, {\bf r}_j ) + \sum _{i<j<k} U_3( {\bf r}_i, {\bf r}_j, {\bf r}_k ) + \cdots + U_N( {\bf r}_i, {\bf r}_j, \dots , {\bf r}_N ), \end{equation} Here, \(U_0\) is some reference potential. \(U_1\) is a “one-body” potential, \(U_2\) is a “two-body” potential, etc. This generally defines a “many-body” potential, and is the heart of the “many-body” problem of statistical physics. Namely, the Boltzmann factor \(e^{-\beta {U}}\) correlates every position with every other position, and the configurational integral (Eq. 39) is not factorizable into many separate easily evaluated integrals.
Analytical potentials are most easily understood by considering model systems which are decomposable into spherically-symmetric, pairwise interactions. Consider then the following total system potential energy: \begin{equation} {U} = \sum _{i=1}^{N}\sum _{j=i+1}^{N} U_{ij}\left (r_{ij}\right ) \end{equation} The double-sum runs over all unique pairs of particles. The function \(U_{ij}\left (r_{ij}\right )\) is called a pair potential, and it is a function of the scalar distance between particle \(i\) and \(j\).
The simplest pair potential is the “hard-sphere”: \begin{equation} \label {eq:hs} U_{HS}\left (r\right ) = \left \{\begin {array}{ll} \infty & r < \sigma \\ 0 & r > \sigma \end {array}\right . \end{equation} Here, \(\sigma \) is an arbitrary interaction range that denotes the “size” of the spheres. When the distance between two spheres is less than \(\sigma \), the energy is infinite; otherwise it is 0. This is the simplest way to enforce excluded volume among spherical particles in an MC simulation. We will see an implementation of a hard-sphere liquid in the next section.
The most celebrated pair potential is the Lennard-Jones potential (Fig. 5): \begin{equation} \label {eq:lj} U_{LJ}\left (r\right ) = 4\epsilon \left [\left (\frac {\sigma }{r}\right )^{12}-\left (\frac {\sigma }{r}\right )^{6}\right ] \end{equation} \(\epsilon \) is the unit of energy, and is the well-depth of the pair potential (that is, it is the minimum value of the potential). \(\sigma \) is the unit of length, and is the separation distance when the potential is identically zero. Unlike the hard-sphere potential (Eq. 68), the Lennard-Jones potential is “smooth” (it has a continuous first derivative). This smoothness makes it useful in Molecular Dynamics simulations, as we will see later in the course.
We will encounter more complex potentials than Lennard-Jones, but it serves as a useful tool for introductory molecular simulation.
Reduced Units. Because the Lennard-Jones potential is so prevalent in molecular simulations, it is essential that we understand the unit system most often chosen for simulations using this potential. For computational simplicity, energy in a Lennard-Jones system is measured in units of \(\epsilon \) and length in \(\sigma \). This means that everywhere in the code you would expect to see \(\epsilon \) or \(\sigma \), you find a 1 (often implied).
Now, to compute the total potential \(U\) for a system of particles, the simplest algorithm is to loop over all unique pairs of particles. Here is a simple pair search C function (the so-called \(N^2\) algorithm because its complexity scales like \(N^2\)) to compute the total energy of a system of Lennard-Jones particles, in reduced units:
double total_e ( double * rx, double * ry, double * rz, int n ) {
int i,j;
double dx, dy, dz, r2, r6, r12;
double e = 0.0;
for (i=0;i<(n-1);i++) {
for (j=i+1;j<n;j++) {
dx = rx[i]-rx[j];
dy = ry[i]-ry[j];
dz = rz[i]-rz[j];
r2 = dx*dx + dy*dy + dz*dz;
r6 = r2*r2*r2;
r12 = r6*r6;
e += 4*(1.0/r12 - 1.0/r6);
}
}
return e;
}
Although it is strictly correct, the \(N^2\) pair search algorithm is inefficient if there is a finite interaction range in the pair potential. Typically in dense liquid simulations, a Lennard-Jones pair potential is truncated at \(r = 2.5\ \sigma \). There are some potentials that are cut off at even shorter distances. The point is that, when the maximum interaction distance is finite, each particle has a finite maximum number of interaction partners. (This assumes number density is bounded, which is a reasonable assumption.) More advanced techniques (which we discuss later) can be invoked to make the pair search much more efficient in this case. The two most common are (1) the Verlet list, and (2) the link-cell list. For now, and to keep the presentation simple, we will stick to the inefficient, brute-force \(N^2\) algorithm.
4.3.3 Trial Moves
Particle Displacement. The most common trial move in continuous-space MC is a particle displacement. First, a small number \(\Delta R\), representing a maximum displacement, is set. A trial move consists of
- 1.
- Randomly select a particle, \(i\).
- 2.
- Displace x-position coordinate of particle \(i\) by a random amount, \(\delta x\), which is given by \begin{equation} \delta x = \Delta R \xi _x \end{equation} where \(\xi _x\) is a uniform random variate on the interval [-0.5:0.5].
- 3.
- Repeat for the \(y\) and \(z\) coordinates, if applicable.
This move guarantees detailed balance, provided that the random particle selection is uniform; for any given move, selection of all possible particles is equally likely. This means that probability of suggesting a move that displaces a particle, going from a state \(n\) to a new state \(m\), has the same probability of selecting the same particle while in state \(m\) and giving it a displacement that will return the configuration to state \(n\). (Do you think such sequential moves ever actually happen?)
For a system of simple particles, random displacements are the only necessary trial moves; thus, \(\alpha _{nm}\) is always unity. For more complicated systems, there are zoos of trial moves all over the literature. We will consider some more complicated systems and trial moves later in the course; one that we consider next is rigid rotation.
The question at this point is, how does one choose an appropriate value for \(\Delta R\)? If \(\Delta R\) is too small, the system will not explore phase space given a reasonable amount of computational effort. If it is too large, displacements will rarely result in new configurations which will be accepted in a Metropolis MC scheme. So it takes a bit of trial and error to find a good value for \(\Delta R\), and the rule of thumb is to set \(\Delta R\) such that the average probability of accepting a new configuration during a run is about 30%. This is termed “tuning \(\Delta R\) to achieve a 30% acceptance ratio.” We will go through the exercise of determining such an appropriate value for \(\Delta R\) for a simple continuous-space system; namely, 2D hard disks confined to a circle.
Rigid rotation. A second common type of trial move is used in systems of more structured molecules than just simple, single-center spheres. Consider a diatomic with a rigid bond length \(r_0\). Clearly, attempting to move one of the two members of the diatomic by a random displacement is likely to result in a new bond length with may be significantly different from \(r_0\). So, for a system of diatomics, a reasonable set of trial moves would include
- 1.
- Small displacement of molecule center of mass; and
- 2.
- Small rotation around molecule center of mass.
With more than one kind of move, an attempt to generate a new state must be preceded by a random selection of the trial move. We can weight each kind of move and then use a random number to decide which move to attempt. For example, let’s say that we choose that 80% of all trial moves be displacements, and the balance rotations (we will see later whether or not this is a good choice). Prior to an attempted move, we select a uniform random variate, \(\xi \), on the interval \([0,1]\). If \(\xi < 0.8\), which it will be 80% of the time, we execute a displacement of a randomly chosen molecule; otherwise, we execute a rotation of a randomly chosen molecule.
4.4 Case Study 2: MC of Hard Disks
Change directory into your instructional-codes repository and issue a pull if you don’t already see hdisk.c
there. This code simulates disks confined to a circle. The Hamiltonian for this system may be expressed as \begin{equation} \mathscr {H} = \sum _{i=1}^N H_1({\bf r}_i) + \sum _{i=1}^{N}\sum _{j=i+1}^{N} H_2({\bf r}_i,{\bf r}_j) \end{equation}
where \begin{equation} H_1({\bf r}_i) = \left \{\begin {array}{ll} 0 & r_i < R\\ \infty & r_i > R \end {array}\right . \end{equation} and \begin{equation} H_2({\bf r}_i,{\bf r}_j) = \left \{\begin {array}{ll} 0 & \sqrt {\left ({\bf r}_i-{\bf r}_j\right )^2} > \sigma \\ \infty & \sqrt {\left ({\bf r}_i-{\bf r}_j\right )^2} < \sigma \end {array}\right . \end{equation}
\(H_1\) acts to keep the particles confined, and \(H_2\) prevents them from overlapping. One nice thing about using hard-disk Hamiltonians is that there is never a reason to evaluate a Boltzmann factor. Any trial move that results in an overlap or a particle crossing the boundary gives an “infinite” \(\Delta {U}\), so \(e^{-\beta \Delta {U}}\) is identically 0 and the trial is unconditionally rejected.
hdisk.c requires as user input any two of the following three parameters: \(R\), the radius of the circle in \(\sigma \); \(\rho \), the
areal number density of particles (# per square \(\sigma \)); and \(N\), the number of particles. The user may also specify \(\delta r\), the
scalar displacement, and nc, the number of MC cycles, where one cycle is \(N\) attempted particle displacements.
Optionally, the code can generate configurational samples as simple text files or as a single XYZ-formatted
trajectory (which looks like a concatenation of XYZ files), and in order to generate these samples,
traj_samp must be set greater than zero. The code reports the acceptance ratio, among other things:
\[ \left [\begin {array}{c} \mbox {acceptance}\\ \mbox {ratio} \end {array}\right ] = \left [\frac {\mbox {number of successful trials}}{\mbox {number of trials}}\right ] \]
One important aspect of any MC simulation code is how the particle positions are initialized. Here, it is best
to assign initial positions to the particles such that the initial energy is 0 (i.e., there are no overlaps nor particles
out of bounds.) Try to figure out how the function init() in the program hdisk.c accomplishes
this.
As a suggested further exercise, use hdisk.c to determine a reasonable displacement to achieve a 30%
acceptance ratio at a density of 0.5. Compare your results across differently sized systems and runs with
different numbers of cycles. For fewer than 10\(^6\) cycles, you will have large acceptance ratios because the initial
condition is not yet fully destroyed.
Below is a plot of acceptance ratio vs. \(\Delta R\) for densities \(\rho \) of 0.2, 0.4, 0.6, from a simulations of 200 particles. 2,000 cycles were performed for each run, and each had a unique seed. Are your results consistent with this data?
One advantage of the XYZ format is that we can use VMD to visualize our configurations, and even to make animations of our simulations. For example, suppose we generate a short 1,000-cycle MC trajectory of the hard-disk system at \(\rho \) = 0.7 and \(N\) = 100:
cd cd cheT580-202035/instructional-codes/my_work mkdir hdisk_run cd hdisk_run gcc -O3 -o hdisk ../../originals/hdisk.c -lm -lgsl ./hdisk -xyz traj.xyz -traj_samp 1 -nc 1000 -rho 0.7 -N 100 -dr 0.5 # R = 6.74336; rho = 0.70000; N = 100; seed = 23410981 Results: Number of Trial Moves: 100000 Maximum Displacement Length: 0.50000 Acceptance Ratio: 0.48499 Reject Fraction Out-of-bounds: 0.09707 Reject Fraction Overlap: 0.90293 ls hdisk traj.xyz
Notice I provide the name of the trajectory file, indicating I want one sample per cycle. I also set the magnitude of the maximum displacement at 0.5 \(\sigma \).
After the run finishes, the file traj.xyz appears. It looks like this:
head -10 traj.xyz 100 Generated by hdisk.c; all z-components are zero; all elements are H H 2.48776 1.78242 0.00000 H 5.80283 0.16579 0.00000 H 2.57789 -4.28715 0.00000 H -4.25852 0.68353 0.00000 H -5.95301 2.69712 0.00000 H 0.28903 -1.04400 0.00000 H 1.86426 -0.93249 0.00000 H -4.16984 1.87756 0.00000
We can use VMD to visualize this trajectory. If you are on a Windows 10 machine and you installed the Windows version of VMD, you can simply launch it from the start menu, and then navigate to
\\wsl\$\<distro>\home\<username>\cheT580-202035\instructional-codes\my_work\hdisk_run\
and then just click on traj.xyz to read it in. If you are on a Mac, you should be able to navigate
straight to the file. In Fig. 7, I show three snapshots from this simulation, at cycle 0, 500, and
1000.
As another suggested exercise, consider generating a trial move in the following way:
- 1.
- Randomly select a particle \(i\).
- 2.
- Randomly choose a direction. In 2D, this is an angle \(\phi \) chosen uniformly from the interval \([0,2\pi ]\). In 3D, this is two angles \(\theta \) and \(\phi \), where \(\cos \theta \) is chosen uniformly from the interval \([-1,1]\), and \(\phi \) is chosen uniformly from the interval \([0,2\pi ]\).
- 3.
- Compute the displacement vector components. In 2D: \(dx = \Delta R\cos \theta \) and \(dy = \Delta R\sin \theta \). In 3D: \(dx = \Delta R\sin \theta \cos \phi \), \(dy = \Delta R\sin \theta \sin \phi \), \(dz = \Delta R\cos \theta \).
- 4.
- Execute the move by adding the displacement vector components to the position of particle \(i\).
The main difference of this scheme with the first one is that, in 2D, only one random number must be chosen per displacement attempt. Explore the acceptance ratio vs. \(\Delta R\) relationship with this new trial move. Does this scheme generate a 30% acceptance ratio for a higher or lower maximum displacement that does the original scheme?
4.5 Case Study 3: Hard-Disk Dumbbells in 2D
In this case study, we consider a slightly different system than the hard-disk system in the previous case study. Here, we imagine that pairs of disks are tethered together to form dumbbells. The “bond-length” of a dumbbell is a constant parameter, \(r_0\). We will again confine the dumbbells to a circle.
What is new is how we have to consider the trial moves for this system. We cannot simply select a random particle and try to displace it, because this is likely to violate the constant bond-length of the dumbbell that particle belongs to. How then do we generate new configurations? A simple idea is to use two kinds of trial moves, translation of entire dumbbells and rotation of dumbbells around their centers of mass. This was originally presented in Sec. 4.3.3. In order to implement an MC code with more than one trial move, we must include a “trial move selection rule” which randomly selects a trial move based on their user-defined “weights”.
The code hdisk-dumbbells.c implements a MC simulation of hard-disk dumbbells. One specifies a
number of particles using -N # at the command line, and the number of dimers is assumed to be
\(N\)/2. One can specify two of \(N\), \(\rho \) (areal particle density), or \(R\) (confining domain diameter). One can
also specify \(\delta r\) (maximum dimer displacement distance) and \(\delta \phi \) (maximum dimer rotation angle), as
well as the dimer bond length \(r_0\). An XYZ-format trajectory can be saved every -fs cycles using
-traj my_traj.xyz.
Let’s run this program for \(N\) = 200, \(\rho \) = 0.6, \(\delta r\) = 1, \(\delta \phi = \pi \), for 10,000 cycles, saving 1000 snapshots in
traj.xyz:
cd cd cheT580-202035/instructional-codes/my_work mkdir hddb_run cd hddb_run gcc -O3 -o hddb../../originals/hdisk-dumbbells.c -lm -lgsl ./hddb -rho 0.6 -dr 1 -da 3.1 -dw 0.5 -N 200 -traj traj.xyz -fs 10 -nc 10000 # R = 10.30; rho = 0.60; N = 200; r_0 = 1.00; s = 1.00; disp_wt = 0.50, seed = 23410981 Results: Number of trial moves: 2000000 Maximum displacement length: 1.000 Number of displacement attempts: 999499 Maximum rotation angle (radians): 3.100 Number of rotation attempts: 1000501 Displacement acceptance ratio: 0.344 Rotation acceptance ratio: 0.405 Reject Fraction Out-of-bounds: 0.10037 Reject Fraction Overlap: 0.89963 Trajectory saved to: traj.xyz ls hddb traj.xyz
As with the hard disks, we can use VMD to visualize frames in this trajectory. Fig. 8 shows three representative snapshots from this simulation, along with a special view showing the histories of four particular dimers. This last rendering serves to indicate that dimers have mostly explored the entire domain for this number of cycles (is this the case for 1,000 cycles?).
Consider the following question: Does the acceptance ratio of rotational moves depend upon the weight given to displacement moves? Why or why not? Below is a plot of the acceptance ratio vs. the maximum displacement for a system of 100 dumbbells at a density of 0.5, for various displacement move weights between 0.1 and 0.9. As you can see, there appears to be no effect on the acceptance of trial displacements if we change how frequently we perform them relative to trial rotations.
Let’s consider computing an observable function that describes how the molecules order themselves in the circular domain. One such meaningful function we’ll call \(\langle P[\theta (r)]\rangle \), where \begin{equation} P(\theta ) = \cos ^2\theta -\frac {1}{2} \end{equation} and \(\theta \) is defined as the angle between a dumbbell’s bond and a vector joining its midpoint to the domain origin. \(\langle P[\theta (r)]\rangle \) is therefore the expectation of \(P\) averaged over all dumbbells located between \(r\) and \(r+\delta r\) from the origin. Why is this a meaningful measure of order for this system? Consider that dumbbells near the periphery like to align so their bonds are tangent to the periphery itself and therefore perpendicular to the radius, and the further away from the periphery you look, the less likely this order is to appear. When two vectors are perpendicular, their angle betwen them is \(\pi \) and since \(\cos \pi =0\), we expect \(\langle P\rangle \) to be -0.5. Is this what we observe for our system?
To explore this function, the code hdisk_dumbbells_order.c was copied from hdisk_dumbbells.c. In
hdisk_dumbbells_order.c, I introduce the arrays theta[] to hold a tally of \(\cos ^2\theta \) values and Rcount[] to hold a
count of hits in each radial bin, so that \(\langle P(\theta (r)) \rangle \) = theta[i]/Rcount[i] - 0.5, where i is the radial bin index set by \(R\)
and a specified number of bins. Fig. 10 shows \(P\) vs \(r\) for a hard-disk dumbbell system at density 0.66 with 400
particles (200 dumbbells), run for 500,000 cycles. We can see some interesting structure here,
notably that it looks like, on average, there is very little ordering at the periphery and it becomes
substatially higher as we get closer to the origin, though still not very strong. Notably, at about 1 \(\sigma \)
from the periphery, we see a dip in \(P\) that might indicate an enrichment in tangentially-oriented
dumbbells.
A final note: Since this is a 2D system, we don’t use the second Legendre polynomial of \(\cos \theta \): \begin{equation}\label {eq:slp} P_2(x)=\frac 12(3x^2-1) \end{equation} The reason for this is clear. Suppose \(f(\theta )\) is the probability distribution of the angle \(\theta \). Now, since the molecules are not polar (they have no head or tail), \(\theta \) is meaningful only on the domain \([0,\frac {\pi }{2}]\). Suppose further that there is no preferred angle; i.e., \(f\) is a constant. In 2D, normalization requires \begin{equation} 1=f\int _0^{\pi /2}d\theta \ \rightarrow \ f=\frac {2}{\pi }, \end{equation} and in this case, with no preferred orientation, the average of \(\cos ^2\theta \) is \begin{equation} \left <\cos ^2\theta \right > = \frac {2}{\pi }\int _0^{\pi /2} \cos ^2\theta d\theta = \frac {2}{\pi }\frac {\pi }{4} = \frac 12. \end{equation} And Eq. 75 will have the described behavior. The second Legendre polynomial evaluates to zero when \(\cos ^2\theta \) is \(\frac 13\), which is indeed what \(\left <\cos ^2\theta \right >\) evaluates to when \(f\) is a constant in three dimensions.
4.6 Case Study 4: Equation of State of the Lennard-Jones Fluid
The final case study we will consider in this unit on Monte Carlo simulation is the prototypical system for continuous-space, 3D liquids: The Lennard-Jones fluid. (This is detailed in Sec. 3.4, “Case Study 1” in Frenkel & Smit [1].) The primary objective of the MC code is to predict the pressure of a sample of Lennard-Jonesium at a given density and temperature; that is, we can use MC to map out (in principle) the phase diagram of a material. We will use this case study to introduce and discuss another important element of a large number of molecular simulations: periodic boundary conditions.
We would like to simulate bulk fluid. The apparently simplest way to approximate bulk behavior in a finite number of particles is to employ periodic boundaries. That is, we imagine the box of length \(L\) is embedded in an infinite space tiled with replicas of the central box. If we focus on the central box, and watch as one particle is displaced “out” of the box, it will reappear in the box at the opposite face. Moreover, particles interact with “images” of other particles in all replica boxes. Periodic boundaries thus allow us to mimic the infinite extent of bulk fluid.
Periodic boundaries require the use of the minimum image convention (MIC) when computing
inter-particle contributions to the total energy. Below is a modified \(N^2\) loop for a 3-D system of point
particles obeying the Lennard-Jones pair potential with periodic boundaries. The function e_i()
computes the sum of all pair interactions between a stipulated particle i and all particles from i0
to N-1. It is called from total_e() such that a sum of pairwise energies for all unique pairs is
computed.
double e_i ( int i, double * rx, double * ry, double * rz, int N,
double L, double rc2, int tailcorr, double ecor,
int shift, double ecut, double * vir, int i0 ) {
int j;
double dx, dy, dz, r2, r6i;
double e = 0.0, hL=L/2.0;
*vir=0.0;
for (j=i0;j<N;j++) {
if (i!=j) {
dx = rx[i]-rx[j];
dy = ry[i]-ry[j];
dz = rz[i]-rz[j];
if (dx>hL) dx-=L;
else if (dx<-hL) dx+=L;
if (dy>hL) dy-=L;
else if (dy<-hL) dy+=L;
if (dz>hL) dz-=L;
else if (dz<-hL) dz+=L;
r2 = dx*dx + dy*dy + dz*dz;
if (r2<rc2) {
r6i = 1.0/(r2*r2*r2);
e += 4*(r6i*r6i - r6i) - (shift?ecut:0.0);
*vir += 48*(r6i*r6i-0.5*r6i);
}
}
}
return e+(tailcorr?ecor:0.0);
}
double total_e ( double * rx, double * ry, double * rz,
int N, double L,
double rc2, int tailcorr, double ecor,
int shift, double ecut, double * vir ) {
int i;
double tvir;
double e = 0.0;
*vir=0.0;
for (i=0;i<N-1;i++) {
e += e_i(i,rx,ry,rz,N,L,rc2,
tailcorr,ecor,shift,ecut,&tvir,i+1);
*vir += tvir;
}
return e;
}
The function e_i() is also useful when determining whether or not to accept a trial displacement move,
since the change in energy \(\Delta \mathscr {U}\) associated with moving particle \(i\) can be found by calling e_i() for particle \(i\) before
and after the move; the latter energy minus the former is \(\Delta \mathscr {U}\), since no other particles are displaced. This is much
faster than just doing a full-blown \(N^2\) calculation to evaluate each trial move, but it forces you to do careful
accounting to keep the total energy up-to-date. If the move is accepted the total energy must be incremented
by \(\Delta \mathscr {U}\); if the virial is also being tallied, the virial must also be similarly incremented by the change in the virial upon
particle displacement.
Note that each \(i\)-\(j\) displacement component is subject to the MIC; if \(\Delta x\) is greater than \(L/2\), we subtract \(L\) from it; if it is less than \(-L/2\), we add \(L\) to it; same for \(\Delta y\) and \(\Delta z\). Many caveats come with using periodic boundaries. (A thorough discussion appears in Sec. 3.2.2 of F&S.) The first thing to realize is that the total potential per particle (as a sum of pair potentials) in principle diverges in an infinite periodic system. This can be circumvented by introducing a finite interaction range to the pair potential. We usually work with systems large enough such that the cutoff of the pair potential, \(r_c\), is less than one-half the box-length, \(L\), in a cubic box. This means that the “image” interactions involve only immediately neighboring replicas.
Truncation of a pair potential is an important idea to understand. The major point is that the cutoff must be spherically symmetric; that is, we can’t simply cut off interactions beyond a box length in each direction, because this results in a directional bias in the interaction range of the potential. So, a hard cutoff radius \(r_c\) is required, and it should be less than half the box length. The secondary point is that, once \(r_c\) is chosen, if you wish to mimic a potential with infinite range, you must use the correction terms for energy and pressure described below.
The system we consider is made of \(N\) particles which interact via the Lennard-Jones pair potential
(Eq. 69). The particles are confined in a cubic box with side-length \(L\). Length is measured in units of \(\sigma \)
and energy in \(\epsilon \), and we consider particles with 1 \(\sigma \) diameters. A code is provided for simulating
this system using Metropolis Monte Carlo: mclj.c). mjlc.c will compute the pressure given a
temperature and density in the manner discussed in the text. If a cutoff radius is chosen by the user,
then a truncated and shifted pair potential is used, and the following tail corrections are applied:
The pressure is computed from \begin{equation} \label {eq:mc:P} P = \rho T + {\rm vir}/V \end{equation} where \(\rm vir\) is the virial: \begin{equation} {\rm vir} = \frac {1}{3}\sum _{i>j} {\bf f}\left (r_{ij}\right )\cdot {\bf r}_{ij} \end{equation} and \(V\) is the system volume. \({\bf f}\left (r_{ij}\right )\) is the force exerted on particle \(i\) by particle \(j\), defined as the negative gradient of the pair potential \(u\left (r_{ij}\right )\) with respect to the position of particle \(i\):
Here we have made use of the fact that \begin{equation} \frac {\partial X}{\partial {\bf r}} = \frac {\bf r}{r}\frac {\partial X}{\partial r} \end{equation} when operating on a function \(X\) which depends on relative particle separations. For the Lennard-Jones pair potential (Eq. 69), we see that \begin{equation} \frac {\partial u\left (r_{ij}\right )}{\partial r_{ij}} = 4\epsilon \left [-12\frac {\sigma ^{12}}{r_{ij}^{13}} + 6\frac {\sigma ^{6}}{r_{ij}^7}\right ] \end{equation} So, \begin{equation} {\bf f}_{ij}\left (r_{ij}\right ) = \frac {{\bf r}_{ij}}{r_{ij}^2}\left \{48\epsilon \left [\left (\frac {\sigma }{r_{ij}}\right )^{12} - \frac {1}{2}\left (\frac {\sigma }{r_{ij}}\right )^6\right ]\right \}, \end{equation} and therefore, \begin{equation} {\rm vir} = \frac {1}{3}\sum _{i>j} \left \{48\epsilon \left [\left (\frac {\sigma }{r_{ij}}\right )^{12} - \frac {1}{2}\left (\frac {\sigma }{r_{ij}}\right )^6\right ]\right \} \end{equation} Notice that any particular pair’s contribution to the total virial is positive if the members of the pair are repelling one another (positive \(f\) along \({\bf r}_{ij})\)), and negative if the particles are attracting one another.
If you read the code mclj.c, you should see that the initialization of positions is accomplished by putting
the particles on cubic lattice sites such that an overall density is achieved. It is therefore convenient to run
simulations with numbers of particles that are perfect cubes, such as 128, 216, 512, etc, so that the initial state
uniformly fills the box.
Another consideration is that a certain number of cycles should be “burned” prior to gathering statistics so
this initial state is fully erased. The flag -ne allows the user to specify how many equilibration cycles are to be
performed before switching to “production” mode.
Fig. 12 shows two snapshots made with VMD of the Lennard-Jones systems with 512 particles.
As a suggested exercise, you can use mclj.c to try to reproduce Figure 3.5 in F&S, which shows \(P\) vs. \(\rho \) at
both \(T\) = 2.0 and \(T\) = 0.9. How many cycles do you need? How many equilibration cycles? What maximum
displacement did you choose?
Below are some of my preliminary results using the code mclj.c. I used only 5,000 cycles for 512 particles
for each point, and each point is the result of a single run. These numbers appear to compare well with
those in Figure 3.5 in F&S, for which we have no idea how many cycles or independent runs were
performed.
