Since the seminal paper by Panagiotopoulos [Mol. Phys. **61**, 813 (1997)], the Gibbs ensemble Monte Carlo (GEMC) method has been the most popular particle-based simulation approach for the computation of vapor–liquid phase equilibria. However, the validity of GEMC simulations in the near-critical region has been questioned because rigorous finite-size scaling approaches cannot be applied to simulations with fluctuating volume. Valleau [Mol. Simul. **29**, 627 (2003)] has argued that GEMC simulations would lead to a spurious overestimation of the critical temperature. More recently, Patel *et al.* [J. Chem. Phys. **134**, 024101 (2011)] opined that the use of analytical tail corrections would be problematic in the near-critical region. To address these issues, we perform extensive GEMC simulations for Lennard-Jones particles in the near-critical region varying the system size, the overall system density, and the cutoff distance. For a system with *N* = 5500 particles, potential truncation at 8*σ* and analytical tail corrections, an extrapolation of GEMC simulation data at temperatures in the range from 1.27 to 1.305 yields *T*_{c} = 1.3128 ± 0.0016, *ρ*_{c} = 0.316 ± 0.004, and *p*_{c} = 0.1274 ± 0.0013 in excellent agreement with the thermodynamic limit determined by Potoff and Panagiotopoulos [J. Chem. Phys. **109**, 10914 (1998)] using grand canonical Monte Carlo simulations and finite-size scaling. Critical properties estimated using GEMC simulations with different overall system densities (0.296 ≤ *ρ*_{t} ≤ 0.336) agree to within the statistical uncertainties. For simulations with tail corrections, data obtained using *r*_{cut} = 3.5*σ* yield *T*_{c} and *p*_{c} that are higher by 0.2% and 1.4% than simulations with *r*_{cut} = 5 and 8*σ* but still with overlapping 95% confidence intervals. In contrast, GEMC simulations with a truncated and shifted potential show that *r*_{cut} = 8*σ* is insufficient to obtain accurate results. Additional GEMC simulations for hard-core square-well particles with various ranges of the attractive well and for *n*-decane molecules represented by the TraPPE force field yield data that support the trends observed for Lennard-Jones particles. The finite-size dependence of the critical properties obtained from GEMC simulations is significantly smaller than those from grand-canonical ensemble simulations. Thus, when resources are not available for a rigorous finite-size scaling study, GEMC simulations provide a straightforward route to determine fairly accurate critical properties using relatively small system sizes.

## I. INTRODUCTION

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^{3} or to a state of known free energy.^{4} 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^{18} to extrapolate the critical properties to the thermodynamic limit. Potoff and Panagiotopoulos^{15} 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^{14} 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^{8} 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.^{11} Close to the critical point, Smit *et al.*^{11} 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^{19} and the law of rectilinear diameters^{20} to determine the critical properties,^{21} Smit^{22} 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.*^{23} 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* < 2*r*_{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.

## II. THEORETICAL BACKGROUND

The Gibbs ensemble was originally introduced by Panagiotopoulos^{9} as a combination of the canonical, isobaric isothermal, and grand-canonical ensembles. Smit *et al.*^{11} 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.^{11} 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.^{11}

## III. SIMULATION DETAILS

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^{5} and 9.6 × 10^{5} 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).^{26} The fraction of swap moves is set at each temperature so that between 1 and 3 particle transfers are accepted per MCC,^{26} and the multi-range configurational-bias Monte Carlo (CBMC) algorithm^{24} 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,^{25} 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.*^{26} explored different GEMC simulation protocols and recommended that, for temperatures below 0.9*T*_{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^{5} 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^{5} and 2 × 10^{6} 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^{31} in the Gibbs ensemble are also carried out for *n*-decane molecules described by the TraPPE–UA force field.^{32} 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.^{26} The equilibration and production periods of these GEMC simulations for *n*-decane consist of 2 × 10^{4} MCCs for both system sizes and of 2 × 10^{5} and 8 × 10^{4} 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.

## IV. RESULTS AND DISCUSSION

### A. Density probability distributions

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.*,^{11} 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^{6} 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^{7} 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^{5} 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.^{47} 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^{33} from actual GEMC simulations and from analysis of temperature-and-density-scaling Monte Carlo simulations. Valleau^{33} 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.

### B. Critical properties for 5500-particle system with 8*σ* cutoff

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^{35} 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^{15} 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^{15} 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.

Method (N_{t}/ρ_{tot}/r_{cut})
. | T_{c}
. | ρ_{c}
. | p_{c}
. |
---|---|---|---|

Potoff & Panagiotopoulos | 1.3120 ± 0.0007 | 0.316 ± 0.001 | 0.1279 ± 0.0018 |

Pérez-Pellitero et al. | 1.3126 ± 0.0006 | 0.3174 ± 0.0006 | N/A |

5500/0.316/8σ | 1.3128 ± 0.0016 | 0.316 ± 0.004 | 0.1274 ± 0.0013 |

5500/0.316/5σ | 1.3132 ± 0.0017 | 0.316 ± 0.004 | 0.1276 ± 0.0013 |

5500/0.316/3.5σ | 1.3154 ± 0.0016 | 0.317 ± 0.003 | 0.1292 ± 0.0015 |

5500/0.316/8σ/SP | 1.3068 ± 0.0011 | 0.316 ± 0.005 | 0.1275 ± 0.0006 |

5500/0.296/5σ | 1.3131 ± 0.0016 | 0.317 ± 0.004 | 0.1278 ± 0.0012 |

5500/0.336/5σ | 1.3138 ± 0.0015 | 0.315 ± 0.004 | 0.1281 ± 0.0011 |

1800/0.316/3.5σ | 1.3149 ± 0.0022 | 0.318 ± 0.005 | 0.1286 ± 0.0016 |

800/0.316/3.5σ | 1.3139 ± 0.0025 | 0.316 ± 0.006 | 0.1275 ± 0.0017 |

216/0.316/2.5σ | 1.3114 ± 0.0016 | 0.311 ± 0.005 | 0.1232 ± 0.0009 |

216/0.316/2.5σ/VC | 1.310 ± 0.003 | 0.311 ± 0.006 | 0.122 ± 0.003 |

216/0.316/2.5σ/MIM | 1.3112 ± 0.0012 | 0.313 ± 0.004 | 0.1231 ± 0.0009 |

Method (N_{t}/ρ_{tot}/r_{cut})
. | T_{c}
. | ρ_{c}
. | p_{c}
. |
---|---|---|---|

Potoff & Panagiotopoulos | 1.3120 ± 0.0007 | 0.316 ± 0.001 | 0.1279 ± 0.0018 |

Pérez-Pellitero et al. | 1.3126 ± 0.0006 | 0.3174 ± 0.0006 | N/A |

5500/0.316/8σ | 1.3128 ± 0.0016 | 0.316 ± 0.004 | 0.1274 ± 0.0013 |

5500/0.316/5σ | 1.3132 ± 0.0017 | 0.316 ± 0.004 | 0.1276 ± 0.0013 |

5500/0.316/3.5σ | 1.3154 ± 0.0016 | 0.317 ± 0.003 | 0.1292 ± 0.0015 |

5500/0.316/8σ/SP | 1.3068 ± 0.0011 | 0.316 ± 0.005 | 0.1275 ± 0.0006 |

5500/0.296/5σ | 1.3131 ± 0.0016 | 0.317 ± 0.004 | 0.1278 ± 0.0012 |

5500/0.336/5σ | 1.3138 ± 0.0015 | 0.315 ± 0.004 | 0.1281 ± 0.0011 |

1800/0.316/3.5σ | 1.3149 ± 0.0022 | 0.318 ± 0.005 | 0.1286 ± 0.0016 |

800/0.316/3.5σ | 1.3139 ± 0.0025 | 0.316 ± 0.006 | 0.1275 ± 0.0017 |

216/0.316/2.5σ | 1.3114 ± 0.0016 | 0.311 ± 0.005 | 0.1232 ± 0.0009 |

216/0.316/2.5σ/VC | 1.310 ± 0.003 | 0.311 ± 0.006 | 0.122 ± 0.003 |

216/0.316/2.5σ/MIM | 1.3112 ± 0.0012 | 0.313 ± 0.004 | 0.1231 ± 0.0009 |

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.

β
. | T_{c}
. | ρ_{c}
. | p_{c}
. |
---|---|---|---|

0.320 | 1.3124 ± 0.0015 | 0.316 ± 0.004 | 0.1272 ± 0.0012 |

0.326 | 1.3128 ± 0.0016 | 0.316 ± 0.004 | 0.1274 ± 0.0013 |

0.333 | 1.3133 ± 0.0016 | 0.316 ± 0.004 | 0.1277 ± 0.0013 |

0.500 | 1.3263 ± 0.0024 | 0.314 ± 0.005 | 0.1345 ± 0.0027 |

β
. | T_{c}
. | ρ_{c}
. | p_{c}
. |
---|---|---|---|

0.320 | 1.3124 ± 0.0015 | 0.316 ± 0.004 | 0.1272 ± 0.0012 |

0.326 | 1.3128 ± 0.0016 | 0.316 ± 0.004 | 0.1274 ± 0.0013 |

0.333 | 1.3133 ± 0.0016 | 0.316 ± 0.004 | 0.1277 ± 0.0013 |

0.500 | 1.3263 ± 0.0024 | 0.314 ± 0.005 | 0.1345 ± 0.0027 |

### C. Potential truncation effects

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.*^{22} 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.*^{17} and Patel *et al.*^{38} 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.^{25} 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.^{39} Based on extensive two-phase molecular dynamics simulations, Fern *et al.*^{17} and Patel *et al.*^{38} 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^{25} 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.

### D. Total density effects

Valleau^{33} 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^{32} and, more recently, Cortes-Morales *et al.*^{26} 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^{47}). 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.

### E. Finite-size effects

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^{11} 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^{43} 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^{14} and Okumura and Yonezawa.^{7} The latter used simulations in the canonical ensemble and larger finite-size effects are expected for a closed ensemble. Caillol,^{14} Potoff and Panagiotopoulos,^{15} and Pérez-Pellitero *et al.*^{34} 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^{15} and Pérez-Pellitero *et al.*^{34} using FSS. Overall, as already suggested by Smit *et al.*,^{11} 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^{15} and Pérez-Pellitero *et al.*^{34} appear to yield slopes with opposite sign. For *p*_{c}, both the data by Potoff and Panagiotopoulos^{15} 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^{15} and Pérez-Pellitero *et al.*^{34}

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 2*r*_{cut}. This constraint limits the volume integral in Eq. (1) to [(2*r*_{cut})^{3}, *V* − (2*r*_{cut})^{3}]. 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 \u2032 = V 1 / 3 /2$ when volume moves are attempted or accepted that yield configurations with *V*_{I} < (2*r*_{cut})^{3} or *V*_{I} > *V*_{t} − (2*r*_{cut})^{3}, 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,^{11} 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^{44} do not observe a cusp in the restricted Gibbs ensemble where volume moves are not included.

### F. Vapor–liquid equilibria for hard-core square-well particles

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^{45} 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.

λ/σ
. | Reference . | T_{c}/ϵ
. | ρ_{c}/σ^{3}
. | $ T c 400 / T c 800 $ . | $ T c 800 / T c \u221e $ . |
---|---|---|---|---|---|

1.25 | This work | 0.762 ± 0.004 | 0.389 ± 0.026 | 1.001 ± 0.009 | 1.007 ± 0.006 |

Li et al.^{30} | 0.757 0 ± 0.000 4 | 0.3978 ± 0.0022 | |||

1.50 | This work | 1.225 ± 0.007 | 0.308 ± 0.020 | 1.005 ± 0.010 | 1.006 ± 0.006 |

Orkoulas et al.^{27} | 1.218 0 ± 0.000 2 | 0.310 ± 0.001 | |||

Orkoulas et al.^{28} | 1.217 9 ± 0.000 3 | 0.3067 ± 0.0004 | |||

Singh et al.^{29} | 1.217 2 ± 0.000 7 | 0.3079 ± 0.0002 | |||

1.75 | This work | 1.810 ± 0.007 | 0.268 ± 0.015 | 1.003 ± 0.007 | 1.003 ± 0.004 |

Singh et al.^{29} | 1.809 ± 0.002 | 0.2653 ± 0.0018 | |||

Li et al.^{30} | 1.800 37 ± 0.000 08 | 0.2642 ± 0.0001 | |||

2.00 | This work | 2.689 ± 0.015 | 0.265 ± 0.015 | 0.996 ± 0.007 | 1.004 ± 0.006 |

Singh et al.^{29} | 2.68 ± 0.01 | 0.251 ± 0.026 | |||

Li et al.^{30} | 2.676 4 ± 0.001 4 | 0.2573 ± 0.0003 |

λ/σ
. | Reference . | T_{c}/ϵ
. | ρ_{c}/σ^{3}
. | $ T c 400 / T c 800 $ . | $ T c 800 / T c \u221e $ . |
---|---|---|---|---|---|

1.25 | This work | 0.762 ± 0.004 | 0.389 ± 0.026 | 1.001 ± 0.009 | 1.007 ± 0.006 |

Li et al.^{30} | 0.757 0 ± 0.000 4 | 0.3978 ± 0.0022 | |||

1.50 | This work | 1.225 ± 0.007 | 0.308 ± 0.020 | 1.005 ± 0.010 | 1.006 ± 0.006 |

Orkoulas et al.^{27} | 1.218 0 ± 0.000 2 | 0.310 ± 0.001 | |||

Orkoulas et al.^{28} | 1.217 9 ± 0.000 3 | 0.3067 ± 0.0004 | |||

Singh et al.^{29} | 1.217 2 ± 0.000 7 | 0.3079 ± 0.0002 | |||

1.75 | This work | 1.810 ± 0.007 | 0.268 ± 0.015 | 1.003 ± 0.007 | 1.003 ± 0.004 |

Singh et al.^{29} | 1.809 ± 0.002 | 0.2653 ± 0.0018 | |||

Li et al.^{30} | 1.800 37 ± 0.000 08 | 0.2642 ± 0.0001 | |||

2.00 | This work | 2.689 ± 0.015 | 0.265 ± 0.015 | 0.996 ± 0.007 | 1.004 ± 0.006 |

Singh et al.^{29} | 2.68 ± 0.01 | 0.251 ± 0.026 | |||

Li et al.^{30} | 2.676 4 ± 0.001 4 | 0.2573 ± 0.0003 |

### G. Vapor–liquid equilibrium for n-decane

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.94*T*_{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^{46} 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.

. | T_{c}
. | ρ_{c}
. | p_{c}
. |
---|---|---|---|

Method (N_{t}/r_{cut})
. | [K] . | [kg/m^{3}]
. | [kPa] . |

Ambrose and Tsonopoulos | 617.7 ± 0.6 | 228 ± 5 | 2110 ± 50 |

1600/30 Å | 617.8 ± 1.4 | 236 ± 3 | 2300 ± 40 |

1600/20 Å | 618.3 ± 1.7 | 236 ± 3 | 2300 ± 60 |

1600/14 Å | 619.8 ± 0.7 | 237 ± 2 | 2380 ± 30 |

1600/30 Å/SP | 611.3 ± 1.6 | 239 ± 3 | 2280 ± 50 |

200/14 Å | 617.9 ± 1.2 | 237 ± 2 | 2300 ± 50 |

. | T_{c}
. | ρ_{c}
. | p_{c}
. |
---|---|---|---|

Method (N_{t}/r_{cut})
. | [K] . | [kg/m^{3}]
. | [kPa] . |

Ambrose and Tsonopoulos | 617.7 ± 0.6 | 228 ± 5 | 2110 ± 50 |

1600/30 Å | 617.8 ± 1.4 | 236 ± 3 | 2300 ± 40 |

1600/20 Å | 618.3 ± 1.7 | 236 ± 3 | 2300 ± 60 |

1600/14 Å | 619.8 ± 0.7 | 237 ± 2 | 2380 ± 30 |

1600/30 Å/SP | 611.3 ± 1.6 | 239 ± 3 | 2280 ± 50 |

200/14 Å | 617.9 ± 1.2 | 237 ± 2 | 2300 ± 50 |

## V. CONCLUSIONS

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^{15} using mixed-field theory and by Pérez-Pellitero *et al.*^{34} 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.^{43}

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.

## Acknowledgments

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.

### APPENDIX: DETERMINATION OF CRITICAL PROPERTIES AND THEIR UNCERTAINTIES

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,^{21} 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,^{21} the uncertainties in $ \sigma \rho liq + \rho vap 2 $ and $ \sigma \rho liq \u2212 \rho 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}*, $ \sigma a r 2 $,

*b*, $ \sigma b r 2 $, cov(

_{r}*a*,

_{r}*b*),

_{r}*a*, $ \sigma a s 2 $,

_{s}*b*, $ \sigma b s 2 $, and cov(

_{s}*a*,

_{s}*b*) that then give

_{s}Note that the expression for $ \sigma B 2 $ is updated here to eliminate minor typographical mistakes in Ref. 21.