Phase equilibria and critical behavior of square‐well fluids of variable width by Gibbs ensemble Monte Carlo simulation

Phase equilibria and critical behavior of square-well fluids of variable width by Gibbs ensemble Monte Carlo simulation Lourdes Vega, Enrique de Miguel, and Luis F. Ruil Departamento de fisica Atbmica Molecular y Nuclear, Universidad de SeviUa, Aptdo. 1065, Sevilla 41080, Spain George Jackson and Ian A. McLure Department of Chemistry, University of Shefieeld, Shefield, S3 7HF, United Kingdom (Received 20 May 199 1; accepted 11 October 199 1) The vapor-liquid phase equilibria of square-well systems with hard-sphere diameters o, welldepths E, and ranges il = 1.25, 1.375, 1.5, 1.75, and 2 are determined by Monte Carlo simulation. The two bulk phases in coexistence are simulated simultaneously using the Gibbs ensemble technique. Vapor-liquid coexistence curves are obtained for a series of reduced temperatures between about T,. = T/T, = 0.8 and 1, where T, is the critical temperature. The radial pair distribution functions g(r) of the two phases are calculated during the simulation, and the results extrapolated to give the appropriate contact values g(a), g(/Za- ), and g(;la -I- ). These are used to calculate the vapor-pressure curves of each system and to test for equality of pressure in the coexisting vapor and liquid phases. The critical points of the squarewell fluids are determined by analyzing the density-temperature coexistence data using the first term of a Wegner expansion. The dependence of the reduced critical temperature Tr = kT,/q pressure P r = P,o-‘/E, number density p: = pE d, and compressibility factor Z = P /(pkT), on the potential range il, is established. These results are compared with existing data obtained from perturbation theories. The shapes of the coexistence curves and the approach to criticality are described in terms of an apparent critical exponent 8. The curves for the square-well systems with il = 1.25, 1.375, 1.5, and 1.75 are very nearly cubic in shape corresponding to near-universal values ofp (flzO.325). This is not the case for the system with a longer potential range; when ;1 = 2, the coexistence curve is closer to quadratic in shape with a nearclassical value of p (PzO.5). These results seem to confirm the view that the departure of fl from a mean-field or classical value for temperatures well below critical is unrelated to longrange, near-critical fluctuations. I. INTRODUCTION Systems of particles interacting via the square-well potential have been extensively studied by statistical mechanical methods. Such interest is invariably due to the fact that the square-well potential is the simplest model which includes the presence of attractive and repulsive forces. The potential energy u(r) for a pair of square-well particles separated by a distance r is given by i + 03, if r<q u(r) = -E, if a<r</2o; (1) 0, if r&la, where u is the hard-sphere diameter of the particle, il is the reduced range of the potential well, and E is its depth. The model is of substantial theoretical importance in studies of systems with a varying potential range since it can represent three limiting cases, namely, the hard-sphere, the shortrange sticky-sphere, and the long-range van der Waals fluids. There has been a series ofstudies in which integral equations were solved numerically for the structural and thermodynamic properties of square-well systems. These have included solutions of the Percus-Yevick (PY),‘-’ hypernetted chain (HNC) ,‘*’ and mean-spherical approximations (MSA) ‘v9 of the Omstein-Zemike (OZ) integral equation, and of the Yvon-Born-Green (YBG) equation’s2’ based on the superposition approximation. In the context of the present study, it is worth discussing the results of the YBG theory in more detail. Numerical calculations for square-well systems with a potential range of il = 1.85 have shown that the YBG equation appears to have a region of long-range solutions characterized by near-universal values of the critical exponents y = 1.24, /3 = 0.33, and S = 4.4 for the compressibility on the critical isochore, the coexistence curve, and the critical isotherm, respectively.‘7-2’ This is in direct contrast with the results of PY theory in which the critical exponents seem to be classical.22 More accurate studies suggest that the YBG equation does not exhibit a true critical point.23-26 It is this confused state of affairs which, in part, motivated the present study. Perturbation theory has often been used to calculate the thermodynamics properties of square-well fluids.27-33 In this paper we compare the critical constants obtained from the second-order theory3’ for systems with a range of values of ;1 with estimates from computer simulation. A perturbation theory for fluids which uses the square-well fluid as the reference has also been developed by de1 Rio et a1.34V35 They have extended the approach to deal with square-well fluids 2296 J. Chem. Phys. 96 (3), 1 February 1992 0021-9606/92/032296-i 0$06.00 0 1992 American Institute of Physics of variable well width.3c39 The theory gives reasonable agreement with existing simulation data for short and long ranges of the potential and most of the fluid range. In addition to theory, a number of Monte Carlo simulation studies of square-well fluids have been undertaken.4043p8*30 The first molecular dynamics simulations of square-well systems were performed by Alder et al.* By and large, only systems with a potential range of II = 1.5 were examined although some results are available for other potential ranges ( 1.125 < ;I < 2) .30,43 Many temperatures and densities in the fluid state have been investigated, but only sketches of the systems’ phase equilibria have been simulated. Direct simulation techniques can be applied to the study of phase equilibria. 45+46 The simplest method involves simulating the two coexisting phases separated by an interface in a single simulation box. Most of the simulation studies focusing on the properties of the interfacial region have been reviewed by Rowlinson and Widom.47 For sufficiently large system sizes, the methods can be used to study bulk phase coexistence. Chapela er aL4’ have used the molecular dynamics method to study the phase equilibria of the squarewell system with /z = 1.5 by confining the fluid in a box between parallel hard walls that are perpendicular to the z direction. Initially the fluid is homogeneous and is given a density close to the critical value, but a fast spinodal decomposition occurs and a slab of liquid condenses at the center of the box with vapor on either side of it. More recently, Benavides et a1.39 have performed a similar study for the squarewell fluid with/z = 3. However, these types of direct simulation studies of phase equilibria possess significant drawbacks because the method is restricted to a film of liquid confined between parallel plates. Unless the size of the system is very large the confinement causes the coexistence properties of the fluid to be different from the bulk coexistence properties of interest. In fact, the method breaks down altogether for temperatures close to the critical point. The effect of confining a fluid between parallel plates is discussed in Sec. IV with reference to the square-well fluid with /z = 1.5. A more recent approach which is used to directly simulate phase equilibria is the so-called Gibbs ensemble Monte Carlo technique. 49*50 In this method the two coexisting homogeneous phases in thermodynamic equilibrium are simulated in separate boxes without the presence of an interface. The simulation involves three distinct Monte Carlo moves: phase space is sampled randomly in each box by moving the molecules following the usual Metropolis Monte Carlo scheme; a random change in volume is made so that an increase in volume of one box corresponds to a decrease in volume in the other, and vice versa; finally, molecules are interchanged between the two boxes. Since the interface is not included in Gibbs ensemble simulations, the coexistence properties are expected to be a good representation of bulk behavior, and in contrast to the former method, the confinement problem does not arise. It has also been suggested51.52 that because the method allows for fluctuations in the numbers of particles and volumes of the subsystems, the Gibbs ensemble method leads to better estimates of the critical constants of the infinite system and a closer proximity to the critical point can be achieved. Gibbs ensemble Monte Carlo simulations have now been performed for a number of potential models, including pure Lennard-Jones fluids and mixtures,49’50*53954 the symmetrical nonadditive hard sphere system,55 the hard-core two-Yukawa tluid,56 the Lennard-Jones fluid with a quadrupole interaction, ” the nonspherical Gay-Berne fluid,58 square-well diatomics,59 and more realistic models for alkanes and water.60v6’ In this work we present the first results of Gibbs ensemble Monte Carlo simulations for the vapor-liquid phase equilibria of square-well fluids with potential ranges of il = 1.25, 1.375, 1.5, 1.75, and 2. The critical constants for the square-well systems are obtained by analyzing the data using the first term of a Wegner expansion. The purpose of the study is twofold: first, to provide accurate vapor-liquid coexistence data for the square-well systems for comparison with mean-field and perturbation theories; second, to determine the dependence of the critical constants on ;1 and to analyze the overall shape of the coexistence curves in terms of an apparent critical exponent 0. In view of the success of the Gibbs method in estimating the critical point, the approach to the critical point is examined with particular attention placed on the shape of two-phase envelope for the various ranges of the potential. Details of the simulations and the data analysis are discussed in Sets. II and III, and the results are presented in Sec. IV. II. COMPUTER SIMULATIONS In Gibbs ensemble Monte Carlo simulations the coexisting vapor and liquid phases are monitored simultaneously as separate subsystems without the presence of an interface. Here, we present a brief summary of the Gibbs ensemble algorithm of Panagiotopoulos49*50 as applied to square-well fluids with variable range. For a clearer and more detailed description of the technique, the reader is referred to the original papers. A thorough statistical mechanical derivation of the algorithm has been presented by Smit ef a1.,51v52 and the method has also been described in a couple of recent reviews.46*62 At the start of the simulation N, = 256 square-well particles representing one of the coexisting phases are placed on a face-centred-cubic (fee) lattice in a cubic box 1 of volume V, with the usual periodic boundary conditions (PBC) and minimum image convention (MIC) .45 Similarly, N, = 256 particles of the other phase are placed in an equivalent box 2 of volume V2. One box represents the vapor phase and the other the liquid phase so that the total number of particles under study is N = N, + N2 = 5 12. During the course of the simulation, the volumes and the numbers of particles of the two subsystems are allowed to vary in such a way that the temperature T, the total volume V= V, + V, and the total number of particles N = N, + N2 remain constant. Since the simulation ensures that the pressures P, and P2, the chemical potentials p, and p2, and the temperature are equal, random fluctuations will force the subsystems into regions of phase space representing the two coexisting phases. A simulation cycle in the Gibbs ensemble Monte Carlo Vega et&: Phase equilibria of the square-well fluid 2297 J. Chem. Phys., Vol. 96, No. 3, 1 February 1992 2298 Vega eta/.: Phase equilibria of the square-well fluid method comprises three distinct types of moves: N, and N2 trial particle displacements within each box, one trial volume change, and a number of trial particle interchanges which depend on the density of the system. In the case of displacement moves, the particles are chosen and displaced randomly within the boxes following the well established canonical NVT Metropolis scheme.63 The acceptance probability for each subsystem is adjusted to -40% by varying the maximum extent of particle displacements. After the N, and N2 trial displacements in each box, a random change in volume is attempted. The volume changes of the two boxes are coupled, with the overall volume V remaining constant in order to satisfy the condition of equality of pressures in the two regions. The acceptance probability for the volume change is essentially given by the prescription of Wooda for simulations in the isothermal-isobaric NPT ensemble. The maximum allowable volume displacement is adjusted to give a-40% acceptance rate for the new trial configurations. The last stage of the Gibbs ensemble method involves particle interchanges between the two subsystems in order to ensure equality of chemical potential in the two phases. This is achieved by creating a particle at a random position in one subsystem and annihilating a randomly chosen particle in the other region. The creation/annihilation moves are reminiscent of Monte Carlo simulations in the grand canonical ,U VT ensemble. The number of interchange moves that are attempted is adjusted so that - l%-5% of the total number of particles are interchanged per cycle; typically 250-3000 interchange attempts are required for the square-well fluid in the density range of interest. model, the pressure cannot be obtained directly from simulation. The virial equation can be used, however, if the radial pair-distribution function g( r) is known or, more specifically, if its contact values g(a), g(/Za- ), and g(/Za+ ) corresponding to the discontinuities of the potential are known. By determining the distances between pairs of particles in each subsystem once every cycle, histograms are constructed which when normalized give g(r) for the two coexisting phases. The contact values are obtained by extrapolating the data for g(r) .& The equation of state can be obtained from the pressure equation’ P -= 1 +p-p+%m --A3[gbw) -g(Ao+)lI, PkT (2) wherep = N/Vis the number density. Hence, if the contact values of g(r) are known for a given number density, the compressibility factor Z and the pressure of the square-well fluid can be calculated from Eq. (2). The pressures of the coexisting vapor and liquid phases are determined in this way for all of the systems studied. The uncertainties in the values of the pressures are difficult to calculate, but can be estimated from the errors in the mean densities of each phase and the errors incurred when g( r) is extrapolated to give the contact values. The simulation of phase equilibria for a given state point requires a total of 2-7 X 10’ such cycles. After equilibrating the system for 1-2X IO5 cycles, a further l-5 X lo5 cycles are performed to accumulate the averages for the properties of the two coexisting phases. The long runs are required for the higher-density configurations. Verlet neighbor 1istP are used for each subsystem to speed up the calculation of pair interactions, the lists being recalculated every cycle. The mean numbers of particles N, and N, , volumes V, and I’,, and energies E, and E, for the two coexisting phases 1 and 2 are obtained as configurational averages over the accumulation stage of the simulation. In this way the coexistence vapor and liquid densities for the square-well system at a fixed temperature are determined. The errors in the average properties of interest are estimated by calculating the standard deviations for blocks of 25 cycles. The initial configurations used in the studies at low temperature are fee lattices. The final equilibrium configurations obtained from the first run are then used as the starting point of simulations at a higher temperature. By repeating the simulations for a series of temperatures and potential ranges, the vapor-liquid coexistence curves are determined between the critical temperature and - 80% of its value for systems withal = 1.25, 1.375, 1.5, 1.75, and 2. In the following section we discuss how the critical constants are obtained by analyzing the phase equilibrium data and how the approach to the critical point can be described in terms of a Wegner expansion with the corresponding critical exponents. The results of the Gibbs-ensemble simulations for the phase-equilibria and critical behavior of the square-well systems are presented in Sec. IV. III. CRITICAL POINT An advantage of the Gibbs ensemble Monte Carlo technique over other computer simulation methods is that a closer proximity to the critical point can be achieved. By a careful analysis of the surface contributions to the free energy density, Smit et al. ” have shown that vapor-liquid coexistence cannot be observed close to but below the critical point. Since Metropolis Monte Carlo simulations for a canonical ensemble overestimate the value of the critical temperature, the Gibbs ensemble leads to a better estimate of its value. They explain this unexpected observation in terms of fluctuations in the volumes and numbers of particles of the two subsystems allowed in Gibbs ensemble simulations. The pressures of the coexisting phases have also been calculated for each one of these systems in order to determine the corresponding vapor-pressure curves and to verify the equality of pressure for given coexistence points. In the case of a discontinuous potential such as the square-well In a study of particular interest, the available data for a number of systems exhibiting vapor-liquid and liquid-liquid phase equilibria have been re-examined.67 The approach to the critical point was described in terms of an effective critical exponent and a Wegner expansion to account for the correction to scaling outside the range of the pure power law behavior. It was suggested that the mean-field or classical critical exponents never accurately describe the shape of a vapor-liquid coexistence curve. We shall question this conclusion for the coexistence data of the square-well systems with ;1 = 2 and 3 in Sec. IV. The study also showed that the departure of the effective critical exponents from their classical values over a wide range of reduced temperatures well J. Chem. Phys., Vol. 96, No. 3, 1 February 1992 below the critical point is unrelated to long-range, near-criti- cal fluctuations. This observation was confirmed for the Gibbs ensemble simulation data of the Lennard-Jones system49 in which long-range fluctuations could not be represented since the maximum system size investigated was N = 500, and it also ties in with the conclusions of Smit et al.” The overall consequence of this, as far as simulations in the Gibbs ensemble are concerned, is that although the correlation lengths of fluctuations are limited by the finite size of the box, the effective critical exponent has a value close to the universal value obtained by the renormalization group (RG) theory for infinite correlation lengths. The effective critical exponent was first defined and used by Verschaffelt6* in 1896 as a sensitive measure of the shape of a coexistence curve. For a vapor of density pU in coexistence with a liquid of density p,, it is defined as Be = JMp, -pv) alnIt ’ where t = 1 - T/T, and T, is the critical temperature. In the limit of small values oft corresponding to temperatures just below the critical point, (p, -pv) =WIS9 (4) where B, is the leading amplitude term. The universal value determined from RG theory of p = 0.325 implies a cubic shape for the coexistence curve, whilst the classical meanfield value of p = 0.5 represents a quadratic coexistence curve. (3) RG theory also allows corrections to scaling to be calculated outside this asymptotic critical region. Wegner69’70 showed that away from the critical point, Eq. (4) can be written as an expansion of the form (pi-p,) =BoItla+B,ItIP+A’+B*ItIs+2A’+..., (5) where A, is one of the so-called gap exponents taking the RG value of A, = 0.5 for the vapor-liquid systems of interest, and Bj are the correction amplitudes or coefficients. Equation (5) can be used with vapor-liquid coexistence data to estimate the system’s critical temperature T, and critical exponent fi. In order to estimate the critical density pc, an equation for the “diameters” (pu + p, )/2 of the coexistence curve must be used,” The anomaly in the diameter of the coexistence curve characterized by tC, is weak and difficult to observe. By combining Eqs. (5) and (6) an equation for the density on each branch can be obtained, p* =pc+C,ItI~+C21tI+C31tI~+A1+~.. *~(BoItI~+B,ItIs+A~+B21tIP+2A~+...). (7) Here, p _ and p + represent the vapor and liquid phase densities, respectively. Originally, it was hoped that most of the terms in the above expression would be used to analyze the data obtained in this simulation study; however, since the results obtained Vega &al.: Phase equilibria of the square-well fluid 2299 were of insufficient precision for a very accurate data analysis, some of the terms in the expression were neglected. More specifically, we excluded the terms C, and C, in the expression for the diameter, and all the higher terms Bi for i > 0 in the Wegner expansion. In our study, temperatures moderately close to the critical point were investigated corresponding to 0 < It I < 0.2, and the gap exponent terms which describe behavior far from the critical point were expected to be small compared with the leading terms. The errors in the Gibbs ensemble coexistence densities were too large to determine whether or not the extended scaling corrections had to be included. Equation (7) then simplifies to P* =pc +C2ltl *PoIv~ (8) in which we are effectively assuming rectilinear diameters and using only the leading amplitude term (4) of the Wegner expansion. This expression allows us to fit the coexistence data and obtain estimates for pc, T,, p, and the amplitude terms C, and B, . The apparent critical exponent obtained in this way is expected to be very similar to the effective critical exponent defined by Eq. (3) and we shall use it to describe the overall shape of the coexistence curves obtained for the square-well fluids with different values of ;1. In order to estimate the critical pressure PC, the vaporpressure curve of the fluid must be determined. By fitting the vapor-pressure data obtained from the simulation to an equation of the Clausius-Clapeyron” form, the value of PC corresponding to the value of T, obtained from Eq. (8) is calculated. The critical constants for the square-well systems were estimated in this way by fitting the Gibbs ensemble coexistence data to Eqs. (8) and (9) using a nonlinear leastsquares procedure. 66 In the following section we report results for the critical points and exponents obtained for the square-well fluids with varying il, and we re-examine the existing data for the Lennard-Jones, hard-core two- Yukawa, and Gay-Beme potential models. IV. RESULTS AND DISCUSSION The phase equilibria are determined using the Gibbs ensemble Monte Carlo technique described in Sec. II for square-well fluids with potential ranges characterized by il = 1.25, 1.375, 1.5, 1.75, and 2. The resulting vapor-liquid coexistence curves are shown in Figs. l-6, the vapor pressure curves are shown in Fig. 7, and the ;1 dependence of the critical temperature, pressure, density, compressibility factor, and exponent are given in Table VI. In the following discussion, it is useful to reduce the temperature and energy with respect to the square-well deptheas T* = kT/eandE * = E/E, thepressureiswritten in terms of E and the hard-sphere diameter cr of the particles, i.e., P * = P~/E, the density is reduced with respect to (T as p* =pd, and th e compressibility factor is given by Z= P/(pkT) = P*/(p*T*).Thesubscriptconanyofthe variables denotes the critical point values. The vapor-liquid coexistence curves for the square-well J. Chem. Phys., Vol. 96, No. 3,1 February 1992 2300 Vega eta/.: Phase equilibria of the square-well fluid T” 0.72 0.70 0.66 0.66 0.64 I 0.0 FIG. 1. The temperature-density vapor-liquid coexistence curve for a square-well fluid with a potential range of/2 = 1.25. The densities ofcoexisting vapor p: and liquid p: phases (squares) and diameters (p: + p:)/Z (triangles) are obtained from Gibbs ensemble Monte Carlo simulations; the error bars represent one standard deviation. The solid curve and line of rettilineardiameters are obtained by fitting Eq. (8) to the simulation data. The estimated critical point (circle) is also shown. fluid with ;1 = 1.25, 1.375, 1.5, and 1.75 are shown in Figs. l-4, respectively, as T - p projections of the PVT surfaces. The data points represent the results of Gibbs ensemble Monte Carlo simulations, and the continuous solid curves are obtained by a least-squares fit of Eq. (8) to the simulation data. The results of the Gibbs ensemble simulations for the densities, energies, and pressures of the coexisting vapor and liquid phases are summarized in Tables I-IV, and the values of the critical constants are given in Table VI. The values of p as listed in Table VI indicate coexistence curves which are closer to cubic than quadratic in shape as expected from the RG theory value ofp = 0.325. The apparent critical exponent p = 0.34 f 0.02 found for il = 1.75 is significantly larger than those of the systems with a shorter potential range, but the overall shape of the coexistence curve is still near-cubic. Also shown in Fig. 3 as solid triangles is the vaporliquid coexistence simulation data obtained by Chapela et a1.48 They used the molecular dynamics method to study the T* 0.1 0.2 0.3 0.4 0.5 0.6 0.7 c P* FIG. 2. The temperature-density vapor-liquid coexistence curve for a FIG. 4. The temperature-density vapor-liquid coexistence curve for a square-well fluid with a potential range of L = 1.375. See the caption of Fig. square-well fluid with a potential range of/2 = 1.75. See the caption of Fig. 1 1 for more details. for more details. 1.25- 1.20. 1.15. T* 1.10- 1.05. l.OO- 0.95 i- 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.6 P* FIG. 3. The temperature-density vapor-liquid coexistence curve for a square-well fluid with a potential range of A= 1.5. See the caption of Fig. 1 for more details. Also shown on this figure as solid triangles are the corresponding data obtained from molecular dynamics simulations (Ref. 48). phase equilibria of the system with 2 = 1.5 by confining the fluid in a box between parallel hard walls. It is clear that the presence of the interface and the confining walls shifts the phase behavior away for the bulk coexistence values. Because of the large critical fluctuations present in such simulations, the method should not be used to predict phase equilibria for temperatures close to the critical point (0.85 < T/T, < 1.0). The Gibbs ensemble Monte Carlo technique allows a closer proximity to the critical point to be achieved and provides an accurate estimate of bulk phase equilibria. The coexistence curve obtained for the system with /z = 2 from simulations in the Gibbs ensemble are shown in Fig. 5 (also see Tables V and VI). In contrast to the results for the shorter potential ranges, the critical exponent fl= 0.53 f 0.11 found for this system indicates a classical critical behavior with a quadratic shape for the coexistence curve. Also shown in the figure as a dashed curve is the theoretical phase equilibria determined by using the accu- 1.85 1.60 1.75 1.70 T* 1.65 -I 0.1 0.2 0.3 0.4 0.5 0.6 J. Chem. Phys., Vol. 96, No. 3,i February 1992 Vega eta/.: Phase equilibria of the square-well fluid 2301 TABLE I. Vapor-liquid coexistence data from Gibbs ensemble Monte Carlo simulations of N = 512 squarewell molecules withapotential rangeofA = 1.25. Thedensitiesp*, energies E *, andpressuresP* ofthecoexisting vapor and liquid phases are labeled u and I, respectively. The errors in p* and E * represent one standard deviation over blocks of 25 cycles, and the error in P* is estimated from the errors in the density and the extrapolated contact values of g( r). T* P: P: ES W Pf p: 0.66 0.062 f 0.001 0.823 f 0.010 - 0.63 f 0.07 - 4.36 f 0.07 0.027 f 0.002 0.025 f 0.008 0.68 0.060 f 0.003 0.792 f0.018 - 0.55 f 0.06 - 4.16 fO.ll 0.033 f 0.002 0.033 f 0.013 0.70 0.085 f 0.007 0.761 f 0.014 - 0.74 f 0.11 - 3.97 f 0.07 0.037 f 0.005 0.045 f 0.016 0.72 0.103 f 0.006 0.708 f 0.014 - 0.83 f 0.12 - 3.69 f 0.07 0.049 f 0.006 0.054 f 0.015 0.73 0.127 f 0.011 0.684 kO.022 - 0.99 f 0.12 - 3.57 f 0.11 0.046 f0.012 0.059 f0.022 0.75 0.161 f0.016 0.597kO.028 - 1.17kO.14 - 3.17kO.12 0.067*0.016 0.079 i-0.025 0.76 0.208 f 0.021 0.516 f 0.039 - 1.46 f 0.14 - 2.84 f. 0.16 0.072 f 0.026 0.098 f 0.034 TABLE II. Vapor-liquid coexistence data from Gibbs ensemble Monte Carlo simulations of N = 512 squarewell molecules with a potential range of A = 1.375. See Table I for details. T* P: P: E: E: p: p: 0.860 0.092 * 0.004 0.722 f 0.009 - 1.13 f 0.11 - 4.67 f 0.07 0.019 f 0.005 0.022 * 0.010 0.880 0.072 & 0.004 0.687 f 0.008 - 0.78 f 0.06 - 4.43 f 0.05 0.027 f 0.004 0.029 f. 0.015 0.900 0.104 0.008 f 0.681 f 0.016 - 1.08 f 0.15 - 4.38 f 0.10 0.042 f 0.010 0.038 f 0.018 0.920 0.109f0.004 0.641&-0.011 -l.lO*O.lO -4.llf0.09 0.048&0.011 0.051+0.016 0.930 0.110 f 0.003 0.620 * 0.013 - 1.11 f 0.10 - 3.99 f 0.09 0.055 f. 0.008 0.059 + 0.019 0.940 0.150~0.011 0.625 f0.021 - 1.40f0.16 -4.00*0.13 0.062&0.016 0.068*0.024 0.945 0.142 & 0.012 0.608 f 0.021 - 1.34 f 0.17 - 3.90 f 0.13 0.064 f 0.016 0.078 f 0.023 0.950 0.151 *0.015 0.593 kO.032 - 1.46kO.17 -3.81 j10.18 0.083 kO.015 0.077f0.031 0.960 0.171 f0.024 0.551 kO.054 - 1.54kO.19 - 3.58kO.28 0.092*0.021 0.085&0.034 0.970 0.202*0.048 0.476kO.137 - 1.76*0.33 - 3.17kO.77 0.096f0.033 0.101 kO.045 TABLE III. Vapor-liquid coexistence data from Gibbs ensemble Monte Carlo simulations of N = 5 12 squarewell molecules with a potential range of/z = 1.5. See Table I for details. T* P: P: E: ET p: p: 1.00 0.038 f 0.001 0.659 f 0.006 - 0.55 f 0.02 - 5.18 f 0.05 0.031 f 0.001 0.031 f 0.009 1.05 0.053 f 0.002 0.632 f 0.007 - 0.72 f 0.04 - 4.96 * 0.06 0.042 f 0.002 0.044 f 0.012 1.08 0.051 f 0.004 0.599 f 0.015 - 0.66 f 0.08 - 4.63 f 0.11 0.043 f 0.003 0.052 f 0.017 1.10 0.058 * 0.003 0.579 f 0.017 - 0.68 f 0.08 - 4.42 f 0.12 0.047 f 0.003 0.056 & 0.021 1.12 0.078 f 0.004 0.567 f 0.018 - 0.95 f 0.08 - 4.45 f 0.13 0.069 f 0.005 0.063 f 0.024 1.15 0.092 f 0.009 0.535 f 0.022 - 1.12 f 0.11 - 4.23 f 0.14 0.074 f 0.009 0.078 f 0.028 1.18 0.133 f. 0.010 0.503 f 0.022 - 1.52 f 0.12 - 4.01 f 0.14 0.090 + 0.014 0.087 f 0.035 1.20 0.147 & 0.017 0.448 f 0.040 - 1.62 f 0.17 - 3.66 f 0.24 0.111 f 0.020 0.096 f 0.041 TABLE IV. Vapor-liquid coexistence data from Gibbs ensemble Monte Carlo simulations of N = 5 12 squarewell molecules with a potential range of/z = 1.75. See Table I for details. T* P: P: E: E: p: p: 1.55 0.048 0.004 f 0.537 f 0.013 - 0.86 + 0.11 - 6.14 + 0.13 0.056 f. 0.004 0.058 f 0.017 1.57 0.053 f 0.005 0.529 f 0.014 - 0.93 f 0.12 - 6.05 + 0.15 0.051 of: 0.007 0.062 f. 0.024 1.60 0.059 f 0.004 0.514 f 0.018 - 1.04 f 0.11 - 5.89 f 0.18 0.060 ~fr 0.007 0.064 f 0.027 1.62 0.062~0.004 0.503 kO.018 - 1.06-j=O.l1 - 5.78kO.18 0.070~0.006 0.075 *0.031 1.65 0.071 * 0.005 0.483 f 0.019 - 1.18 f 0.12 - 5.57 f 0.18 0.071 f 0.010 0.086 f 0.033 1.68 0.104 k 0.005 0.489 * 0.008 - 1.73 f 0.09 - 5.62 f 0.09 0.095 f 0.016 0.105 f 0.028 1.70 0.106~0.007 0.474-&0.010 - 1.73kO.12 -5.46&0.10 0.112~0.018 0.111 kO.029 1.72 0.116 f 0.004 0.457 f 0.012 - 1.86 f 0.09 - 5.30 f 0.11 0.128 f 0.018 0.126 f 0.032 1.73 0.118~0.008 0.445*0.013 -1.87f0.12 -5.17&0.13 0.131*0.020 0.132*0.035 J. Chem. Phys., Vol. 96, No. 3, 1 February 1992 2302 Vega eta/: Phase equilibria of the square-well fluid 2.7 2.6 T* 2.5 2.4 2.3 ( FIG. 5. The temperature-density vapor-liquid coexistence curve for a square-well fluid with a potential range of/2 = 2. See the caption of Fig. 1 for more details. The dashed curve is obtained from VDW mean-field theory (Ref. 39). rate hard-sphere equation of state of Carnahan and Starling’* together with a mean-field attractive term; this is essentially the van der Waals (VDW) equation of state for the square-well fluid of de1 Rio and co-workers.3~39 The VDW approximation underestimates the critical temperature at T: = 2.64 and overestimates the critical density at p: = 0.249. Although the mean-field coexistence curve is in poor agreement with the simulation data close to the critical temperature it provides a reasonable description of the vapor-liquid equilibria for temperatures far enough from the critical point ( T * < 2.4). This is not unexpected for a squarewell fluid with long-range interactions such as ;1 = 2 since in the limit of infinite range the square-well potential is well described by the VDW approximation. The validity of this approximation for the coexistence curve of the system with ;1= 3 has also been noted in a separate study by Benavides et a1.39 The discrepancy between the VDW theory and simulation is more marked for the systems with the smaller values of R. By using higher-order perturbation theories (see for example Refs. 36-39), a better representation of the vaporliquid coexistence curves of short-range square-well fluids can be achieved. In order to verify the adequacy of the principle of corresponding states for the square-well fluids with different values of& the coexistence curves of Figs. l-5 have been replotted in terms of the reduced variables T, = T/T, and pr = p/p= in Fig. 6. It is clear that the principle holds moderately well for systems with small values of il; the coexistence curves for II = 1.25, 1.375, 1.5, and 1.75 are shown as the solid curves with progressively larger widths. It breaks down, however, for the longer range potential with /z = 2 (dotted curve) corresponding to a more pointed curve in the proximity of the critical point. This suggests that the principle of corresponding states can describe the separate behaviors of systems with small and large A but cannot simultaneously describe both. Also shown as the dashed curve in Fig. 6 is the vapor-liquid equilibria obtained from meanfield theory (cf. Fig. 5). The overall shapes of the coexistence curves change from the cubic character of the system with ;1= 1.25 to the quadratic character of the system with A = 2. The critical pressures of the square-well fluids can be obtained from the simulation results by fitting Eq. (9) to the vapor-pressure data. The logarithm of the vapor pressure is plotted as a function of the reciprocal temperature for the systems with;1 = 1.25, 1.375, 1.5, 1.75, and 2 in Fig. 7. Here, the pressures of the vapor (squares) and liquid (solid squares) are shown; the pressures of the coexisting phases calculated from the Gibbs ensemble simulations are close but not exactly equal to each other. The straight lines shown on the figure represent a least-squares fit of Eq. (9) to the data and from these the critical pressures can be estimated. The results obtained for the systems with different values of ;I are summarized in Table VI. For the system with R = 2, close proximity to the critical point could not be simulated and the error in the estimated value of the critical pressure is expected to be quite large. In order to check the value of critical pressure obtained for this system, an additional NVT Monte Carlo simulation was performed for Tr = 2.764 and pr = 0.225 corresponding to the critical point. A critical pressure of P, = 0.210 f 0.006 was obtained by determining the contact values of the radial distribution function; this result is in good agreement with the value of P r = 0.197 & 0.026 obtained from the Clausius-Clapeyron plot (Fig. 7). The values of the critical constants estimated from the simulation data are summarized in Table VI, and a compari- TABLE V. Vapor-liquid coexistence data from Gibbs ensemble Monte Carlo simulations of N = 512 squarewell molecules with a potential range of 1 = 2. See Table I for details. T* P: P: E: E: p: pt 2.35 0.085 f 0.005 0.571 f 0.007 - 1.90 f 0.08 - 9.42 0.11 f 0.110 0.010 f 0.120 0.019 f 2.45 0.106 & 0.005 0.521 f 0.008 - 2.32 f 0.05 - 8.62 0.14 f 0.128 0.016 f 0.135 0.024 f 2.50 0.108 f 0.010 0.489 f 0.008 - 2.31 f 0.10 - 8.12 0.13 f 0.141 0.018 f 0.132 0.025 f 2.52 0.107 f 0.010 0.472 f 0.019 - 2.21 f 0.15 - 7.87 0.18 f 0.144 0.017 * 0.150 0.032 & 2.53 0.107 f 0.008 0.466 f 0.011 - 2.29 f 0.09 - 7.77 0.17 f 0.147 0.015 f 0.149 0.029 f 2.54 0.098 f 0.011 0.447 f 0.014 - 2.00 f 0.21 - 7.49 0.22 f 0.146 0.016 f 0.155 0.038 f 2.56 0.106 f 0.010 0.457 f 0.017 - 2.16 & 0.16 - 7.63 0.27 f 0.150 0.017 f 0.160 0.041 f 2.57 0.110~0.010 0.431 *to.014 -2.31 kO.11 -7.23*0.20 0.153*0.020 0.158~0.040 2.58 0.120 -&0.008 0.428 kO.011 -2.51 *to.09 - 7.19f0.17 0.156f0.025 0.162~0.&0 J. Chem. Phys., Vol. 96, No. 3,1 February 1992 Vega &al.: Phase equilibria of the square-well fluid 2303 TABLE VI. The critical temperature Ty, pressure Pr, density p:, compressibility factor Z,, and exponent p estimated from the Gibbs ensemble Monte Carlo (MC) data for square-well fluids with variable potential ranges R. The fitted values of 4, C, , A, and Bare also given. The errors are estimated from the respective errors in the densities of the coexisting vapour and liquid phases. Also shown in the table are the corresponding values obtained from second-order perturbation theory (PT) (Ref. 30). R T: p: PF ZC B BO G A B 1.250 MC 0.764 f 0.004 0.081 f 0.015 0.370 0.023 f 0.29 0.07 f 0.28 f 0.04 1.35 f 0.13 0.56 f 0.16 4.95 - 5.70 PT 0.913 0.133 0.34 0.43 1.375 MC 0.974 f 0.010 0.105 f 0.023 0.355 0.045 f 0.30 0.11 f 0.25 f 0.07 1.10 f 0.17 0.40 f 0.48 9.99 - 11.93 PT 1.11 0.148 0.34 0.39 I.500 MC 1.219*0.08 0.108 f0.016 0.299*0.023 0.30f0.07 0.30 f 0.02 1.04 f 0.04 0.27 f 0.13 3.67 - 7.19 PT 1.35 0.153 0.31 0.36 1.750 MC 1.811 f 0.013 0.179 f 0.020 0.284 0.009 f 0.35 0.05 f 0.34 f 0.02 0.95 f 0.01 0.04 f 0.03 5.79 - 13.60 PT 2.04 0.196 0.25 0.38 2.ooO MC 2.764 f 0.023 0.197 f 0.026 0.225 f 0.018 0.32 0.07 j, 0.53 f 0.11 1.31 f0.12 0.72 f 0.08 1.50 - 8.63 PT 2.88 0.255 0.24 0.37 son is made with the corresponding results obtained from second-order thermodynamic perturbation theory.30 The perturbation theory provides only a reasonable description of the simulation values with improving agreement as /2 is increased. Carley 33 has used radial distribution functions calculated from integral equation theory together with firstorder perturbation theory to solve for phase coexistence in the system with R = 1.5, and has obtained critical point values of Tr = 1.35 and &’ = 0.30 f 0.023. The value for the critical density compares favorably with the Gibbs ensemble estimate of pr = 0.299 f 0.023, but the result obtained for the critical temperature is considerably higher than the simulation value of TT = 1.219 f 0.008. Perhaps the most interesting result is the change in the shape of the coexistence curve as the range is increased. The systems withal = 1.25,1.375,1.5, and 1.75 haveanear-cubic shape with apparent critical exponents close to the universal value of/3 = 0.325. A small increase in the critical exponent is, however, apparent. These results are in good agreement with the value of p = 0.33 obtained for the system with R = 1.85 from numerical studies using the YBG integral equation 17-*’ although a more accurate analysis has suggested that the equation does not exhibit a true critical point.23-26 For the system with R = 2 the shape of the coexistence curve is nearly quadratic corresponding to the classical mean-field value of p = 0.5. In order to ensure that this result is not due the small system size of N = 5 12, a larger system of N= 1000 was also examined; the resulting data were very similar to that of the smaller system. In the case of the square-well fluid with il = 3 studied by Benavides et a1.,39 the appare n t c ritical exponent estimated from the simulation data has increased top = 0.77 f 0.15. If their data is analyzed using Eq. (8) the best-fitted parameters are T~=11.68~0.15,p~=0.181f0.010,~=0.77f0.15, Be = 1.32 + 0.06, and C2 = 0.43 f 0.03. It must be noted, however, that since the temperatures investigated for this system are far from the critical value (0.5 < T/T, < 0.8)) it is difficult to estimate the critical constants and exponents. Existing Gibbs ensemble simulation data for the phase equilibria of the Lennard-Jones,49r’o hard-core twol.oM) 0 075 0.950 0.925 0.900 0.075 0.650 ~ - _^ ” “” ".ZS ".SO 0,;s l.bO l.iS 1.50 1,;s 2.bo 2.i5 2.;0 PI FIG. 6. The vapor-liquid coexistence curves for square-well fluids with variable potential range A plotted in terms of the reduced temperature T, = T/T, and density p, =p/pc. The solid curves with progressively larger widths represent the systems with ,l = 1.25, 1.375, 1.5, and 1.75; the dotted curve represents the coexistence curve for the system with R = 2. Also shown as a dashed curve is the coexistence curve obtained from VDW mean-field theory (Ref. 39). -1.25 2 1.75 1.5 1.375 1.25 -1.75 1.25 -2.25 IIlP * -2.75 -3.25 -3.75 ..,,I_\.\; 0.25 0.50 0.75 1 .oo 1.25 1.50 -4.25 0.25 0.50 0.75 1 .oo 1.25 1.50 5 i/r- FIG. 7. The vapor-pressure curves for square-well fluids with variable potential range /2. The curves are labeled with the appropriate values of,%, and the pressures of the coexisting vapor (solid squares) and liquid phases (open squares) are shown together with the estimated critical points (circles) . J. Chem. Phys., Vol. 96, No. 3,1 February 1992 2304 Vega eta/: Phase equilibria of the square-well fluid Yukawa,56 and the Gay-Berne5* fluids have been reanalyzed using Eq. (8) as described in Sec. III. Best fits to the simulation data were obtained with the following values of the parameters: T: = 1.321 f 0.005, p: = 0.321 f 0.017, fl= 0.36 f 0.05, B0 = 1.11 f 0.05, and C, = 0.20 f 0.07 for the Lennard-Jones fluid; T$ = 1.294 f 0.009, /.I: = 0.342 f 0.021, /? = 0.33 f 0.08, B, = 1.07 f 0.07, and C, = 0.12 & 0.04 for the hard-core two-Yukawa fluid; and T,* = 0.488 f 0.004, p: = 0.101 f 0.009, fl= 0.32 f 0.03, &, = 0.43 f 0.03, and C, = 0.16 f 0.05 for the Gay-Berne fluid. The estimated values of the apparent critical exponent indicate that the shapes of the coexistence curves were nearly cubic for these three systems. A considerable amount of care must be taken in determining values of the critical points and critical exponents; a small change in the critical exponent can make quite a significant difference on the estimated values of the critical constants. For example, in obtaining the critical point of water from Gibbs ensemble simulation data, de Pablo et a1.6’ state that a classical critical exponent with a fixed value of /3 = 0.5 was used. Their results would have been quite different had they used the correct universal value of p = 0.325. *D. Henderson, W. G. Madden, and D. D. Fitts, J. Chem. Phys. 64,5026 (1976). 9W. R. Smith, D. Henderson, and Y. Tago, J. Chem. Phys. 67, 5308 (1977). “I. B. Schrodt and K. D. Luks, J. Chem. Phys. 57,200 (1972). ” J. J. Kozak. I. B. Schrodt. and K. D. Luks. J. Chem. Phvs. 57,207 ( 1972). i21, B. Schroht, J. S. Ku, and K. D. Luks, J. Chem. Phyd. 57,4589 ( 1972). l3 I. B. Schrodt, J. J. Kozak, and K. D. Luks, J. Chem. Phys. 60,170 ( 1974). i’ W. Lincoln, Y. Tago, and K. D. Luks, J. Chem. Phys. 61,4129 (1974). ‘$W. Lincoln, J. J. Kozak, and K. D. Luks, J. Chem. Phys. 62,1116 (1975). i6W. Lincoln, J. J.Kozak,andK. D. Luks, J.Chem.Phys.62,2171(1975). “K. U. Co, J. J. Kozak, and K. D. Luks, J. Chem. Phys. 64,2197 ( 1976). 18K. U. Co, K. D. Luks, and J. J. Kozak, Mol. Phys. 36,1883 (1978). 19K. A. Green, K. D. Luks, and J. J. Kozak, Phys. Rev. Lett. 42, 985 (1979). “‘K. A. Green, K. D. Luks, E. Lee, and J. J. Kozak, Phys. Rev. A 21,356 (1980). 2’ K. A. Green, K. D. Luks, and J. J. Kozak, Phys. Rev. A24,2093 ( 1981). 22 S. Fishman and M. E. Fisher, Physica 108A, 1 ( 198 1). 23 G. L, Jones, J. J. Kozak, E. Lee, S. Fishman, and M. E. Fisher, Phys. Rev. Lett.46,795 (1981). “S. Fishman and M. E. Fisher, Phys. Rev. Lett. 47,421 (1981). 25G. L. Jones, E. Lee, and J. J. Kozak, Phys. Rev. Lett. 48,447 ( 1982). r6K. A. Green, K. D. Luks, G. L. Jones, E. Lee, and J. J. Kozak, Phys. Rev. A 25, 1060 (1982). “W. R. Smith, D. Henderson, and J. A. Barker, J. Chem. Phys. 53, 508 (1970). Before we conclude, it is important to point out that the apparent critical exponents estimated here only give the general qualitative shape of the two-phase envelope, and that the shape of some of the coexistence curves may change quite dramatically close to the critical point. Since for square-well systems with larger values of il it is increasingly difficult to simulate state points close to the critical point, the errors in the estimated values of the critical exponents may be considerable. This having been said, the behavior of the square-well fluid is in stark contrast with that observed for the Lcnnard- Jones system.67 The effective critical exponent of the LennardJones fluid is fairly constant and close to that obtained by the RG method (/3 = 0.325) over a wide range of temperatures well below the critical point. This does not seem to be the case for the square-well fluids with long-range interactions in which the critical exponents would have to decrease quite dramatically to the universal value as the critical point is approached, assuming, of course, that these systems exhibit universal critical behavior. r* W. R. Smith, D. Henderson, and J. A. Barker, J. Chem. Phys. 55,4027 (1971). “D. Henderson, J. A. Barker, and W. R. Smith, J. Chem. Phys. 64.4244 (1976). “D. Henderson, 0. H. Scalise, and W. R. Smith, J. Chem. Phys. 72,243l (1980). 3’ D. D. Carley, J. Chem. Phys. 67, 1267 (1977). 32D. D. Carley and A. C. Dotson, Phys. Rev. A 23, 1411 (1981). “D. D. Carley, J. Chem. Phys. 78,5776 (1983). “D. A. de Lonngi and F. de1 Rio, Mol. Phys. 48,293 (1983). 35 D. A. de Lonngi and F. de1 Rio, Mol. Phys. 56,691 (1985). 36F. de1 Rio and L. Lira, Mol. Phys. 61,275 ( 1987). “F. de1 Rio and L. Lira, J. Chem. Phys. 87,7179 (1987). ‘*A. L. Benavides and F. de1 Rio, Mol. Phys. 68,983 (1989). 39A. L. Benavides, J. Alejandre, and F. de1 Rio, Mol. Phys. 74,32 1 ( 199 1). “‘A. Rote&erg, J. Chem. Phys. 43, 1198 ( 1965). 4’ F. Lado and W. W. Wood, J. Chem. Phys. 49,4244 (1968). 42 Y. Rosenfeld and R. Thieberger, J. Chem. Phys. 63, 1875 ( 1975 ). 43K. D. Scarfe, I. L. McLaughlin, and A. F. Collings, J. Chem. Phys. 65, 2991 (1976). *B. J. Alder, D. A. Young, and M. A. Mark, J. Chem. Phys. 56, 3013 (1972). ACKNOWLEDGMENTS We wish to thank Fernando de1 Rio for useful discussions and Ian L. Pegg for suggesting the form of Eq. (7). This research was supported in part by grants from the Acciones Integradas progzamme of the British Council and from the Consejeria de Educaci6n y Ciencia of the Junta de Andalucia. Generous allocations of computer time at the Centro Informltico Cientifico de Andalucia (CICA) are also acknowledged. ‘D. Levesque, Physica 32, 1985 (1966). *L. Verlet and D. Levesque, Physica 36,254 ( 1967). 3 Y. Tago, J. Chem. Phys. 58,2096 ( 1973). ‘Y. Tago, Phys. Lett. 44A, 43 (1973). ‘Y. Tago, J. Chem. Phys. 60, 1528 (1974). 6 W. R. Smith, D. Henderson, and R. D. Murphy, J. Chem. Phys. 61,29 11 (1974). ‘W. R. Smith and D. Henderson, J. Chem. Phys. 69,319 (1978). ” M. P. Allen and D. J. Tildesley, Computer Simulations of Liquids (Clarendon, Oxford, 1987). 46K. E. Gubbins, Mol. Simul. 2,223 ( 1989). “J S Rowlinson and B. Widom, Mokcular Theory ofcapillarity (Claren- . . don, Oxford, 1982). “G A. Chapela, S. E. Martmez-Casas, and C. Varca, J. Chem. Phys. 86, 5683 (1987). “A. Z. Panagiotopoulos, Mol. Phys. 61, 813 (1987). %A. Z. Panagiotopoulos, N. Quirke, M. Stapleton, and D. J. Tildesley, Mol. Phys. 63,527 (1988). r’ B. Smit, Ph. de Smedt, and D. Frenkel, Mol. Phys. 68,931 ( 1989). ‘*B. Smit and D. Frenkel, Mol. Phys. 68,95 1 ( 1989). s3 A. Z. Panagiotopoulos, Mol. Phys. 62,701 ( 1987). “R R Singh, K. S. Pitzer, J. J. de Pablo, and J. M. Prausnitz, J. Chem. . . Phys. 92,5463 ( 1990). “5. G. Amar, Mol. Phys. 67,739 ( 1989). s6E. N. Rudisill and P. T. Cummings, Mol. Phys. 68, 629 (1989). “M. Stapleton and A. Z. Panagiotopoulos, J. Chem. Phys. 92, 1285 (1990). ‘*E. de Miguel, L. F. Rull, M. K. Chalam, and K. E. Gubbins, Mol. Phys. 71, 1223 (1990). 59A. Yethiraj and C. K. Hall, Mol. Phys. 72,619 (1991). 6oJ. J. de Pablo and J. M. Prausnitz, Fluid Phase Equil. 53, 177 ( 1989). 6’J. J. de Pablo, J. M. Prausnitz, H. J. Strauch, and P. T. Cummings, J. Chem. Phys. 93,7355 (1990). J. Chem. Phys., Vol. 96, No. 3, 1 February 1992