Accurate and precise determination of critical properties from Gibbs ensemble Monte Carlo simulations

The Lennard-Jones (LJ) particle constitutes the prototypical model for realistic atomic fluids and has been used in the development of sophisticated theories of liquids.^1,2 The accurate determination of the vapor–liquid phase envelope for LJ particles has been the focus of many computational studies using numerous simulation techniques. Early works used simulations at a large number of state points to construct a reversible path either through the two-phase region by applying constraints³ or to a state of known free energy.⁴ Another approach also uses numerous simulations in the NVT or NpT ensemble and test particle insertions to find state points that satisfy the conditions for phase coexistence.^5–8 Most commonly, Monte Carlo simulations with fluctuating particle numbers in the Gibbs ensemble^9–12 or the grand-canonical ensemble^13–15 are used for the determination of fluid phase equilibria. All of the above techniques involve simulations where the periodically replicated system is homogeneous, but increases in computational resources have also led to the use of molecular dynamics simulations for systems with explicit vapor–liquid interfaces.^16,17

For approaches utilizing ensembles with constant volume, one can perform simulations for multiple system volumes and use rigorous mixed-field finite-size scaling (FSS) approaches¹⁸ to extrapolate the critical properties to the thermodynamic limit. Potoff and Panagiotopoulos¹⁵ performed a FSS study using grand-canonical Monte Carlo (GCMC) simulations with standard cubic simulation boxes for the LJ fluid with tail corrections for system sizes smaller than about 5000 particles and determined the critical properties as T_c = 1.3120(7), ρ_c = 0.316(1), and p_c = 0.1279(6), where the numbers in parentheses give the statistical uncertainties in the last digit reported by the authors; in this case, these are standard deviations obtained by performing three independent sets of simulations with different random number seeds. Caillol¹⁴ performed GCMC simulations for system sizes smaller than about 2000 particles for the full LJ potential on the surface of a four-dimensional hypersphere, a technique that does not require long-range corrections, and obtained the critical properties by a FSS analysis as T_c = 1.326(2) and ρ_c = 0.316(2), where the statistical uncertainties are the standard deviations obtained by a block analysis. Okumura and Yonezawa⁸ used FSS and molecular dynamics simulations with the NVT test particle insertion technique for the LJ fluid with tail corrections for system sizes smaller than 2000 particles; they estimated T_c = 1.3207(4), ρ_c = 0.316(1), and p_c = 0.1288(5), where the statistical uncertainties are the standard deviations reported by the authors.

The Gibbs ensemble Monte Carlo (GEMC) technique pioneered by Panagiotopoulos^9,10 is another appropriate method to determine the vapor–liquid coexistence curve (VLCC) and critical properties. This technique utilizes two distinct simulation boxes which are coupled via Monte Carlo moves that allow for particle transfers and volume exchanges. The acceptance rules are chosen such that the subsystems for liquid and vapor phases are in thermodynamic equilibrium, i.e., they have the same chemical potential, pressure, and temperature. The GEMC method allows for fluctuations of energy, volume, and the number of particles in each phase. Since density fluctuations diverge near the critical point, GEMC simulations cannot be used to directly observe the critical point. Furthermore, near the critical point, the phase identity may change, i.e., throughout the simulation trajectory, a given simulation box can represent the vapor phase during certain periods and the liquid phase during others. As a consequence, the orthobaric densities can only be obtained from density probability histograms.¹¹ Close to the critical point, Smit et al.¹¹ observed a third peak in the density histogram that was tentatively explained in terms of the diminishing interfacial free energy near the critical point that would allow for the short-lived formation of liquid droplets or vapor bubbles or even double layers. Using the scaling law¹⁹ and the law of rectilinear diameters²⁰ to determine the critical properties,²¹ Smit²² obtained T_c = 1.316(6) and ρ_c = 0.304(6) for the LJ potential with r_cut = 2.5σ and long-range corrections for system sizes smaller than about 500. In 2013, Mick et al.²³ carried out GEMC simulations for much larger system sizes and found T_c = 1.312(2) and ρ_c = 0.316(3) for a system containing 10 000 LJ particles.

Given the popularity of the GEMC technique, but also motivated by the disagreement in the critical properties determined from various FSS studies,^8,14,15 we perform an extensive set of GEMC simulations to assess the influences of system size (800, 1800, and 5500 LJ particles), of potential truncation (r_cut = 3.5, 5, and 8σ with analytical tail corrections and a shifted LJ potential with r_cut = 8σ), and overall system density (0.296 ≤ ρ_t ≤ 0.336). We also explore the effect of two different ways to treat volume moves that result in a violation of the cutoff distance (i.e., L < 2r_cut), which frequently occur in simulations for small systems (N = 216). In addition, we present some data for n-decane to address whether the conclusions also hold for more complex articulated particles.

The Gibbs ensemble was originally introduced by Panagiotopoulos⁹ as a combination of the canonical, isobaric isothermal, and grand-canonical ensembles. Smit et al.¹¹ took a different point of view and considered it as a new ensemble. They derived the Gibbs ensemble partition function and demonstrated that, in the thermodynamic limit, it is equal to the canonical partition function when the system volume, V_t, is held constant. The Gibbs ensemble partition function considers the distributions of N_t particles and V_t over the two subsystems and all possible configurations as follows (for a one-component system):

where N_I and V_I denote the number of particles and volume of box I, respectively; r_NI and r_Nt−NI are the positions of the particles in boxes I and II, respectively; U_I and U_II are the potential energies of boxes I and II, respectively; and Λ and β are the thermal de Broglie wavelength and inverse temperature, respectively.

From the free energy density for the Gibbs ensemble and accounting for surface effects, one can prove that for an infinite system in the two-phase region, one box must contain the liquid phase and the other the gas phase.¹¹ Nevertheless, for a finite system, one needs to consider the surface contributions to the free energy in a different way to obtain the probability of finding the density distribution over the subsystems.¹¹ 

Extensive GEMC simulations are carried out in order to compute the orthobaric vapor and liquid densities and the saturated vapor pressures for the LJ fluid. The data are reported in reduced units. Unlike previous GEMC simulation studies, the present simulations are performed only at temperatures very close to the critical point in order to reduce any spurious effects that may originate from using the scaling and rectilinear laws over a temperature range beyond the validity applicable to their truncated Taylor series expansions. Although this reduces any systematic error, the tradeoff is that large fluctuations in the near-critical region require very long simulations to reduce the statistical uncertainties. For LJ particles, T = 1.27, 1.28, 1.29, 1.295, 1.30, and 1.305 are considered. The equilibration and production periods for each independent simulation consist of (0.6 to 1.2) × 10⁵ and 9.6 × 10⁵ Monte Carlo cycles (MCCs), where each MCC consists of N_t randomly selected trial moves. The equilibration period is divided into three parts: (i) starting from an initial configuration with equal volume for both boxes, but 80% of the particles in box I, the system is pre-equilibrated for 6000 MCCs using only translational moves, (ii) another period of 12 000 MCCs using only translational and volume exchange moves, and (iii) for the remainder of the equilibration period, translational, volume exchange, and particle swap moves are all used. For the production period, the fraction of volume moves is set to 2/N with a maximum displacement selected so that about 50% of the moves are accepted (i.e., about one volume move per MCC is accepted).²⁶ The fraction of swap moves is set at each temperature so that between 1 and 3 particle transfers are accepted per MCC,²⁶ and the multi-range configurational-bias Monte Carlo (CBMC) algorithm²⁴ with a hard cutoff at 0.5σ (to define an overlap), an inner cutoff at 1.5σ for the calculation of the Rosenbluth weights, and 10 trial sites is used. Translational moves of randomly selected particles make up the remaining fraction of the moves. The pressure is computed via the virial formula after every 5 MCCs.

Our most expensive (and accurate) set of GEMC simulations utilizes a system size consisting of 5500 particles, a spherical potential truncation at r_cut = 8σ with analytical tail corrections for energy and pressure,²⁵ and an overall density of ρ_t = N_t/V_t = 0.316. Since the lengths of boxes I and II (with L_I = V_I^1/3 and L_II = (V_t − V_I)^1/3) fluctuate in a GEMC simulation, the common practice is to use a fixed r_cut value that is somewhat smaller than half of the average box length for the liquid phase (often achieved by adjusting the system size so that it is compatible with the r_cut value tied to a specific force field). Recently, Cortes et al.²⁶ explored different GEMC simulation protocols and recommended that, for temperatures below 0.9T_c, V_t should be adjusted so that about 10%-20% of the particles are found in the vapor phase; this setup leads to a very large volume for the vapor phase at temperature far from T_c, and different r_cut values should be used for the liquid and vapor phases when the Ewald summation technique is used for the Coulomb interactions. The influence of truncation distance on the VLCC and critical properties is investigated by making use of various cutoff distances, namely, 3.5σ, 5σ, and 8σ for N_t = 5500 and ρ_t = 0.316. In addition, we perform GEMC simulations using a shifted potential to assess whether r_cut = 8σ is sufficient to obtain accurate vapor liquid equilibrium (VLE) data. To investigate the effect of overall density on the VLCC and critical properties, GEMC simulations are performed with ρ_t = 0.296, 0.316, and 0.336 for N_t = 5500 and r_cut = 5σ. Furthermore, in order to assess the effect of system size, GEMC simulations for N_t = 5500, 1800, and 800 for r_cut = 3.5σ and ρ_t = 0.316 are used and also for an even smaller system containing only 216 particles with r_cut = 2.5σ and ρ_t = 0.316.

For comparison of density probability distributions (DPDs), single-phase NVT simulations are performed at ρ = 0.316 and T = 1.30, 1.32, and 1.36 for N = 5500 and r_cut = 5σ. The production periods of these NVT simulations consist of 2.4 × 10⁵ MCCs.

In addition to this large set of simulations for the LJ fluid, GEMC simulations are also carried out to compute VLCC and critical properties for hard-core square-well (HCSW) particles. It should be noted here that simulations for models including a hard-core potential in ensembles with fluctuating volume are particularly challenging because of the inefficiency of volume moves. The acceptance probability of trial moves attempting a decrease in volume (note due to the coupling of volume changes in the Gibbs ensemble, every volume move involves a decrease of the volume of one box) tends toward zero as the system size increases because the separation of the closest pair limits the extent of the acceptable volume displacement. Four different ratios of the square-well to hard-core truncation distances are explored (λ/σ = 1.25, 1.50, 1.75, and 2.00) with overall densities (ρ_t = 0.398, 0.307, 0.264, and 0.257, respectively) falling close to the critical densities reported by other authors.^27–30 GEMC simulations are performed for systems containing 400 and 800 HCSW particles. Due to the faster computation of energy changes, the simulations for the HCSW systems are extended to 2 × 10⁵ and 2 × 10⁶ MCCs for equilibration and production with the distribution of moves following the recipe described above for the LJ systems. The pressure is not computed for the HCSW particles because the usual virial formula does not hold.

To explore whether the observations made for LJ and HCSW particles also hold for molecules with articulated structure, coupled-decoupled CBMC simulations³¹ in the Gibbs ensemble are also carried out for n-decane molecules described by the TraPPE–UA force field.³² Two system sizes (N_t = 200 and 1600) and, only for the larger system, four potential truncations schemes (r_cut = 14, 20, and 30 Å with analytical tail corrections and a shifted potential with r_cut = 30 Å) are explored. In these simulations, the total system volume is adjusted at each temperature so that about 20% of the molecules are found in the vapor phase.²⁶ The equilibration and production periods of these GEMC simulations for n-decane consist of 2 × 10⁴ MCCs for both system sizes and of 2 × 10⁵ and 8 × 10⁴ MCCs for the 200- and 1600-molecule systems, respectively. The fraction of volume and swap moves is adjusted to yield about one accepted move of each type per MCC, and the remainder is distributed evenly over translational, rotational, and CBMC conformational moves.

To estimate the statistical uncertainties in the coexistence properties, eight independent simulations are performed for all systems and state points considered in this work. The data from these independent simulations are used to estimate the standard error of the mean (SEM), and these are multiplied by 2.365 to obtain the 95% confidence intervals reported for all numerical values in tables and figures. The procedure used for the determination of the critical properties and the associated error propagation is described in the  Appendix.

One practical difficulty that inevitably arises for the analysis of GEMC simulations in the near-critical region is the so-called box-identity switch where a given simulation box may contain the lower-density gas phase for part of the trajectory and the higher-density liquid phase for another part of the trajectory with the remaining part of the trajectory spent transforming from low to high density and vice versa. In regions of N_t, V_t, T space where such box-identity switches occur, simply taking the running averages of the density for each of the boxes would lead to systematically incorrect answers (that in the limit of infinite simulation length would converge to ρ_t because there is no gravity that pins the liquid phase to the “lower” box).

As suggested by Smit et al.,¹¹ a more appropriate way to analyze the GEMC results is to use DPDs and ascribe the peaks to the densities of the two coexisting phases. In this work, we use the combined DPDs obtained from N_I/V_I and N_II/V_II data at the end of every MCC (total of 1.92 × 10⁶ data points) to determine the vapor and liquid coexistence densities for each independent simulation. These are then averaged to calculate the sample mean and the corresponding standard error of the mean. As a further reliability test, we also determine the cumulative DPD by combining the data from the eight independent runs (total of 1.536 × 10⁷ data points) and compare to the sample mean computed from the independent simulations. To determine the vapor and liquid densities from a given DPD, we fit a Gaussian peak to only the data above 0.75 of maximum height for the vapor and liquid peaks. This procedure removes any spurious effects from the non-Gaussian parts of the DPD that become prevalent very close to the critical point. We also assessed that using data above 0.5 of maximum height does not yield any significant change in the coexistence densities (see also Section IV B).

In Fig. 1, the results from this DPD-based analysis method are demonstrated for the simulations with N_t = 5500, ρ_t = 0.316, and r_cut = 8σ with tail corrections. It is clearly evident that at T = 1.30, the DPDs for the independent simulations still show considerable differences despite that each independent simulation consisted of 9.6 × 10⁵ MCCs. The corresponding sample means and 95% confidence intervals are ρ_vap = 0.2084 ± 0.0078 and ρ_liq = 0.4303 ± 0.0040, i.e., the relative uncertainties are 3.7% and 0.9%, respectively. The relative uncertainties decrease to 1.3% and 0.2% at T = 1.27 but are 6.4% and 1.6%, respectively, at T = 1.305. The numerical data for all systems and state points are provided in the supplementary material.⁴⁷ The DPDs in Fig. 1 also indicate zero probability for densities near 0.316 for the cumulative distributions at T ≤ 1.30, i.e., box-identity switches did not occur in any of the independent simulations over this temperature range for N_t = 5500. In contrast, at T = 1.305, the vapor and liquid regions are connected by non-zero probabilities.

In the DPD analysis presented above, the data points obtained for each simulation box are combined to determine the DPD that is then fitted to a Gaussian. There are two other options for the computation of the DPD, termed here as separated and sorted that both yield misleading answers as will be shown below. For the separated DPD, the instantaneous density values are collected separately for boxes I and II from a GEMC simulation and for the left and right halves of a box from a canonical ensemble simulation. For the sorted DPD, each pair of instantaneous density values obtained from the two GEMC boxes or the two halves of the canonical box at each MCC are sorted into the low (vapor) and high (liquid) DPD. Since the number of data points is equal for all three options, it is obvious that summing the two separated or the two sorted DPDs yields the combined DPD.

Fig. 2 illustrates the three types of DPDs obtained from GEMC and canonical ensemble simulations at T = 1.30 and 1.36 (below and above T_c, respectively) for N_t = 5500 and ρ_t = 0.316 (using cumulative data from all eight independent simulations). As can be seen at T = 1.30, the separated and sorted DPDs for the GEMC simulations are identical and coincide with part of the combined DPD because no box-identity switches occur in any of the eight simulations. Thus, the Gaussian peak analysis would in these cases lead to consistent answers. However, the peaks are somewhat asymmetric with a more extended tail toward intermediate densities. Thus, straightforward computation of the mean gives vapor and liquid densities that are larger by 2.9% and smaller by 1.2%, i.e., similar in magnitude as the 95% confidence intervals obtained from the independent simulations. The peak asymmetry is a finite-system effect,^11,13,18 and we recommend to use the Gaussian peak fits to the combined DPDs. For the canonical ensemble at T = 1.30, the fallacies of the separated and sorted analyses are exposed. Since the simulations are not of infinite length, the separated DPD for the left and right halves does not match and yield different Gaussian fit centroids and mean values that differ by 1.1%. Clearly, the sorted DPDs must be different; individually, they do not resemble a Gaussian peak and their mean values differ by 7%. The combined DPD is symmetric by construction (because the left and right halves have the same volume, a condition that does not hold for the two boxes in the GEMC simulations); its shape is Gaussian and the centroid is at ρ_t = 0.316 because of the maximum in the number of accessible states at ρ_t.

Qualitatively similar behavior as for the canonical ensemble at T = 1.30 is observed for the GEMC simulations and canonical ensemble simulations at T = 1.36 (see Fig. 2). The combined DPDs yield a single Gaussian-type peak that is somewhat broader for the GEMC simulations, but also nearly symmetric. The separated DPDs do not match and yield mean values that differ by 0.6% and 0.3% for the GEMC and canonical ensemble simulations, respectively, but it is important to note that above T_c, the distributions for boxes I and II become very similar (and identical in the limit of infinite simulation length). Of course, the differences are magnified for the sorted DPDs that yield low and high mean densities that differ by 9% and 5% for GEMC and canonical ensemble simulations, respectively. Based on the data shown, it is clear that a sorted (one-dimensional) DPD yields very non-Gaussian peaks and density differences even when the system is far above the critical point. This provides an explanation for the density cusp observed for the phase dome by Valleau³³ from actual GEMC simulations and from analysis of temperature-and-density-scaling Monte Carlo simulations. Valleau³³ utilizes a two-dimensional DPD that involves sorting each pair of data points into low and high density. Although this two-dimensional DPD yields a single Gaussian-type peak, its centroid reflects the differences observed here for the sorted one-dimensional DPDs. The differences in these sorted mean densities will decrease with increasing system size at near-critical temperatures because the DPDs become narrower due to decreasing fluctuations. Importantly, the differences in the sorted DPD densities persist at temperatures well above T_c. However, the differences between the sorted, separated, and combined DPDs vanish for GEMC simulations at temperatures well below T_c where the probability for box-identity switches becomes negligible. Thus, the cusp observed by Valleau is caused by a problem with the data analysis and not by an inherent problem of GEMC simulations. In Section IV E, we will show that allowing one box to vanish will also lead to formation of a third peak in the DPD.

Here, we present first the results for our most expensive set of GEMC simulations carried out for N_t = 5500, ρ_t = 0.316, and r_cut = 8σ with tail corrections that should also be the most accurate data. Using the sample mean densities obtained from the DPD analysis with Gaussian fits described in Sec. IV A (see Fig. 1), the orthobaric vapor and liquid densities are obtained at six temperatures falling in the range of 1.27 ≤ T ≤ 1.305. These temperatures are sufficiently close to T_c that extrapolation to the critical point using the scaling law and the law of rectilinear diameter^19,20 should hold very well. The error propagation procedure is described in the  Appendix. For this system, we find T_c = 1.3128 ± 0.0016 and ρ_c = 0.316 ± 0.004. Once T_c is known, the critical pressure is obtained by an Antoine fit³⁵ to the saturated vapor pressures, and we find p_c = 0.1274 ± 0.0013. These results are in excellent agreement with those obtained by Potoff and Panagiotopoulos¹⁵ from a mixed-field finite-size scaling analysis of grand-canonical ensemble data with histogram reweighting methods. Numerical data are listed in Table I. In fact, the percentage differences in T_c, ρ_c, and p_c are only about 0.06%, 0.0%, and 0.4%, respectively. It should be noted that the statistical uncertainties provided by Potoff and Panagiotopoulos¹⁵ are standard deviations from three independent simulations, and conversion to the 95% confidence interval of the sample mean would increase them by a factor of 2.5. Thus, these two simulation estimates of the critical properties of the LJ fluid have similar precision, but it needs to be emphasized that the current GEMC simulation data are for a finite system with N_t = 5500.

In terms of a sensitivity analysis, the critical properties are also obtained using the same VLCC data but different critical exponents β for the scaling law. We recognize a value of 0.326 as the most accurate estimate of the critical exponent for the Ising universality class.^36,37 However, many earlier GEMC studies used β = 0.32. Critical properties calculated for β = 0.320, 0.326, and 0.333 are presented in Table II. Using VLCC data over our admittedly very narrow temperature range, the differences in T_c and p_c are about a factor of two smaller than the 95% confidence intervals, but there is obviously a trend of systematically increasing T_c and p_c with increasing β. Since the line describing the law of rectilinear diameters is nearly vertical, the differences in ρ_c values obtained for this range of β values are much smaller than the 95% confidence intervals. On the other hand, when one tries to fit the VLCC data using the mean-field critical exponent of β = 0.5, then T_c and p_c increase by 1.0% and 5.6%, respectively, and ρ_c decreases by 0.6%. Although these differences may not appear very large, the quality of the fits decreases significantly as is evident from the larger 95% confidence intervals. In addition, we also explore the sensitivity of the critical properties to the range of values considered for the Gaussian peak fits. When selecting data above 0.5 of the maximum heights of the DPDs, then we obtain (using β = 0.326) the following values: T_c = 1.3124 ± 0.0016, ρ_c = 0.316 ± 0.004, and p_c = 0.1272 ± 0.0013, i.e., differences to those using data only above 0.75 of maximum peak height that are at least four times smaller than the 95% confidence intervals.

Numerous other studies have investigated the influence of truncation effects on the VLE properties of the pure LJ fluid and other fluids. Smit et al.²² showed that truncation and shifting of the potential have a considerable effect on the location of the critical point. In recent studies, Fern et al.¹⁷ and Patel et al.³⁸ made use of the assumption that r_cut = 6σ is sufficiently large to eliminate effects of the truncated pair potential on the phase dome loci. To our knowledge, there are no previous studies systematically exploring the effect of r_cut when tail corrections are used. In molecular simulations of homogeneous systems with the full LJ potential, the attractive and repulsive parts are explicitly calculated up to the cutoff distance and the contributions beyond r_cut are commonly approximated by analytical tail corrections assuming that the pair distribution functions have converged to unity beyond the cutoff distance.²⁵ As r_cut becomes larger, more and more interactions are calculated explicitly and the contribution of the analytical tail correction to the potential energy becomes smaller. Therefore, more accurate results are obtained for larger r_cut at the cost of increasing computational load.

Fig. 3 illustrates the cumulative DPDs for various approaches at T = 1.30. For the large system (N_t = 5500 and ρ_t = 0.316), data are also presented for r_cut values of 5 and 3.5σ with tail corrections in addition to the r_cut = 8σ data introduced in Sec. IV B. The DPDs for the two larger r_cut values match extremely well, but the DPD for the liquid phase with r_cut = 3.5σ is shifted slightly to larger densities (still well within the 95% confidence intervals). The corresponding VLCC is depicted in Fig. 4 (top); and it is evident that the data for r_cut = 8 and 5σ match extremely well over the entire temperature range investigated here, whereas the phase dome is slightly wider for r_cut = 3.5σ. On the other hand, the saturated vapor pressures for these three simulations with tail corrections are in close agreement (see Fig. 5 (top)). With regards to the critical properties (see Table I), using r_cut = 3.5σ yields increases of 0.2%, 0.3%, and 1.4% for T_c, ρ_c, and p_c, respectively, compared to the data for r_cut = 8σ. Although their 95% confidence intervals overlap, the differences in T_c and p_c are about 1.5 times larger than the individual 95% confidence intervals. Thus, if very accurate data are essential, then it appears advisable to use an r_cut value of at least 5σ with tail corrections.

When the system exhibits large-scale inhomogeneities, then the analytical tail corrections should not be applied. In single-box, two-phase molecular dynamics simulations,^16,17,38 the density is by construction inhomogeneous and the vapor– liquid interfaces may not be planar. There are methods available in the literature to add long-range corrections in an inhomogeneous system but these methods are applicable only to a planar interface.³⁹ Based on extensive two-phase molecular dynamics simulations, Fern et al.¹⁷ and Patel et al.³⁸ concluded that r_cut ≥ 6σ is sufficient to eliminate truncation effects. Here, we also carry out GEMC simulations for the large system (N_t = 5500 and ρ_t = 0.316) using a LJ potential that is truncated and shifted at 8σ. Shifting of the potential is standard practice in molecular dynamics simulations²⁵ because it yields an energy that is continuous for all distances and, hence, allows for better energy conservation along the trajectory (but we do not use a shifted-force potential that would tilt the shifted LJ potential to yield a zero force at r_cut). From Figs. 4 (top) and 5 (top), it is immediately obvious that the simulations for the shifted potential yield data that systematically and to a large extent deviate from all the other simulations that employ tail corrections. For the shifted potential, the phase dome is narrower, the saturated vapor pressures are elevated by about 2.5%, and T_c is lowered by about 0.5%, well outside the combined 95% confidence intervals. It turns out that GEMC simulations using N_t = 5500 and a shifted potential with r_cut = 8σ yield less accurate results than those using N = 216 and r_cut = 2.5σ but with tail corrections, where the former are more expensive by a factor of about 500.

Valleau³³ also investigated whether the overall density, ρ_t, would influence the VLCC properties obtained by GEMC simulations, but his analysis is based on data from temperature-and-density-scaling Monte Carlo simulations and not from actual GEMC simulations. Valleau concluded that the densities of both phases are shifted by selecting a different ρ_t (data given for a system with N_t = 512 and ρ_t = 0.301, 0.320, and 0.341) and that predictions of the critical density are influenced by changes in ρ_t as the temperature is increased. Martin and Siepmann³² and, more recently, Cortes-Morales et al.²⁶ advocated to adjust ρ_t in a manner that about 10%-20% of the particles are found in the vapor phase because this phase ratio will allow for efficient simulations yielding relatively low statistical uncertainties for the saturated vapor pressure at temperatures far below T_c while not requiring extremely large systems. In the near-critical region, however, the phase dome narrows and finding sufficient particles in the vapor phase (to yield small uncertainties for the saturated vapor pressure) ceases to be an issue.

To investigate the effects of the total system density on the critical properties, we compare GEMC simulation data for ρ_t = 0.296, 0.316, and 0.336 using N_t = 5500 and r_cut = 5σ with tail corrections, i.e., a similar range of ρ_t values but a larger N_t than investigated by Valleau. Different choices of ρ_t affect the average numbers of particles in each of the simulation boxes as can be deduced from the lever rule. Here, the vapor phases at T = 1.27 contain about 21%, 26%, and 31% of the particles for ρ_t = 0.336, 0.316, and 0.296, respectively; that is, all fall well above the 10%-20% range needed for the computation of the saturated vapor pressure with good precision. At T = 1.27, ρ_t > 0.40 would be needed to shift the phase ratio so that fewer than 10% of the particles are found in the vapor phase.

The DPD data in Fig. 3 indicate small shifts to higher densities for both the vapor and liquid peaks obtained for the simulations at T = 1.30 with ρ_t = 0.296 compared to those with ρ_t = 0.336. These small shifts are in the opposite direction as those deduced by Valleau who argued that an increase in ρ_t also increases the orthobaric densities. More importantly, the shifts obtained here are not statistically significant and the numerical values are ρ_vap = 0.2099 ± 0.0059 and ρ_liq = 0.4332 ± 0.0048 for ρ_t = 0.296 and ρ_vap = 0.2037 ± 0.0064 and ρ_liq = 0.4296 ± 0.0034 for ρ_t = 0.336 (see also the supplementary material⁴⁷). Figs. 4 (middle) and 5 (middle) illustrate that the excellent agreement for the GEMC simulations with the three different ρ_t values persists over the entire temperature range considered here. The differences in the critical properties determined for the simulations with ρ_t = 0.296 and 0.336 (see Table I) are 0.05%, 0.6%, and 0.3% for T_c, ρ_c, and p_c, respectively, and are at least a factor of two smaller than the individual 95% confidence intervals.

Finite-size effects are the most important concern for the determination of the critical properties from particle-based simulations. From statistical-mechanical arguments, it is clear that fluctuations increase as the system size is decreased. The DPDs for three system sizes from GEMC simulations at T = 1.30 using ρ_t = 0.316 and r_cut = 3.5σ with tail corrections are shown in Fig. 3. The simulations for the largest system with N_t = 5500 yield two Gaussian-type peaks that are both fairly symmetric and the probability to observe an instantaneous density value equal to ρ_t is vanishingly small. When N_t is decreased to 1800, then the peaks broaden considerably and there is a non-negligible probability for densities near ρ_t connecting the two peaks and allowing for box-identity switches. These trends continue when N_t is further reduced to 800; now the peaks are much broader and well connected at intermediate densities with the probability at ρ_t = 0.316 being only three times smaller than the peak heights. However, the orthobaric densities obtained for these three system sizes still yield overlapping 95% confidence intervals. Nevertheless, as indicated by the data in Figs. 4 and 5, the shape of the phase dome appears to be slightly different where the dome for N_t = 5500 exhibits a slightly larger curvature than found for N_t = 800, and the p_sat values are shifted downward with decreasing system size (they are about 0.5% lower for N_t = 800 than for 5500). The critical temperature for the 5500-particle system is 0.1% higher than for the 800-particle system, a difference that is just smaller than the individual 95% confidence intervals. Concommitantly, the critical pressure for the former system is larger by 1.3% than for the latter system. GEMC simulations are also used to explore a system with N_t = 216, but this requires a reduction of r_cut to 2.5σ. The data for the 216-particle system follow the trends of lower-curvature phase dome, downward shift in p_sat, and decrease in T_c and p_c. However, the effects are more dramatic and it remains an unresolved question whether this is caused by more severe finite-size effects as the system size is further reduced or by the change in the Hamiltonian caused by the different r_cut values.

As mentioned in the Introduction, rigorous mixed-field finite-size scaling approaches (using scaling operators comprised of linear combinations of particle and energy densities) can be applied to simulations in ensembles with fixed volume.^18,37,40,41 For comparison purposes, we plot the apparent critical temperature (for a given system size) as a function of L^−(θ+1)/ν, where L is the box length for simulations in the canonical and grand-canonical ensembles using cubic simulation boxes, and scaling parameters are set to θ = 0.54 and ν = 0.629.^40,42 A comparison of the present GEMC simulation data to the literature data is shown in Fig. 6. It is important to note that the error bars are the ones reported by the respective authors and their values would, in general, increase upon conversion to the 95% confidence intervals used in this work. Analyzing GEMC simulations using DPDs¹¹ and Gaussian peak fits and extrapolation to the critical point via the scaling and rectilinear laws yields T_c values that appear to be much less sensitive to changes in the system size than those obtained from other methods. In earlier GEMC simulations, Recht and Panagiotopoulos⁴³ did not observe significant finite-size effects for the three-dimensional LJ fluid for 200 ≤ N_t ≤ 400, but critical properties were not reported.

Comparatively large system-size dependencies are observed for the studies by Caillol¹⁴ and Okumura and Yonezawa.⁷ The latter used simulations in the canonical ensemble and larger finite-size effects are expected for a closed ensemble. Caillol,¹⁴ Potoff and Panagiotopoulos,¹⁵ and Pérez-Pellitero et al.³⁴ used simulations in the grand-canonical ensemble and the latter two yield consistent finite-size effects, but Caillol’s simulations on the 4-dimensional hypersphere deviate significantly. For the system sizes investigated here, the T_c values obtained from GEMC simulations are always higher than those from grand-canonical ensemble simulations in cubic boxes for a comparable system size, but the difference narrows as the system size increases. The GEMC data for the N = 5500 system using r_cut = 8σ with tail corrections are in extremely good agreement with the thermodynamic limit obtained by Potoff and Panagiotopoulos¹⁵ and Pérez-Pellitero et al.³⁴ using FSS. Overall, as already suggested by Smit et al.,¹¹ the Gibbs ensemble with its fluctuation in the number of particles as well as the volume allows one to obtain critical properties that match well with those for the infinite system.

Data for the apparent critical density and pressure as a function of L^{−(d−1/ν)}, where d is the dimensionality of the system, are illustrated in Fig. 6. For all simulation techniques, the relative statistical uncertainties are about an order of magnitude larger for ρ_c and p_c than those for T_c. Thus, with the exception of the canonical ensemble simulations, there are no significant trends. For ρ_c, the data from Potoff and Panagiotopoulos¹⁵ and Pérez-Pellitero et al.³⁴ appear to yield slopes with opposite sign. For p_c, both the data by Potoff and Panagiotopoulos¹⁵ and our GEMC data indicate a weak increase as the system size increases. Again, our most accurate data (N_t = 5500 and r_cut = 8σ with tail corrections) agree very well with the estimates at the thermodynamic limit obtained by Potoff and Panagiotopoulos¹⁵ and Pérez-Pellitero et al.³⁴ 

The extremely long simulations conducted in this work allow us to collect enough data to construct two-dimensional density histograms in N–V space, as shown in Figs. 7 and 8. The warmer colors indicate regions of higher phase-space density and these N–V maps resemble parallelograms because changes in the number of particles and volume of the two subsystems are coupled. Comparing data for system sizes from 5500 to 216 particles at T = 1.305 ≈ 0.994 T_c, it is found that the smaller systems exhibit much larger fluctuations in both particle number and volume (noting the different color scales) and access much more frequently the region of densities intermediate between vapor and liquid phases. In other words, simulations with small system sizes are most prone to the box-identity switch problem mentioned earlier.

Another feature of these N–V maps that is most pronounced for the smallest system is the missing corners along the long diagonal of the parallelograms (see Fig. 7). We view r_cut as an integral part of the Hamiltonian (in the same manner as the number of neighbors considered in the Ising model is not a function of system size) and, hence, in our standard GEMC protocol, we reject any volume move that attempts to decrease the linear dimension of one of the simulation boxes to less than 2r_cut. This constraint limits the volume integral in Eq. (1) to [(2r_cut)³, V − (2r_cut)³]. Since the corners along the long diagonal are associated with states near ρ_t, those regions of the DPD are depressed. Although the data for N_t = 5500 are not affected by the volume constraint, the missing corners emerge for N_t = 1800 and are significant for N_t = 800; this is also evident from the DPD for this system size at T = 1.30 where a depression is found at ρ_t = 0.316 (see Fig. 3).

There are two other approaches to handle fluctuations resulting in one of the boxes becoming very small: (i) a variable cutoff that adjusts to $r cut ′ = V 1 / 3 /2$ when volume moves are attempted or accepted that yield configurations with V_I < (2r_cut)³ or V_I > V_t − (2r_cut)³, and restoring the original r_cut when the volume split between the boxes becomes more balanced; or (ii) allowing such attempts to proceed without adjusting r_cut, but continuing to apply the minimum image convention, which would artificially lower the particle density within the truncation radius and lead to a loss of spherical symmetry for the explicit-interaction shells that are not compensated by the tail corrections. Either of these two approaches leads to a modification of the Hamiltonian. It is apparent that none of these three schemes are perfect, but one must be chosen to enable such simulations, as large volume fluctuations are an intrinsic problem at temperatures sufficiently close to critical point, especially for small system sizes.

The critical properties estimated with these three schemes for N_t = 216, ρ_t = 0.316, and r_cut = 2.5σ with tail corrections are provided in Table I. T_c, ρ_c, and p_c agree well with each other but fall 0.1%, 1.3%, and 4% below those for the most accurate GEMC simulations. Although this agreement may appear comforting, it is only achieved through analysis involving DPDs, Gaussian peak fits, and removal of data points above T_c that yield very non-Gaussian peaks. It is clearly surprising that simulations of such a small system size can already predict T_c and ρ_c with reasonable accuracy. As discussed previously, however, the smaller r_cut and N_t cause some noticeable differences in the shape of the phase envelope and in p_sat. The pathological behavior of these small systems is most evident at temperatures just above T_c. As illustrated in Fig. 8 for simulations with N_t = 216 at T = 1.315, constraining the volume range to maintain the Hamiltonian results in the removal of the corners in the N–V maps that would otherwise contain a significant part of the unconstrained trajectory. As a consequence, densities near ρ_t are very much depressed in the one-dimensional DPD leading to a spurious minimum and a double peak with rather non-Gaussian shape. The scheme with the variable r_cut yields a single peak, but with a strong cusp centered at ρ_t = 0.316. The two-dimensional N–V density distribution indicates that this cusp is not caused by the two subsystems sampling intermediate densities associated with the formation of vapor bubbles in the liquid phase and liquid droplets in the vapor phase,¹¹ but rather from one subsystem vanishing and the other assuming (by default) the total density. The small peak at ρ = 0 and the long tail at very high densities in the DPD are both due to the box with vanishing particle number and volume.

If T_c is not known, some care must be taken to examine DPDs carefully; whenever in doubt (e.g., because the peak shapes cease to be Gaussian), lower-temperature data should be used to provide an estimate of T_c. In any case, it should be emphasized that the problems related to small system size can be alleviated only by using a system size that is sufficient with respect to r_cut, so that the corners in N–V space become irrelevant. This also explains why Mon and Binder⁴⁴ do not observe a cusp in the restricted Gibbs ensemble where volume moves are not included.

The hard-core square-well fluid is another system of wide-spread theoretical interest because it serves as a reference for statistical-mechanical perturbation theories. Simulations of particles with hard-core interactions are troublesome in ensembles with fluctuating volume and, unfortunately, early GEMC studies for the HCSW fluid⁴⁵ turned out to be less accurate than those for the LJ fluid. Here, we focus on the data for the 800-particle system because the differences in the saturated densities and the critical properties obtained from GEMC with 400 and 800 particles are smaller than the statistical uncertainties for the length of the production periods used here (as an example, the ratio of critical temperatures is provided in Table III). The critical temperatures computed from the 800-particle GEMC simulations are found to be slightly higher (by 0.3%–0.7%) than those obtained from FSS studies using the grand canonical ensemble.^27–30 Although this difference appears to be systematic, the lower bound of the 95% confidence intervals is close to or below the unweighted mean value obtained from the FSS studies. The relative uncertainties for the critical density are about an order of magnitude larger and it is not possible to discern any systematic difference. Although the GEMC approach again yields critical properties for a modest system size that appear to be in satisfactory agreement with elaborate FSS studies, we hesitate to recommend GEMC for the HCSW fluid because of the rather large statistical uncertainties that are related to the inefficiency of the sampling of volume fluctuations required to establish mechanical equilibrium.

To assess whether our findings for the LJ and HCSW fluids also hold for a more complex molecule with conformational degrees of freedom, the VLE for n-decane is studied. The r_cut values of 14, 20, and 30 Å correspond to 3.5σ_CH2, 5.1σ_CH2, and 7.6σ_CH2, respectively. The simulation data for n-decane cover a larger range of reduced temperatures, but only the highest four temperatures are used for the determination of the critical properties. The VLCCs and p_sat curves are shown in Fig. 9. The numerical data for the critical properties are summarized in Table IV. Again, VLE data for the larger system (N_t = 1600) and the two larger r_cut values agree extremely well. Reducing r_cut to 14 Å (the standard value for the TraPPE force field) leads to a 2 K (or 0.3%) increase in T_c and an 80 kPa (or 3.5%) increase in p_c; these differences are about equal to the combined 95% confidence intervals. A reduction of the linear system dimension by a factor of two (N_t = 200 with r_cut = 14 Å) yields a system with critical properties that are remarkably similar to those for N_t = 1600 and r_cut = 30 Å, that is, these reflect a similar fortuitous cancellation of finite-size and truncation errors as observed for the LJ fluid (compare data for N_t = 5500 and r_cut = 8σ to those for N_t = 800 and r_cut = 3.5σ). One should also note that T = 580 K is the highest temperature used here for the fitting of the data for the 200-molecule n-decane system. At this temperature (0.94T_c), problems with having to reject moves because the volume of one box becomes too small are not encountered; note that the ratio of L to r_cut is somewhat more favorable for a molecule consisting of ten connected LJ sites than for a one-site LJ fluid with the same number of particles. Not at all unexpected and consistent with the data for the LJ fluid, a shifted potential with r_cut = 30 Å yields a narrower phase dome, a significant increase in p_sat, and a downward shift of T_c by 7 K (or 1.1%) that falls well outside the combined 95% confidence intervals. The predicted T_c, ρ_c, and p_c values with the standard TraPPE r_cut value of 14 Å and averaged over the two system sizes deviate from the recommended experimental critical properties⁴⁶ by 0.2%, 4%, and 11%, respectively. Thus, the accuracy and precision of these GEMC simulations are better than the accuracy of the force field.

In this study, we determine the VLE and critical properties of the three-dimensional LJ fluid by means of GEMC simulations for a system consisting of 5500 particles with a cutoff distance of 8σ together with tail corrections. We process the probability distributions of vapor and liquid densities by fitting Gaussian peaks. This analysis leads to very precise estimates of the critical properties with T_c = 1.3128 ± 0.0016, ρ_c = 0.316 ± 0.004, and p_c = 0.1274 ± 0.0013 that are in excellent agreement with finite-size scaling studies by Potoff and Panagiotopoulos¹⁵ using mixed-field theory and by Pérez-Pellitero et al.³⁴ using the Binder cumulant parameter.

We investigate the influence of cutoff distance, total density, and system size on the VLE and critical properties of the three-dimensional LJ fluid and do not observe any significant changes in these properties for simulations with different cutoff distances (3.5σ ≤ r_cut ≤ 8σ) and total densities (0.296 ≤ ρ_t ≤ 0.336). We also notice that the critical properties obtained for different system sizes (800 ≤ N_t ≤ 5500) vary relatively little, with the differences smaller than the 95% confidence intervals. This indicates that finite-size effects are relatively small in GEMC simulations compared to other techniques. Thus, when resources are not available for a rigorous finite-size scaling study, then GEMC simulations may provide a suitable route to determine fairly accurate critical properties using relatively small system sizes but still yielding results close to the thermodynamic limit. However, this may not hold for all systems.⁴³ 

We demonstrate that sorted density probability distributions yield erroneous results and a cusp in the VLCC that extends above T_c. From analysis of two-dimensional probability histograms in the N–V plane, we attribute the appearance of a third peak near ρ_t at T ≈ T_c to states where one of the subsystems vanishes and highlight some potential problems associated with volume moves in GEMC simulations for small system sizes.

GEMC simulations for hard-core square-well particles and n-decane molecules indicate that the trends observed for the LJ fluid are also exhibited by systems with discontinuous potentials and more complex molecular systems.

Financial support from the National Science Foundation (No. CBET-1159837) is gratefully acknowledged. Part of the computer resources were provided by the Minnesota Supercomputing Institute.

In this section, we briefly describe the fitting procedure that is used to obtain the critical properties. All six temperatures investigated for the LJ fluid (T > 0.96 T_c) and only the four highest temperatures for n-decane (T > 0.87 T_c) are used as input data. In Section A.2 of Smit’s thesis,²¹ the orthobaric densities, ρ_liq and ρ_vap, as a function of T are fitted to the law of rectilinear diameters and the scaling law for the density,^19,20

where the three-dimensional critical exponent of β = 0.326 is used.^36,37 These two equations are transformed into linear forms to obtain

with transformed variables

The statistical uncertainties in the transformed variables are then given by uncertainty propagation,

Since ρ_liq = N_I/V_I and ρ_vap = (N_I − N_I)/(V_t − V_I), the orthobaric densities are not actually statistically independent. In Smit’s thesis,²¹ the uncertainties in $σ ρ liq + ρ vap 2$ and $σ ρ liq − ρ vap 2$ are written in terms of N_I and V_I, while ρ_liq and ρ_vap are more often used when reporting the VLE data. We obtain the upper bound by taking corr(ρ_liq, ρ_vap) to be 1,

Fitting to the transformed linear equations yields a_r, $σ a r 2$⁠, b_r, $σ b r 2$⁠, cov(a_r, b_r), a_s, $σ a s 2$⁠, b_s, $σ b s 2$⁠, and cov(a_s, b_s) that then give

Note that the expression for $σ B 2$ is updated here to eliminate minor typographical mistakes in Ref. 21.