Spatial updating in the great grand canonical ensemble Free

In spatial updating grand canonical Monte Carlo algorithms,^1,2 the elementary steps, insertions and removals, are deduced by examining the local environment around a point that is either selected at random or according to a predefined sequence. Spatial updating is more efficient than standard grand canonical updating³ and is ideal for parallel implementation via domain decomposition.² A topic of considerable theoretical and practical importance is simulation of phase transitions involving very dense, nearly-incompressible phases, i.e., fluid-solid phase transitions. In protein crystallization, for instance, it is believed that fluid-fluid demixing is metastable against fluid-solid separation.⁴ Grand canonical simulations of dense phases are difficult because particle transfers have low acceptance probabilities.

Consider spatial updating grand canonical Monte Carlo for a fluid of hard spheres of diameter $σ$⁠. The volume, $V$⁠, is divided into $4n3$ rhombic dodecahedra [the Wigner–Seitz cell of the face-centered cubic (fcc) lattice] of volume $σ3/2$⁠. Consider a configuration of $N$ particles at positions $r1,r2,…,rN$⁠. A point, $r$⁠, is generated with uniform probability within the dodecahedron. Then, spatial updating is implemented by analyzing the local environment at point $r$ in cell $J$ where $J=1,2,…,4n3$⁠. If $|r−rn|>σ/2,∀n=1,2,…,N$⁠, particle insertion is attempted and the position of the particle is obtained at random within a sphere of diameter $σ$⁠. Otherwise, there is one particle, $m$⁠, for which $|r−rm|≤σ/2$ and $|r−rn|>σ/2,∀n≠m$⁠. In such a case, an attempt is made to remove particle $m$⁠. For the hard sphere fluid, the transition from state $i$ to state $j$ is accepted according to $min{1,(zv)±1}$ where $(+)$ corresponds to insertion, (−) to removal, and $v=πσ3/6$ is the volume of a sphere. The activity, $z$⁠, of the particles is defined as $z=Λ−3exp(βμ)$ where $μ$ is the chemical potential, $β=1/kBT$ is the inverse temperature,⁵ $Λ$ is the thermal De Broglie wavelength, and $kB$ is Boltzmann’s constant.

The results of the grand canonical simulations, shown in Fig. 1, are in excellent agreement with the prediction of the Carnahan and Starling equation.⁶ In addition, the data extend to densities, $ρ=N/V$⁠, much higher than possible with standard grand canonical algorithms. For reduced densities $ρ∗=ρσ3≳0.94$⁠, the hard-sphere fluid undergoes a first-order phase transition to an fcc crystal.⁷ The pressure and densities of the two coexisting phases are⁷ $βpσ3≅11.7$⁠, $ρfluidσ3≅0.94$⁠, and $ρfccσ3≅1.04$⁠. The solid phase is an expanded fcc crystal in which the equilibrium distance between neighboring particles is slightly greater than $σ$⁠. In the simulations shown in Fig. 1, the distance between the centers of neighboring dodecahedra is $σ$⁠; thus, the stable solid phase will never be probed.

Suppose that the volume is also allowed to vary in a grand canonical system. In so doing, one forms the so-called great grand canonical ensemble⁸ with partition function, $Z$⁠, defined as

In Eq. (1), $p$ is the pressure, $E$ is the energy, and $W=W(N,V,E)$ is the degeneracy or the density of states or the microcanonical partition function.⁹ In view of the connection of $W$ with entropy, $S=kBlnW$⁠, it can be trivially shown that $Z≡1$⁠. Thus, the corresponding thermodynamic potential is zero. Moreover, the fluctuations of $N$⁠, $V$⁠, and $E$⁠, are unbounded and Monte Carlo simulations of a system with partition function given by Eq. (1) cannot be implemented.

In order to bound the range of fluctuations, a constraint or a restriction must be placed on the previous system. For example, it may be specified that $N$ varies within a lower and an upper bound, $Nmin$ and $Nmax$ respectively. Alternatively, upper and lower bounds may be imposed on the system volume and/or the energy. The choice of the constraint is not unique and other possibilities, such as placing high free energy barriers, also exist. The resulting system corresponds to a restricted or a constrained great grand canonical ensemble in which all intensive variables are fixed and all extensive variables fluctuate within specified bounds.

Even with the imposition of a constraint, it is still not obvious as to how this ensemble can be constructed and simulated since only two of the three intensive variables $μ$⁠, $p$⁠, and $T$⁠, can be specified independently. Suppose that $p$ and $T$ are specified and a constraint is imposed on $N$⁠, i.e., $Nmin≤N≤Nmax$⁠. The partition function of $N$ particles at $p$ and $T$⁠, $Δ(N,p,T)$⁠, is defined as

where

is the canonical partition function of a system of $N$ particles at volume $V$ and temperature $T$⁠. These isobaric systems must be weighted with a Boltzmann type of factor and summed to yield the desired ensemble. Assuming that the Boltzmann factor is of the type $exp(βGN)$⁠, the partition function, $Z∗$⁠, of the constrained ensemble is

The set of $GN$ plays the role of $N$-dependent weights that must somehow be specified and “$∗$” signifies the presence of a constraint. It is easily checked that the assignment

i.e., the Gibbs function for $N$ particles at $p$ and $T$ will result in uniform sampling of all $N$ within $[Nmin,Nmax]$⁠. $μN$ is the $N$-dependent value of the chemical potential.

Alternatively, one may view the restricted great grand canonical ensemble as a series of grand canonical systems at the same chemical potential, $μ$⁠, and temperature, $T$⁠, but at different volumes, i.e.,

with

In Eq. (6), $pV$ are $V$-dependent weights and $Ξ(μ,V,T)$ is the grand partition of an open system at volume $V$⁠, and temperature $T$⁠. The assignment

i.e., the pressure, $pV$⁠, of an open system at $V$ and $T$ will again result in uniform sampling.

Consider simulations of a system with partition function given by Eq. (4). Inspection of Eq. (4) indicates that three types of moves are possible: particle displacements, volume changes, and particle transfers (i.e., insertions and removals). Displacements and volume changes are accepted with ordinary criteria.¹⁰ Using spatial updating, the transfer step is accepted according to

where $(+)$ corresponds to insertion and (−) to removal. In Eq. (9), $ΔE$ is the energy difference between state $i$ with $N$ particles and state $j$ with $N±1$ particles, $ΔG=GN±1−GN$⁠, and $v$ is the volume of a sphere. Transitions that violate the constraint in Eq. (4) are rejected. The implementation of the previous simulation scheme requires specification of the temperature, $T$⁠, the pressure, $p$⁠, the type of the constraint, and a set of free energies, $GN$⁠. The outcome of the simulation is the joint probability of $N$⁠, $V$⁠, and $E$⁠, $P(N,V,E)$⁠, defined as

and obtained as a histogram of observations. If $βGN=−lnΔ$⁠, it follows from Eq. (10) that $P(N,V,E)=const$ within the $(N,V,E)$ range permitted by the constraint. From Eq. (10) it also follows that measurement of $P(N,V,E)$ yields the density of states $W(N,V,E)$ in terms of all its independent variables.

A suitable set of free energies, $GN$⁠, can be found iteratively. An initial guess is made, a short simulation is executed, $W$ is determined from the measured histogram $P(N,V,E)$⁠, $Q$ and $Δ$ are estimated from Eqs. (3) and (2), and finally a refined estimate for $GN$ is obtained from Eq. (5). The process is then repeated using the refined estimates for $GN$ at each iteration until the resulting histograms are nearly flat. Once satisfactory estimates for $GN$ have been obtained, a long simulation is executed to collect good-quality statistics and obtain $W(N,V,E)$⁠.

Once $W(N,V,E)$ has been determined, one can calculate microcanonical, canonical, grand canonical, and isobaric averages using histogram reweighting.¹¹ For example, for $N$ particles at $V$ and $T$⁠, the canonical partition function, $Q(N,V,T)$⁠, is first found from Eq. (3) and the canonical average for the energy from

A single simulation does not explore the whole range of $N$⁠, $V$⁠, and $E$⁠, and it thus yields relative probabilities $P(N,V,E)$⁠. Since the value of $Z∗$ in Eq. (10) cannot be estimated, $W$ and all other partition functions obtained from $W$ contain the same arbitrary multiplicative constant. The corresponding free energies thus contain the same additive constant. However, since $N$⁠, $V$⁠, and $E$ fluctuate, the variation of the partition functions with respect to those variables can be calculated, allowing for determination of absolute free energies and entropies.¹² For example, for $N$ particles at $V$ and $T$⁠, $p$ and $μ$ may be found from $Q(N,V,T)$ via

by approximating the derivatives with difference type of equations. Likewise, the chemical potential of $N$ particles at $p$ and $T$⁠, is found from $Δ(N,p,T)$ via

The previous method was used to find $μ=μ(p)$ for $N=256$ hard sphere particles by imposing the restriction $200≤N≤300$⁠. The volume was divided into 256 dodecahedra and each simulation was executed in sweeps. One sweep comprises 256 elementary steps, each consisting of an attempted particle transfer, a volume change, and a particle displacement. Initial values were assigned to $GN$ and then refined as explained previously. Once a suitable set of $GN$ was determined, a long simulation of $(10–50)×106$ sweeps was executed to collect statistics: see Fig. 2 for a test case. For simulations at different values of $p$⁠, the set of $GN$ was estimated by histogram reweighting¹¹ of the data of the previous simulations.

The results for the hard sphere model, shown in Fig. 3, are in very good agreement with the results of the Carnahan and Starling⁶ equation. For $βpσ3>4$⁠, the range of fluctuations in $N$ was reduced to speed up the calculations $(250≤N≤260)$⁠. The circles correspond to results for $N=256$ hard spheres at fixed $p$⁠. The value of $μ$ was determined by approximating the derivative in Eq. (13) with a centered finite difference equation. Other types of approximations gave nearly identical results. The crosses in Fig. 3 correspond to $N=256$ particles at fixed $V$ and $p$ and $μ$ were estimated from Eq. (12).

A simple way to locate the fluid-solid transition for the hard sphere model is to simulate the solid phase independently and compare the values of $p$ and $μ$ of the fluid with those of the solid. In order to prevent the fluid from solidifying at high densities, an additional constraint was placed on the minimum allowed volume. In this way, the $μ−p$ curve of the fluid can be extended to high densities $(ρσ3≲0.95)$⁠, for which the fluid may be metastable against the solid. Assuming that the solid phase is of the fcc type without vacancies, the properties of the solid may be estimated from a modification of the cell model of Hoover and Ree.⁷ In this model, the simulation volume is divided into $4n3$ dodecahedra. Each dodecahedron may be either empty or occupied by one particle. This model was simulated using the same procedure as the hard sphere fluid and the value of $μ$ that corresponds to full occupancy (i.e., $N=4n3=256$⁠) was determined as a function of $p$⁠: see the inset of Fig. 3. The point at which the two $μ−p$ curves intersect defines coexistence. For $N=256$⁠, the coexisting properties are: $βpσ3≅11.16$⁠, $ln(zσ3)≅15.636$⁠, $ρfluidσ3≅0.932$⁠, and $ρfccσ3≅1.029$⁠.

The present algorithm allows for determination of the density of states, $W(N,V,E)$⁠, in terms of all its independent variables. In addition, since $N$⁠, $V$⁠, and $E$ fluctuate, the variation of the partition functions with respect to those variables may be estimated, allowing for calculation of absolute free energies. The technique described here executes a nearly uniform random walk in $N$⁠, $V$⁠, and $E$⁠. In other techniques of the same type, i.e., flat-histogram^13,14 or tempering^15,16 techniques, one (or more) of the variables $(N,V,E)$ is kept fixed and thus only free energy differences can be calculated. If absolute free energies are desired, a series of simulations toward a state of known free energy (e.g., ideal gas) must be employed.

Shell et al.¹⁷ have proposed a technique that also implements a uniform random walk in $N$⁠, $V$⁠, and $E$⁠. The difference between the present method and that of Shell et al. is the nature of the particle transfer step. In this work, the combination of spatial updating with volume dilations and contractions allows the simulation to explore very high densities and results in enhanced statistical performance. If the instantaneous value of the volume of a given configuration of $N$ particles is greater (less) than the average, then it will rapidly fill with (become depleted of) particles. In the context of fluid-solid and solid-solid transitions, the lattice or phase switch Monte Carlo method of Wilding and Bruce¹⁸ is currently the only method that leads to direct determination of phase equilibria of this type.

The approximation of derivatives with difference equations appears frequently in statistical mechanics. For example, the calculation of $μ$ via test particle insertions¹⁹ is based on this type of approximation. Although there are no inconsistencies in the data shown in Fig. 3, sufficient care must be exercised if precise results are required. In the approximation of a derivative with a centered finite difference, the error of the estimate scales linearly with the square of the spacing of the discrete points. Hence, a linear fit of the estimate versus the square of the spacing will result in a better prediction of the derivative.

The primary factor that determines the length of a simulation is the number of entries in the histogram of $N$⁠, $V$⁠, and $E$⁠, $P(N,V,E)$⁠. Each entry must the visited a sufficient number of times to allow for a precise determination of $W(N,V,E)$ and a consistent calculation of the required derivatives without noise. In this work, it was assumed that the stable solid structure is of the fcc type without vacancies or lattice defects. Future work comprises investigations of the relative stability of other lattice types, such as the hexagonal closed-packed structure, as well as consideration of the effects of vacancies, lattice defects, and impurities. The present technique is ideal in simulations of second-order phase transitions because a single simulation may yield the required information to implement a complete finite-size scaling analysis.²⁰