3 Monte Carlo Simulation

The first simulation technique we will study is the Monte Carlo method, covered in Ch. 3 of Frenkel & Smit [1]. In these pages, I provide some elaboration of the content in the book, but it should be emphasized that a complete presentation of MC involves most of what is presented in Ch. 3. The primary “flavor” of MC most appropriate for this introductory level course is the original Metropolis method, explained in Sec. 3.1.

We will also use MC to explain basic technical aspects of molecular simulation code in Sec. 3.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 system (Sec. 3.2) which does not appear in Frankel & Smit. The second considers hard-disks confined in a circle (Sec. 3.4), and the third is hard-disk-dumbbells (Sec. 3.5). The fourth is “Case Study 1” from the Frankel & Smit, namely using MC to explore the equation of state of a model liquid known as the Lennard-Jones fluid (Sec. 3.6).

3.1 The Metropolis Monte Carlo Method

Monte Carlo is in general a technique for performing 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 average of the quantity \(f/\rho \) on \(\rho \) in the interval \([a,b]\): \begin{equation} I = \left \langle \frac {f\left (x\right )}{\rho \left (x\right )}\right \rangle _\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 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 more traditional 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} \frac {\pi }{4} = \int _0^1 \left (1-x^2\right )^{1/2}dx \end{equation} Allen and Tildesley [2] mention that, in order use Eq. 43 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 required three orders of magnitude fewer points to discretize the interval to obtain an accuracy of one part in 10\(^7\). So the answer is, integral estimation using 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 100 particles having 300 coordinates, the configurational average \(\left \langle G\right \rangle \) (Eq. 35) could be discretized using Simpson’s rule. If we did that, requesting only a modest 10 points per axis in configurational space, we would need to evaluate the integrand \(Ge^{-\beta \mathscr {U}}\)

10\(^{300}\)

times. This is an almost unimaginably large number (a googul cubed). 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 \langle G\right \rangle \) 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 \langle G\rho _{NVT}/\rho _{NVT}\right \rangle \approx \left \langle G\right \rangle , \end{equation} which should give an excellent approximation for \(\left \langle G\right \rangle \). The idea of using the \(\rho _{NVT}\) as the sampling distribution is due to Metropolis et al. [6]. This makes the real work in computing \(\left \langle G\right \rangle \) generating states that randomly sample \(\rho _{NVT}\).

Metropolis et al. [6] 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 called “state space,” 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. There is a large variety of such “trial moves” for any particular system; we will only deal with a few simple ones in this course.

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. 48 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. 49 constrains our choice of \(\bf \pi \). This means there is still more than one way to specify \(\bf \rho \). Metropolis et al. [6] 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; this fact is of little importance in this course.

Metropolis et al. [6] 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 \mathscr {U}\left (\Gamma _m\right )}}{e^{-\beta \mathscr {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 [\mathscr {U}\left (\Gamma _m\right )-\mathscr {U}\left (\Gamma _n\right )\right ]\right \} \equiv \exp \left (-\beta \Delta \mathscr {U}_{nm}\right ) \end{equation} where we have defined the change in potential energy as \begin{equation} \Delta \mathscr {U}_{nm} = \mathscr {U}\left (\Gamma _m\right )-\mathscr {U}\left (\Gamma _n\right ) \end{equation}

There are many choices for \({\rm acc}\left (n \rightarrow m\right )\) that satisfy Eq. 55. 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 \mathscr {U}_{nm}\right ) & \Delta \mathscr {U}_{nm} > 0\\ 1 & \Delta \mathscr {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.

3.2 Case Study 1: The 2D Ising Magnet

3.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 a Ising system of noninteracting spins, because each spin made its own independent contribution to the Hamiltonian.

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:

where the notation \(\left \langle ij\right \rangle \) denotes that we are summing over unique pairs of nearest neighbors, and, as before a spin \(s_i\) is either +1 or -1. 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 neigbhors 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 \langle M\right \rangle = \frac {\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 ]}{\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 intution tells us that as \(T\rightarrow \) 0, \(\left \langle M\right \rangle \rightarrow \) 1, and as \(T\rightarrow \infty \), \(\left \langle M\right \rangle \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 positive 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 lattice 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 a visit to the website of Peter Young at UC Santa Cruz. You will see a Java implentation 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 system. 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.

3.2.2 A C Code for the 2D Ising Magnet

In this section, we will dissect piece-by-piece a small program (written in C) which implements an NVT Metropolis Monte Carlo simulation of a 2D Ising lattice. Click here to download the code. You can compile the code using the command

cfa@abrams01:/home/cfa> gcc -O3 -o ising ising.c -lm -lgsl

and you can then run it at the command line as ./ising. The code will conduct a canoncial Metropolis Monte Carlo simulation of an Ising lattice 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 \langle \epsilon \right \rangle \) and the average spin value, \(\left \langle s_1\right \rangle \). Periodic boundaries are employed in calculating the nearest-neighbor interactions.

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.

_________________________________________________________________________________________

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_rng.h>

int E ( int ** F, int L, int i, int j ) {
  return -2*(F[i][j])*
     (F[i?(i-1):(L-1)][j]+F[(i+1)%L][j]+
             F[i][j?(j-1):(L-1)]+F[i][(j+1)%L]);
}

void samp ( int ** F, int L, double * s, double * e ) {
  int i,j;

  *s=0.0;
  *e=0.0;
  for (i=0;i<L;i++) {
    for (j=0;j<L;j++) {
      *s+=(double)F[i][j];
      *e-=(double)(F[i][j])*
   (F[i][(j+1)%L]+F[(i+1)%L][j]);
    }
  }
  *s/=(L*L);
  *e/=(L*L);
}

void init ( int ** F, int L, gsl_rng * r ) {
  int i,j;
  for (i=0;i<L;i++) {
    for (j=0;j<L;j++) {
      F[i][j]=2*(int)gsl_rng_uniform_int(r,2)-1;
    }
  }
}

int main (int argc, char * argv[]) {

  int ** F;
  int L = 20;
  int N;
  double T = 1.0;

  int nCycles = 1000000;
  int fSamp = 1000;

  int nSamp;
  int de;
  double b;
  double x;
  int i,j,a,c;

  double s=0.0, ssum=0.0;
  double e=0.0, esum=0.0;

  gsl_rng * r = gsl_rng_alloc(gsl_rng_mt19937);
  unsigned long int Seed = 23410981;

  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]);
  }

  gsl_rng_set(r,Seed);

  N=L*L;
  F=(int**)malloc(L*sizeof(int*));
  for (i=0;i<L;i++) F[i]=(int*)malloc(L*sizeof(int));

  init(F,L,r);
  T=1.0/T;

  for (c=0;c<nCycles;c++) {

    for (a=0;a<N;a++) {

      i=(int)gsl_rng_uniform_int(r,L);
      j=(int)gsl_rng_uniform_int(r,L);

      de = E(F,L,i,j);
      b = exp(de*T);
      x = gsl_rng_uniform(r);

      if (x<b) {
 F[i][j]*=-1;
      }

      if (gsl_rng_uniform(r)<exp(E(F,L,i,j)*T)) F[i][j]*=-1;

    }

    if (!(c%fSamp)) {
      samp(F,L,&s,&e);
      ssum+=s;
      esum+=e;
      nSamp++;
    }

  }

  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");
}

3.2.3 Suggested Exercises

1. Run the code for the following values of temperature: 5.0, 4.0, 3.0, 2.0, 1.0. Run the code several times at each \(T\) with a different value for the random number generator seed. Report the average spin and average energy per spin. What is happening near \(T = 2.0\)?
2. Modify the code so that when samples are taken in accumulating statistics for \(\left \langle s_1\right \rangle \) and \(\left \langle \epsilon \right \rangle \), 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);
3. 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?
4. (Advanced) Modify the code ising.c to compute the quantity \(\left \langle s_is_j\right \rangle - \left \langle s_i\right \rangle \left \langle s_j\right \rangle \) as a function of various distances between spins \(i\) and \(j\).

3.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 \(\mathscr {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

3.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 FORTRAN and C, we might declare these arrays for 1,000 particles as

      parameter (N = 1000)
      real*8 rx(N), ry(N), rz(N)

  const 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 code fragment to write a configuration of \(N\) atoms appears below:

 do 100 i=1,N
    write(9,101) i,rx(i),ry(i),rz(i)
100     continue
101 format(i,e13.8,e13.8,e13.8)
! writes to file fort.9

  fp = fopen("file1.dat","w");
  for (i=0;i<N;i++)
    fprintf(fp,"%i %.8le %.8le %.8le\n",
     i,rx[i],ry[i],rz[i]);
/* fp is declared as a pointer of type FILE */

It’s quite easy to read ASCII data in:

 do 100 i=1,N
    read(9,101) i,rx(i),ry(i),rz(i)
100     continue
101 format(i,e13.8,e13.8,e13.8)
! read from file fort.9

  fp = fopen("file1.dat","r");
  for (i=0;i<N;i++)
    fscanf(fp,"%i %.8le %.8le %.8le\n",
     &i,&rx[i],&ry[i],&rz[i]);
/* fp is declared as a pointer of type FILE */

3.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} \mathscr {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 \mathscr {U}}\) correlates every position with every other position, and the configurational integral (Eq. 35) 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} \mathscr {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: \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. 63), 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.

The Lennard-Jones pair potential.

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 simulation, 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.

Now, to compute the total potential \(\mathscr {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.

3.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.

Details of the rigid body rotation implementation are presented in the text, in Sec. 3.3.2. The two techniques described are vector addition for linear molecules, and quaternions for nonlinear molecules. For the purposes of one of our simple 2D case-studies, we will consider rotation by a small angle, which is effectively equivalent to the vector addition method described by Frenkel and Smit.

3.4 Case Study 2: MC of Hard Disks

Download the code hdisk.c by clicking here. 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 \mathscr {U}\), so \(e^{-\beta \Delta \mathscr {U}}\) is identically 0 and the trial is unconditionally rejected.

hdisk.c accepts as user input any two of the following three parameters: \(R\), the radius of the circle; \(\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 nCycle, the number of MC cycles, where one cycle is \(N\) attempted particle displacements. The only output of the code currently is the acceptance ratio: \[ \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. 10,000 cycles were performed for each run. Are your results consistent with this data?

Acceptance ratio vs. \(\Delta R\) for various densities in a 200-particle hard-disk system from 100,000-cycle MC simulations.

As a third 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 \(\theta \) 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?

3.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.

Configuration snapshot of a system of hard-disk-dumbbells. \(N\) = 200 particles, \(rho\) = 0.5, and \(R\) = 11.3.

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. 3.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 hdb.c implements a MC simulation of hard-disk dumbbells. 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.

Acceptance ratio vs. maximum dumbbell displacement for various displacement trial move weights between 0.1 and 0.9. \(N\) = 200 (100 dumbbells), \(\rho \) = 0.5, and \(R\) = 11.3.

As a second suggested exercise (one of an advanced nature, suitable for a course project), we can use this code to explore the liquid crystalline nature of the dumbbell fluid. Because the dumbbells are slightly elongated along one direction, they may tend to line up in a dense liquid. For dumbbell \(i\), let the quantity \(\theta _i\) represent the angle made by the interparticle segment and some global coordinate frame axis (say the \(x\) axis). We could use this code to compute the average orientation, \(\left \langle \theta \right \rangle \) as a function of density and temperature. The problem is that if the system is truly a liquid, even if large numbers of dumbbells are lined up at any one time, the system will slowly evolve so that all dumbbell orientations are eventually realized. So \(\left \langle \theta \right \rangle \) is probably not so interesting to calculate. What would be interesting would be to assess the effect of the confining circle on the spatially resolved orientation field, \(\left \langle \theta \right \rangle (x,y)\). You could construct a 2D histogram of orientation and perform MC to populate it. But because the system is nominally cylindrically symmetric, it would suffice to consider a one-dimensional field, \(\left \langle \theta \right \rangle (r)\). Plot \(\left \langle \theta \right \rangle (r)\) vs. \(r\) for various densities and temperatures.

A third suggested (advanced) exercise would be to measure an orientational correlation function. Here, we let the quantity \(\theta _{ij}\) be the relative angle between the segments of dumbbells \(i\) and \(j\). You can use the code to accumulate statistics on \(\theta _{ij}\) as a function of \(r_{ij}\), where \(r_{ij}\) is the center-of-mass-to-center-of-mass distance of the two dumbbells. It is actually better to accumulate statistics of the following polynomial of the cosine of \(\theta _{ij}\): \begin{equation}\label {eq:p2} P\left (\cos \theta _{ij}\right ) = \cos ^2\theta _{ij}-\frac 12. \end{equation} \(\left \langle P\right \rangle \) as a function of \(r\) will reveal the character of the liquid-crystalline-like ordering of the dumbbells; this is the orientational correlation function. When two dumbbells are aligned, \(\cos \theta _{ij}\) is unity, and therefore \(P\) approaches unity. When two dumbbells are perpendicular, \(P\) is -0.5. The average, \(\left \langle P\right \rangle (r)\) at a particular distance \(r\), will range between -0.5 and 1.0, for perfectly perpendicular to perfectly aligned, and 0 implies no preferred orientation. Try to hypothesize how \(\left \langle P(r)\right \rangle \) behaves as one changes temperature and density.

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 \langle \cos ^2\theta \right \rangle = \frac {2}{\pi }\int _0^{\pi /2} \cos ^2\theta d\theta = \frac {2}{\pi }\frac {\pi }{4} = \frac 12. \end{equation} And Eq. 69 will have the described behavior. The second Legendre polynomial evaluates to zero when \(\cos ^2\theta \) is \(\frac 13\), which is indeed what \(\left \langle \cos ^\theta \right \rangle \) evaluates to when \(f\) is a constant in three dimensions.

3.6 Case Study 4 (F&S Case Study 1): 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.

A schematic representation of periodic boundary conditions in two dimensions. The black particle leaves the central box by leaving a through right-hand boundary, and consequently re-enters through the left-hand boundary. The two white particles interact through the boundary.

Many caveats come with using periodic boundaries. A thorough discussion appears in Sec. 3.2.2 of F&S. The first thing to realize, as pointed out by F&S, is that the total potential (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 nieghboring replicas. In the code mclj.c, the user can specify both the box length and the cutoff radius.

Truncation of a pair potential is an important idea to understand, and F&S devote a significant portion of the MC chapter to it. 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 should be less than half a 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. 64). 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. 64), 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. F&S do not discuss this; we will try to determine this by experimentation. The flag -ne allows the user to specify how many equilibration cycles are to be performed before switching into “production” mode.

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 1,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.

Pressure vs. density in a Lennard-Jones fluid, measured by Metropolis MC simulation, at two temperatures, 0.9 and 2.0. Because \(T\) = 0.9 is below the critical temperature, we see negative pressures; as F&S point out, such a small system can sustain negative pressure without phase separating into a liquid/vapor system (there are simply too few particles to nucleate a stable second phase). 1,000 cycles were performed for each run, and the system contained 512 particles.

3.7 Suggest Project: Lennard-Jones Dumbbells

In Sec. 3.5, MC simulation of 2D disk dumbbells was presented to introduce MC simulation using more than one kind of trial move; in that case, displacement of randomly selected molecules and rotation of randomly selected molecules. For this suggested project, we can imagine performing simulations of dumbbells in 3D in a completely analogous manner. Modify the code mclj.c such that every even numbered particle is bonded to its odd-numbered neighbor one up. Each such pair forms a molecule, and trial moves consist of (1) moves of molecule centers by a distance \(\Delta {\bf r}\), and (2) rotations of molecules by random angular displacements \((\Delta \theta ,\Delta \phi )\). (Recall how one must sample polar angles to maintain uniform density of direction in 3D space: one must sample uniformly from \(\cos \theta \)!) An interesting quantity to measure is the orientational correlation function \(P_2\left (\cos \theta _{ij}\right )\) (Eq. 69) as a function of density.