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. 83). 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:

$$u^{\rm tail} = \frac{8}{3}\pi\rho\epsilon\sigma^3\left[\frac{1}{3}\left(\frac{\sigma}{r_c}\right)^{9}-\left(\frac{\sigma}{r_c}\right)^{3}\right]$$ (92)

$$\Delta P^{\rm tail} = \frac{16}{3}\pi\rho^2\epsilon\sigma^3\left[\frac{2}{3}\left(\frac{\sigma}{r_c}\right)^{9}-\left(\frac{\sigma}{r_c}\right)^{3}\right]$$ (93)

The pressure is computed from

$$P = \rho T + {\rm vir}/V$$ (94)

where ${\rm vir}$ is the virial:

$${\rm vir} = \frac{1}{3}\sum_{i>j} {\bf f}\left(r_{ij}\right)\cdot{\bf r}_{ij}$$ (95)

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$:

$${\bf f}\left(r_{ij}\right) = -\frac{\partial u\left(r_{ij}\right)}{\partial {\bf r}_i}$$ (96)

$$= -\frac{{\bf r}_{ij}}{r_{ij}}\frac{\partial u\left(r_{ij}\right)}{\partial r_{ij}}$$ (97)

Here we have made use of the fact that

$$\frac{\partial X}{\partial {\bf r}} = \frac{\bf r}{r}\frac{\partial X}{\partial r}$$ (98)

when operating on a function $X$ which depends on relative particle separations. For the Lennard-Jones pair potential (Eq. 83), we see that

$$\frac{\partial u\left(r_{ij}\right)}{\partial r_{ij}} = 4\ep... ...sigma^{12}}{r_{ij}^{13}} + 6\frac{\sigma^{6}}{r_{ij}^7}\right]$$ (99)

So,

$${\bf f}_{ij}\left(r_{ij}\right) = \frac{{\bf r}_{ij}}{r_{ij}... ...frac{1}{2}\left(\frac{\sigma}{r_{ij}}\right)^6\right]\right\},$$ (100)

and therefore,

$${\rm vir} = \frac{1}{3}\sum_{i>j} \left\{48\epsilon\left[\le... ...\frac{1}{2}\left(\frac{\sigma}{r_{ij}}\right)^6\right]\right\}$$ (101)

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.