Isothermal-Isobaric

In this section, we consider how to conduct Monte Carlo simulation in ensembles other than the canonical ensemble. In deriving the partition function for the canonical ensemble (Eq. 46), we imagined our sytem of constant $N$, $V$, and $T$ in thermal contact with a large reservoir. This thermal contact allowed the system and reservoir to exchange energy such that the system remained at constant $T$, and what resulted was the Boltzmann factor. In Section 5.4.1, F&S explain the case when we have the reservoir and the system both thermally and mechanically coupled. The mechanical coupling allows the volume of the system to change such that the pressure in the system is the same as the reservoir, which is again considered as an inifinite ideal gas. In addition to thermostatting our system, the reservoir acts as a barostat.

First, for convenience, we express the set of coordinates, ${\bf r}^N$, scaled by the box length, $L$, as ${\bf s}^N$. The partition function in the NPT ensemble is then

$$Q\left(N,P,T\right) = \frac{\beta P}{\Lambda^{3N}N!} \int dV V^N ... ...ght) \int d{\bf s}^N \exp\left[-\beta\mathscr{U}\left({\bf s}^N;L\right)\right]$$ (168)

The free energy associated with this ensemble is the Gibbs free energy:

$$G = -k_BT\ln Q\left(N,P,T\right)$$ (169)

Now, compared to the canonical ensemble, in the NPT ensemble, volume is an additional degree of freedom. We need the probability distribution to include volume:

$$\mathscr{N}\left(V;{\bf s}^N\right) \propto V^N\exp\left(-\beta PV\right)\exp\left[-\beta\mathscr{U}\left({\bf s}^N;L\right)\right]$$ (170)

$$= \exp\left\{-\beta\left[\mathscr{U}\left({\bf s}^N,V\right) + PV - N\beta^{-1}\ln V\right]\right\}$$ (171)

We can use this new Boltzmann factor in an acceptance rule for a Monte Carlo trial move involving a simple volume change from $V$ to $V + \Delta V$, where $\Delta V$ is randomly chosen from $[-\Delta V_{max},\Delta V_{max}]$:

$$\left(o\rightarrow n\right) = \min\left(1,\exp\left\{ -\beta\left... ...t(V-V^\prime\right)-N\beta^{-1}\ln\left(V^\prime/V\right)\right]\right\}\right)$$ (172)

We can also consider trial move that changes the logarithm of the box size from $\ln V$ to $\ln V + \Delta\left(\ln V\right)$. In this case, the integral of $V^N$ over $dV$ is re-expressed as an integral of $V^{N+1}$ over $d\ln V$, and the acceptance rule is the same as the one above except for a factor of $(N+1)$ multiplying $\beta^{-1}$, instead of $N$.

The C-code mclj_npt.c implements an NPT MC simulation of the Lennard-Jones liquid using both particle displacements and log-$V$ displacements. For each cycle, there is a $1/(N+1)$ probability that a trial move is a volume displacement. The trial move generates a random displacement, computes a new box length, rescales all particle positions, scales the cutoff radius, and recomputes the tail corrections and shift, if applicable. If the Metropolis criterion is not met after a random number is selected, then all of these operations are undone. Otherwise, the new box size with the newly scaled particle positions is kept. The particle displacement algorithm is the same as in mclj.c.

As an exercise, you can use the code to regenerate Figure 5.3 in the text, which is again a slice through the phase diagram of the Lennard-Jones fluid at $T$ = 2.0. This temperature is above the critical tempeerature, so we do not anticipate a phase transition at the pressures investigated. However, we saw that when we considered $T$ = 0.9 using the NVT MC simulation, negative pressures were predicted, indicating that the system would have liked to phase separate but couldn't due to its fixed density and finite size. That is, at the density specified, there might not be enough particles to “nucleate” the denser of the two phases. NPT simulations in principle offer a way around that by allowing the system density to fluctuate.

I ran the code with $N$ = 108 particles for 10$^6$ cycles (Note that I have changed my definition of “cycle”. Before, one “cycle” was $N$ moves; now it is a single move. This distinction isn't important for now, but I thought you'd like to be made aware.) The log-volume maximum displacement was set at 0.25, and the maximum particle displacement varied from 0.3 for $P$, to 0.5 at the lowest value of $P$. You can see from Fig. 20 that the data at $T$ = 2.0 is equally well reproduced here as it was using conventional NVT MC (Fig. 13). However, for $T$ = 0.9, we notice that the densities which arise are clearly indicate a high-density phase is prevalent. (Indeed, we saw in NVT simulations that forcing a $T$=0.9 system to exist at densities below about 0.75 resulted in negative pressures!) This code also computes the pressure from the virial, and the measured pressure and imposed pressures agreed, as you can see from the right-hand panel in Fig. 20.

Figure 20: (Left) Pressure vs. density in a Lennard-Jones fluid at two temperatures computed using NPT MC simulation of systems of $N$ = 108 particles. Each point is the average of three independent simulations, all initialized at a density of 0.5. (Right) Measured pressure vs. requested pressure for all simulations.

For temperatures near the critical temperature, we would expect the fluctuations in density to be maximum. As an exercise, you can modify mclj_npt.c to compute the average fluctuations in $\rho$.