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 Tc = 1.3128 ± 0.0016, ρc = 0.316 ± 0.004, and pc = 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 rcut = 3.5σ yield Tc and pc that are higher by 0.2% and 1.4% than simulations with rcut = 5 and 8σ but still with overlapping 95% confidence intervals. In contrast, GEMC simulations with a truncated and shifted potential show that rcut = 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 constraints3 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 ensemble9–12 or the grand-canonical ensemble13–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) approaches18 to extrapolate the critical properties to the thermodynamic limit. Potoff and Panagiotopoulos15 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 Tc = 1.3120(7), ρc = 0.316(1), and pc = 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. Caillol14 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 Tc = 1.326(2) and ρc = 0.316(2), where the statistical uncertainties are the standard deviations obtained by a block analysis. Okumura and Yonezawa8 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 Tc = 1.3207(4), ρc = 0.316(1), and pc = 0.1288(5), where the statistical uncertainties are the standard deviations reported by the authors.
The Gibbs ensemble Monte Carlo (GEMC) technique pioneered by Panagiotopoulos9,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 law19 and the law of rectilinear diameters20 to determine the critical properties,21 Smit22 obtained Tc = 1.316(6) and ρc = 0.304(6) for the LJ potential with rcut = 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 Tc = 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 (rcut = 3.5, 5, and 8σ with analytical tail corrections and a shifted LJ potential with rcut = 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 < 2rcut), 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 Panagiotopoulos9 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, Vt, is held constant. The Gibbs ensemble partition function considers the distributions of Nt particles and Vt over the two subsystems and all possible configurations as follows (for a one-component system):
where NI and VI denote the number of particles and volume of box I, respectively; rNI and rNt−NI are the positions of the particles in boxes I and II, respectively; UI and UII 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) × 105 and 9.6 × 105 Monte Carlo cycles (MCCs), where each MCC consists of Nt 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) algorithm24 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 rcut = 8σ with analytical tail corrections for energy and pressure,25 and an overall density of ρt = Nt/Vt = 0.316. Since the lengths of boxes I and II (with LI = VI1/3 and LII = (Vt − VI)1/3) fluctuate in a GEMC simulation, the common practice is to use a fixed rcut 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 rcut value tied to a specific force field). Recently, Cortes et al.26 explored different GEMC simulation protocols and recommended that, for temperatures below 0.9Tc, Vt 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 Tc, and different rcut 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 Nt = 5500 and ρt = 0.316. In addition, we perform GEMC simulations using a shifted potential to assess whether rcut = 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 Nt = 5500 and rcut = 5σ. Furthermore, in order to assess the effect of system size, GEMC simulations for Nt = 5500, 1800, and 800 for rcut = 3.5σ and ρt = 0.316 are used and also for an even smaller system containing only 216 particles with rcut = 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 rcut = 5σ. The production periods of these NVT simulations consist of 2.4 × 105 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 × 105 and 2 × 106 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 simulations31 in the Gibbs ensemble are also carried out for n-decane molecules described by the TraPPE–UA force field.32 Two system sizes (Nt = 200 and 1600) and, only for the larger system, four potential truncations schemes (rcut = 14, 20, and 30 Å with analytical tail corrections and a shifted potential with rcut = 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 × 104 MCCs for both system sizes and of 2 × 105 and 8 × 104 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 Nt, Vt, 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 NI/VI and NII/VII data at the end of every MCC (total of 1.92 × 106 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 × 107 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 Nt = 5500, ρt = 0.316, and rcut = 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 × 105 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 Nt = 5500. In contrast, at T = 1.305, the vapor and liquid regions are connected by non-zero probabilities.
Density probability distributions (dashed lines) and Gaussian peak fits using only the data above 0.75 of maximum peak height (solid lines) obtained from GEMC simulations for Nt = 5500, ρt = 0.316, and rcut = 8σ with tail corrections: (top) data shown individually for the eight independent simulations at T = 1.30; (bottom) cumulative data averaged over eight independent simulations shown for all six temperatures.
Density probability distributions (dashed lines) and Gaussian peak fits using only the data above 0.75 of maximum peak height (solid lines) obtained from GEMC simulations for Nt = 5500, ρt = 0.316, and rcut = 8σ with tail corrections: (top) data shown individually for the eight independent simulations at T = 1.30; (bottom) cumulative data averaged over eight independent simulations shown for all six temperatures.
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 Tc, respectively) for Nt = 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.
Density probability distributions obtained from Gibbs ensemble (left) and canonical ensemble (right) simulations at T = 1.30 (top) and 1.36 (bottom) for Nt = 5500 and ρt = 0.316 (vertical dashed lines). The unsorted density distributions for boxes I and II in the GEMC simulations or the left and right halves of the box in the canonical simulations are shown as red and blue lines, respectively. The sorted density distributions for the GEMC boxes and the half boxes in the canonical ensemble with the instantaneously lower and higher densities are shown as magenta and cyan lines, respectively. The diamonds of the corresponding colors show the averages obtained from unsorted and sorted density distributions. The black lines and circles show the Gaussian peak fits to the cumulative DPDs and the corresponding centroids.
Density probability distributions obtained from Gibbs ensemble (left) and canonical ensemble (right) simulations at T = 1.30 (top) and 1.36 (bottom) for Nt = 5500 and ρt = 0.316 (vertical dashed lines). The unsorted density distributions for boxes I and II in the GEMC simulations or the left and right halves of the box in the canonical simulations are shown as red and blue lines, respectively. The sorted density distributions for the GEMC boxes and the half boxes in the canonical ensemble with the instantaneously lower and higher densities are shown as magenta and cyan lines, respectively. The diamonds of the corresponding colors show the averages obtained from unsorted and sorted density distributions. The black lines and circles show the Gaussian peak fits to the cumulative DPDs and the corresponding centroids.
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 Tc, 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 Valleau33 from actual GEMC simulations and from analysis of temperature-and-density-scaling Monte Carlo simulations. Valleau33 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 Tc. However, the differences between the sorted, separated, and combined DPDs vanish for GEMC simulations at temperatures well below Tc 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 Nt = 5500, ρt = 0.316, and rcut = 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 Tc that extrapolation to the critical point using the scaling law and the law of rectilinear diameter19,20 should hold very well. The error propagation procedure is described in the Appendix. For this system, we find Tc = 1.3128 ± 0.0016 and ρc = 0.316 ± 0.004. Once Tc is known, the critical pressure is obtained by an Antoine fit35 to the saturated vapor pressures, and we find pc = 0.1274 ± 0.0013. These results are in excellent agreement with those obtained by Potoff and Panagiotopoulos15 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 Tc, ρc, and pc are only about 0.06%, 0.0%, and 0.4%, respectively. It should be noted that the statistical uncertainties provided by Potoff and Panagiotopoulos15 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 Nt = 5500.
Critical temperature, density, and pressure obtained by different sets of GEMC simulations. SP, VC, and MIM denote use of a shifted potential, a variable rcut, and minimum-image-convention for rcut > L/2, respectively. The 95% confidence intervals are reported for all GEMC data. The FSS data from Potoff and Panagiotopoulos15 using mixed-field theory and Pérez-Pellitero et al.34 using the Binder cumulant parameter are provided for comparison.
Method (Nt/ρtot/rcut) . | Tc . | ρc . | pc . |
---|---|---|---|
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 (Nt/ρtot/rcut) . | Tc . | ρc . | pc . |
---|---|---|---|
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 Tc and pc are about a factor of two smaller than the 95% confidence intervals, but there is obviously a trend of systematically increasing Tc and pc 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 Tc and pc 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: Tc = 1.3124 ± 0.0016, ρc = 0.316 ± 0.004, and pc = 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.
Critical temperature, density, and pressure obtained with different critical exponents. Data and 95% confidence intervals are for GEMC simulations using Nt = 5500, ρt = 0.316, and rcut = 8σ with tail corrections.
β . | Tc . | ρc . | pc . |
---|---|---|---|
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 |
β . | Tc . | ρc . | pc . |
---|---|---|---|
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 rcut = 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 rcut 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 rcut are commonly approximated by analytical tail corrections assuming that the pair distribution functions have converged to unity beyond the cutoff distance.25 As rcut 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 rcut at the cost of increasing computational load.
Fig. 3 illustrates the cumulative DPDs for various approaches at T = 1.30. For the large system (Nt = 5500 and ρt = 0.316), data are also presented for rcut values of 5 and 3.5σ with tail corrections in addition to the rcut = 8σ data introduced in Sec. IV B. The DPDs for the two larger rcut values match extremely well, but the DPD for the liquid phase with rcut = 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 rcut = 8 and 5σ match extremely well over the entire temperature range investigated here, whereas the phase dome is slightly wider for rcut = 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 rcut = 3.5σ yields increases of 0.2%, 0.3%, and 1.4% for Tc, ρc, and pc, respectively, compared to the data for rcut = 8σ. Although their 95% confidence intervals overlap, the differences in Tc and pc 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 rcut value of at least 5σ with tail corrections.
Density probability distributions (dashed lines) and Gaussian peak fits (solid lines) obtained from GEMC simulations for different system sizes, overall densities, and cutoff distances (all with tail corrections), at T = 1.30 averaged over eight independent simulations. Data for Nt = 5500, ρt = 0.316, and different cutoff distances of rcut = 8σ, 5σ, and 3.5σ are shown as orange, black, and magenta lines, respectively. Data with different overall densities of ρt = 0.296 and 0.336 are shown as red and green lines, respectively. Data for ρsys = 0.316, rcut = 3.5σ, and different system sizes of Nt = 5500, 1800, and 800 are shown as magenta, cyan, and blue lines, respectively.
Density probability distributions (dashed lines) and Gaussian peak fits (solid lines) obtained from GEMC simulations for different system sizes, overall densities, and cutoff distances (all with tail corrections), at T = 1.30 averaged over eight independent simulations. Data for Nt = 5500, ρt = 0.316, and different cutoff distances of rcut = 8σ, 5σ, and 3.5σ are shown as orange, black, and magenta lines, respectively. Data with different overall densities of ρt = 0.296 and 0.336 are shown as red and green lines, respectively. Data for ρsys = 0.316, rcut = 3.5σ, and different system sizes of Nt = 5500, 1800, and 800 are shown as magenta, cyan, and blue lines, respectively.
Vapor–liquid coexistence curves obtained from GEMC simulations with different Nt, ρt, and rcut parameters. The symbols show the orthobaric and average densities and the critical points. The lines show the fits to the scaling law and the law of rectilinear diameters. (Top) Data probing cutoff effects with values provided in the legend (SP stands for shifted potential) for Nt = 5500 and ρt = 0.316. (Middle) Data probing total density effects with ρt values provided in the legend for Nt = 5500 and rcut = 5σ. (Bottom) Data probing system size effects with Nt values provided in the legend for ρsys = 0.316 and rcut = 3.5σ for the three larger systems and rcut = 2.5σ for the 216-particle system.
Vapor–liquid coexistence curves obtained from GEMC simulations with different Nt, ρt, and rcut parameters. The symbols show the orthobaric and average densities and the critical points. The lines show the fits to the scaling law and the law of rectilinear diameters. (Top) Data probing cutoff effects with values provided in the legend (SP stands for shifted potential) for Nt = 5500 and ρt = 0.316. (Middle) Data probing total density effects with ρt values provided in the legend for Nt = 5500 and rcut = 5σ. (Bottom) Data probing system size effects with Nt values provided in the legend for ρsys = 0.316 and rcut = 3.5σ for the three larger systems and rcut = 2.5σ for the 216-particle system.
Clausius-Clapeyron plots of the logarithm of the saturated vapor pressure as function of inverse temperature. The symbols show the saturated vapor and critical pressures. The solid lines show the Antoine fits. (Top) Data probing cutoff effects with values provided in the legend (SP stands for shifted potential) for Nt = 5500 and ρt = 0.316. (Middle) Data probing total density effects with ρt values provided in the legend for Nt = 5500 and rcut = 5σ. (Bottom) Data probing system size effects with Nt values provided in the legend for ρsys = 0.316 and rcut = 3.5σ for the three larger systems and rcut = 2.5σ for the 216-particle system. The statistical uncertainties in the psat data are similar to the symbol size and for reasons of clarity not shown.
Clausius-Clapeyron plots of the logarithm of the saturated vapor pressure as function of inverse temperature. The symbols show the saturated vapor and critical pressures. The solid lines show the Antoine fits. (Top) Data probing cutoff effects with values provided in the legend (SP stands for shifted potential) for Nt = 5500 and ρt = 0.316. (Middle) Data probing total density effects with ρt values provided in the legend for Nt = 5500 and rcut = 5σ. (Bottom) Data probing system size effects with Nt values provided in the legend for ρsys = 0.316 and rcut = 3.5σ for the three larger systems and rcut = 2.5σ for the 216-particle system. The statistical uncertainties in the psat data are similar to the symbol size and for reasons of clarity not shown.
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 rcut ≥ 6σ is sufficient to eliminate truncation effects. Here, we also carry out GEMC simulations for the large system (Nt = 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 simulations25 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 rcut). 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 Tc is lowered by about 0.5%, well outside the combined 95% confidence intervals. It turns out that GEMC simulations using Nt = 5500 and a shifted potential with rcut = 8σ yield less accurate results than those using N = 216 and rcut = 2.5σ but with tail corrections, where the former are more expensive by a factor of about 500.
D. Total density effects
Valleau33 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 Nt = 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 Siepmann32 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 Tc 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 Nt = 5500 and rcut = 5σ with tail corrections, i.e., a similar range of ρt values but a larger Nt 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 material47). 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 Tc, ρc, and pc, 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 rcut = 3.5σ with tail corrections are shown in Fig. 3. The simulations for the largest system with Nt = 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 Nt 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 Nt 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 Nt = 5500 exhibits a slightly larger curvature than found for Nt = 800, and the psat values are shifted downward with decreasing system size (they are about 0.5% lower for Nt = 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 Nt = 216, but this requires a reduction of rcut to 2.5σ. The data for the 216-particle system follow the trends of lower-curvature phase dome, downward shift in psat, and decrease in Tc and pc. 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 rcut 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 DPDs11 and Gaussian peak fits and extrapolation to the critical point via the scaling and rectilinear laws yields Tc 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 Panagiotopoulos43 did not observe significant finite-size effects for the three-dimensional LJ fluid for 200 ≤ Nt ≤ 400, but critical properties were not reported.
Finite-size behavior of the critical temperature (top), critical density (middle), and critical pressure (bottom). The green right triangles, dark green squares, purple diamonds, and red down triangles denote FSS studies by Potoff and Panagiotopoulos,15 Caillol,14 Okumura and Yonezawa,8 and Pérez-Pellitero et al.,34 respectively. The black star denotes the large-scale GEMC simulation by Mick et al.23 The filled symbols denote data from the current GEMC study with the simulation parameters provided in the legend.
Finite-size behavior of the critical temperature (top), critical density (middle), and critical pressure (bottom). The green right triangles, dark green squares, purple diamonds, and red down triangles denote FSS studies by Potoff and Panagiotopoulos,15 Caillol,14 Okumura and Yonezawa,8 and Pérez-Pellitero et al.,34 respectively. The black star denotes the large-scale GEMC simulation by Mick et al.23 The filled symbols denote data from the current GEMC study with the simulation parameters provided in the legend.
Comparatively large system-size dependencies are observed for the studies by Caillol14 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 Tc 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 rcut = 8σ with tail corrections are in extremely good agreement with the thermodynamic limit obtained by Potoff and Panagiotopoulos15 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 pc than those for Tc. Thus, with the exception of the canonical ensemble simulations, there are no significant trends. For ρc, the data from Potoff and Panagiotopoulos15 and Pérez-Pellitero et al.34 appear to yield slopes with opposite sign. For pc, both the data by Potoff and Panagiotopoulos15 and our GEMC data indicate a weak increase as the system size increases. Again, our most accurate data (Nt = 5500 and rcut = 8σ with tail corrections) agree very well with the estimates at the thermodynamic limit obtained by Potoff and Panagiotopoulos15 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 Tc, 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.
Number–volume probability distributions for different system sizes at T = 1.305 from GEMC simulations with ρt = 0.316. Data are shown for rcut = 3.5σ with tail corrections and Nt = 5500 (top left), 1800 (top right), and 800 (bottom left), and for rcut = 2.5σ with tail corrections and Nt = 216. In these simulations, volume move attempts that would result in VI < (2rcut)3 or VI > Vt − (2rcut)3 are rejected. Please note the different color scales.
Number–volume probability distributions for different system sizes at T = 1.305 from GEMC simulations with ρt = 0.316. Data are shown for rcut = 3.5σ with tail corrections and Nt = 5500 (top left), 1800 (top right), and 800 (bottom left), and for rcut = 2.5σ with tail corrections and Nt = 216. In these simulations, volume move attempts that would result in VI < (2rcut)3 or VI > Vt − (2rcut)3 are rejected. Please note the different color scales.
Number–volume probability distributions (top) and corresponding density probability distributions (bottom) for cutoff treatments at T = 1.315 from GEMC simulations with Nt = 216, ρt = 0.316, and rcut = 2.5σ with tail corrections. Data are shown for simulations where volume move attempts that would result in VI < (2rcut)3 or VI > Vt − (2rcut)3 are rejected (top left) or the value of rcut for the box with the smaller volume is instantaneously adjusted to be equal to half of its linear dimension (top right).
Number–volume probability distributions (top) and corresponding density probability distributions (bottom) for cutoff treatments at T = 1.315 from GEMC simulations with Nt = 216, ρt = 0.316, and rcut = 2.5σ with tail corrections. Data are shown for simulations where volume move attempts that would result in VI < (2rcut)3 or VI > Vt − (2rcut)3 are rejected (top left) or the value of rcut for the box with the smaller volume is instantaneously adjusted to be equal to half of its linear dimension (top right).
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 rcut 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 2rcut. This constraint limits the volume integral in Eq. (1) to [(2rcut)3, V − (2rcut)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 Nt = 5500 are not affected by the volume constraint, the missing corners emerge for Nt = 1800 and are significant for Nt = 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 when volume moves are attempted or accepted that yield configurations with VI < (2rcut)3 or VI > Vt − (2rcut)3, and restoring the original rcut when the volume split between the boxes becomes more balanced; or (ii) allowing such attempts to proceed without adjusting rcut, 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 Nt = 216, ρt = 0.316, and rcut = 2.5σ with tail corrections are provided in Table I. Tc, ρc, and pc 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 Tc that yield very non-Gaussian peaks. It is clearly surprising that simulations of such a small system size can already predict Tc and ρc with reasonable accuracy. As discussed previously, however, the smaller rcut and Nt cause some noticeable differences in the shape of the phase envelope and in psat. The pathological behavior of these small systems is most evident at temperatures just above Tc. As illustrated in Fig. 8 for simulations with Nt = 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 rcut 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 Tc 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 Tc. 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 rcut, so that the corners in N–V space become irrelevant. This also explains why Mon and Binder44 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 fluid45 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.
Critical temperature and density for the hard-core square-well particles. The 95% confidence intervals are reported for all GEMC data. Data from FSS studies are given for comparison. The last two columns provide the ratio of the critical temperatures for the GEMC simulations with 400 over those with 800 particles and for the 800-particle GEMC simulations over the unweighted mean value from the FSS studies.
λ/σ . | Reference . | Tc/ϵ . | ρc/σ3 . | . | . |
---|---|---|---|---|---|
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 . | Tc/ϵ . | ρc/σ3 . | . | . |
---|---|---|---|---|---|
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 rcut 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 psat 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 (Nt = 1600) and the two larger rcut values agree extremely well. Reducing rcut to 14 Å (the standard value for the TraPPE force field) leads to a 2 K (or 0.3%) increase in Tc and an 80 kPa (or 3.5%) increase in pc; these differences are about equal to the combined 95% confidence intervals. A reduction of the linear system dimension by a factor of two (Nt = 200 with rcut = 14 Å) yields a system with critical properties that are remarkably similar to those for Nt = 1600 and rcut = 30 Å, that is, these reflect a similar fortuitous cancellation of finite-size and truncation errors as observed for the LJ fluid (compare data for Nt = 5500 and rcut = 8σ to those for Nt = 800 and rcut = 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.94Tc), 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 rcut 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 rcut = 30 Å yields a narrower phase dome, a significant increase in psat, and a downward shift of Tc by 7 K (or 1.1%) that falls well outside the combined 95% confidence intervals. The predicted Tc, ρc, and pc values with the standard TraPPE rcut value of 14 Å and averaged over the two system sizes deviate from the recommended experimental critical properties46 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.
Vapor–liquid coexistence curves (top) and Clausius-Clapeyron plots (bottom) obtained from GEMC simulations for the TraPPE–UA n-decane model with different Nt and rcut parameters. The symbols show the orthobaric and average densities, saturated vapor pressures, and critical points. The lines show the fits to the scaling and rectilinear laws. The specific Ncut and rcut values are denoted in the legend. Error bars are not shown because symbol size is larger than 95% confidence interval.
Vapor–liquid coexistence curves (top) and Clausius-Clapeyron plots (bottom) obtained from GEMC simulations for the TraPPE–UA n-decane model with different Nt and rcut parameters. The symbols show the orthobaric and average densities, saturated vapor pressures, and critical points. The lines show the fits to the scaling and rectilinear laws. The specific Ncut and rcut values are denoted in the legend. Error bars are not shown because symbol size is larger than 95% confidence interval.
Critical temperature, density, and pressure for TraPPE–UA n-decane obtained by different sets of GEMC simulations. SP denotes use of a shifted potential. The 95% confidence intervals are reported for all GEMC data. The recommended experimental values46 are provided for comparison.
. | Tc . | ρc . | pc . |
---|---|---|---|
Method (Nt/rcut) . | [K] . | [kg/m3] . | [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 |
. | Tc . | ρc . | pc . |
---|---|---|---|
Method (Nt/rcut) . | [K] . | [kg/m3] . | [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 Tc = 1.3128 ± 0.0016, ρc = 0.316 ± 0.004, and pc = 0.1274 ± 0.0013 that are in excellent agreement with finite-size scaling studies by Potoff and Panagiotopoulos15 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σ ≤ rcut ≤ 8σ) and total densities (0.296 ≤ ρt ≤ 0.336). We also notice that the critical properties obtained for different system sizes (800 ≤ Nt ≤ 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 Tc. From analysis of two-dimensional probability histograms in the N–V plane, we attribute the appearance of a third peak near ρt at T ≈ Tc 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 Tc) and only the four highest temperatures for n-decane (T > 0.87 Tc) 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 = NI/VI and ρvap = (NI − NI)/(Vt − VI), the orthobaric densities are not actually statistically independent. In Smit’s thesis,21 the uncertainties in and are written in terms of NI and VI, 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 ar, , br, , cov(ar, br), as, , bs, , and cov(as, bs) that then give
Note that the expression for is updated here to eliminate minor typographical mistakes in Ref. 21.