An improved approach for predicting the critical constants of large molecules with Gibbs Ensemble Monte Carlo simulation

Fluid Phase Equilibria 425 (2016) 432e442

Contents lists available at ScienceDirect

Fluid Phase Equilibria j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / fl u i d

An improved approach for predicting the critical constants of large molecules with Gibbs Ensemble Monte Carlo simulation Richard A. Messerly, Thomas A. Knotts IV, Richard L. Rowley, W. Vincent Wilding* Department of Chemical Engineering, Brigham Young University, Provo, UT 84602, United States

a r t i c l e i n f o

a b s t r a c t

Article history: Received 6 May 2016 Received in revised form 20 June 2016 Accepted 23 June 2016 Available online 25 June 2016

In this work we focus on predicting the critical temperature (Tc), critical density (rc), and critical pressure (Pc) from Gibbs Ensemble Monte Carlo (GEMC) simulations. Our primary objective is to reduce the uncertainty associated with the critical point constants, particularly Pc, for large molecules. To achieve this goal, we demonstrate the advantages of using the Rackett equation to predict Pc compared to the traditional approach of using the Antoine equation. The main difference is that the Rackett equation utilizes liquid density (rL) while the Antoine equation uses vapor pressure (Pv). The Rackett equation yields a better prediction of Pc than the Antoine equation because rL values are more reliable than Pv values when obtained from GEMC simulations for the standard force field models. As either method will yield large uncertainties in Pc if the uncertainties in rc and/or Tc are large, we also develop a statisticallyrigorous experimental design to minimize the uncertainty in Tc, rc, and Pc. The greatest improvement in uncertainty is found for rc and Pc when compared to other contemporary methods. © 2016 Elsevier B.V. All rights reserved.

Keywords: Molecular simulation Rackett equation Experimental design Uncertainty

The critical temperature (Tc), critical density (rc), and critical pressure (Pc) are important thermophysical properties in science and engineering. The critical constants are required to determine equation of state parameters for PVT calculations and to predict other properties based upon the corresponding states principle [1] and other prediction methods. Experimental measurement of critical constants is possible and reliable data exist for a large number of compounds; however, for a number of reasons no experimental data exist for many compounds of interest. For example, obtaining reliable critical constant values becomes difficult for larger compounds because they thermally decompose at temperatures well below Tc. In other cases, the chemical toxicity or reactivity is such that experiments are inadvisable. Situations such as these, where no experimental critical constants are available, render the prediction methods for other properties ineffective for the very molecules for which they are needed. Due to the importance of the critical properties in engineering practice, multiple groups have made efforts to develop techniques which can obtain the values for molecules that do decompose. For example, the Nikitin group has developed a pulse-heating technique which has permitted experimental measurements of Tc and

* Corresponding author. E-mail address: [email protected] (W.V. Wilding). http://dx.doi.org/10.1016/j.fluid.2016.06.041 0378-3812/© 2016 Elsevier B.V. All rights reserved.

Pc for n-alkanes as large as C60 [2]. However, the Pc values reported by Nikitin et al. for C36 differ significantly between their original and more recent publications [3]. This discrepancy is troubling since both the original and more recent values were obtained from the same experimental data and differ only in the data analysis method. In addition, experimental Pc data measured by Teja for nalkanes up to C18 show a strong discontinuity compared with Nikitin’s data for larger compounds [4]. Molecular simulation has the potential to help resolve the disagreement but, unfortunately, the existing simulation data found in the literature are not reliable enough to discern between the different experimental data sets. Historically, predicted critical constant values from Grand Canonical Monte Carlo (GCMC) finite-size scaling methods were considered superior to those obtained from Gibbs Ensemble Monte Carlo (GEMC). There are at least two key limitations to predicting the critical constants with GEMC. First, the GEMC approach suffers from large fluctuations near the critical point. For this reason, it is necessary to perform GEMC simulations at temperatures below Tc and then extrapolate to the critical point. In a previous publication we demonstrated how to rigorously quantify the uncertainty due to extrapolation [5]. In this work, we demonstrate how to reduce this uncertainty. The second supposed limitation is that there is no rigorous approach to correct for finite-size effects with GEMC [6]. However, it was recently shown that finite-size effects for Tc, rc, and Pc are smaller than the corresponding statistical uncertainty for

R.A. Messerly et al. / Fluid Phase Equilibria 425 (2016) 432e442

GEMC systems with 200 molecules of n-decane [7]. These conclusions suggest that GEMC is an attractive option to obtain reliable critical constants because the computational demand is much less than that for GCMC finite-size scaling methods. Unfortunately, the traditional data analysis approach for predicting Pc from GEMC simulation data yields inaccurate results with large statistical uncertainties. The purpose for the work at hand is to demonstrate an improved approach for estimating the critical constants, particularly Pc, from GEMC simulations. Our main objective is to obtain accurate results with correspondingly small uncertainties. The proposed method is equally rigorous as the approach traditionally found in the literature yet it provides some significant advantages, especially for large compounds, as will be shown. The proposed methodology is composed of two key aspects. First, as Vetere proposed, we utilize the Rackett equation with liquid density (rL) to more accurately predict Pc [8]. Second, we implement an experimental design to predict the ideal temperatures for the GEMC simulations that will minimize the uncertainty in Tc, rc, and Pc. We show that only two temperatures are needed in the optimal experimental design while most GEMC studies perform simulations at several temperatures. The optimal reduced temperatures (TR) for performing the GEMC simulations are determined in two different ways. First, a so-called D-optimal experimental design is developed to minimize the uncertainty in Tc, rc, and Pc. A D-optimal design has the advantage that the uncertainties for all three critical constants are reduced simultaneously. Second, a parameter-specific optimization is performed by finding the two temperatures that minimize the variance in just Tc or rc. This approach is ideal in the case where only one critical constant is of interest. In addition to reducing the statistical uncertainty, we show that this experimental design is less susceptible to finite-size effects which have historically posed a serious limitation for predicting the critical point of large molecules (greater than about C26) [9]. Since this experimental design can reduce uncertainty and finite-size effects without increasing computational time, its use makes GEMC simulation a viable option for studying the thermodynamic behavior of large compounds and should become a valuable tool in this regard. The outline for this document is the following. In Section 1 we explain the benefits of utilizing Vetere’s method for predicting Pc from the Rackett equation and rL simulation data. In Section 2 we explain how the uncertainty in Pc can be rigorously estimated from this method. Then, in Section 3 we demonstrate the need for a Doptimal experimental design to minimize the uncertainty in Pc. In Section 4 the experimental design for GEMC simulations is derived in detail. Next, in Section 5 we present a quantitative comparison of the uncertainties in Tc, rc, and Pc between the experimental design proposed in this work and those found in the literature. This is followed in Section 6 by a brief consideration of modifications to the experimental design, the implications of this experimental design on finite-size effects, and possible limitations of this approach. A step-by-step outline of the proposed methodology is then presented in Section 7. Finally, in Section 8 we discuss our conclusions and the impact they have on predicting the critical constants for large compounds. 1. Alternative method to Predict Pc A primary reason for the present work is to develop a more effective method for estimating Pc from GEMC simulations of larger compounds. The traditional approach to predicting Pc from GEMC simulation results makes use of the Antoine equation that relates vapor pressure (Pv) to temperature (T) [7]. The form typically used

433

in the simulation literature is

lnðPv Þ ¼ A0 þ

A1 T

(1)

where A0 and A1 are fitting parameters. The critical pressure is obtained by first regressing the Pv GEMC data to Equation (1). Then, Tc is obtained by regressing the GEMC results to the law of rectilinear diameters [10] and density scaling law [11]. The law of rectilinear diameters can be expressed as

rr ≡

rL þ rv 2

¼ rc þ AðTc  TÞ

(2)

where rr is the rectilinear density, rv is the vapor density, and A is a fitting parameter. The density scaling law is

rs ≡rL  rv ¼ BðTc  TÞb

(3)

where rs is the scaling density, B is a fitting parameter, and b is a constant, typically 0.326 [7]. Finally, Pc is calculated using Tc and the regression to Equation (1). The traditional approach for predicting Pc is not ideal for larger compounds for several reasons. The first is that it is very difficult to achieve accurate vapor pressure results from molecular simulation. This is partially due to the fact that Pv is an inherently noisy property when obtained from GEMC simulations. For better statistics, more vapor phase molecules are required which necessitates an increase in the overall number of molecules. However, this comes at a great computational cost for longer chains. Furthermore, even without increasing the number of molecules in the vapor phase, obtaining Pv becomes computationally intensive for larger molecules because of the “virial” force calculation. The virial forces require an individual calculation for each unique pair of sites in a simulation [12]. Thus, the calculation cost for the virial forces scales as N2, where N is the number of sites, not the number of molecules. Therefore, assuming the same number of molecules, the computational cost for the virial forces is about 242 times more for C48 H98 than ethane when using a united-atom model. By contrast with Molecular Dynamics, the virial forces are not necessarily computed in Monte Carlo simulations [6]. Therefore, the computationally expensive virial calculation can be completely eliminated by not calculating Pv. Another reason accurate Pv values are difficult to obtain is because most intermolecular potential model functions, for example the Lennard-Jones 12-6 model (LJ 12-6), do not have the right form to accurately model both rL and Pv [13]. For this reason, most potential models (e.g. OPLS, TraPPE, NERD, etc. [14e16]) are optimized to reproduce rL data. Because the extant models were not designed to yield accurate Pv values the traditional approach for predicting Pc from Pv is inherently limited. On the other hand, the Exponential-6 model developed by Errington et al. [17] and the Mie 16-6 model developed by Potoff et al. [18] are capable of reliably predicting Pv within 3%. Despite this fact, both the Exp-6 and Mie 16-6 models, widely considered two of the best for n-alkanes, show systematic deviations in Pc, especially for larger compounds. This poor performance for Pc is rationalized by Potoff et al. as a “result of compound errors in the predictions of critical temperatures and vapor pressures.” In other words, because of the exponential relationship between Pv and T in Equation (1), any deviation or uncertainty in Pv and Tc will be magnified for Pc. The final limitation of the traditional method is that Equation (1) is often not reliable over a large temperature range. The traditional approach seeks to alleviate this deficiency by obtaining several Pv values at temperatures near Tc. However, since fluctuations in simulation results increase near the critical point this only

434

R.A. Messerly et al. / Fluid Phase Equilibria 425 (2016) 432e442

exacerbates the uncertainty in Pv. Therefore, although it should be expected that more accurate Pv values from simulation will be achievable in the future, the shortcomings of the Antoine equation extrapolation will remain relevant. Although more flexible models, such as the Riedel equation, utilize additional terms they are not traditionally used considering the noise in the Pv simulation results [19]. Since most potential models do well modeling rL, Vetere’s approach to estimating Pc is a logical alternative. This method rearranges the Rackett equation [20] to solve for Pc according to




Pc ¼

Rg Tc rc rc Mw rL




g

B 2

rL ¼ rc þ AðTc  TÞ þ ðTc  TÞb

(5)

where A, B, and b are the same parameters found in Equations (2) and (3). This expression is beneficial because the uncertainty in rL can be obtained from the uncertainties in the parameters A, B, rc, and Tc, which are all fit at the same time and thus account for the correlation. Specifically, we showed that a rigorous statistical method for assessing the joint uncertainty in these parameters may be done using the fundamental equation

1þ

p Fp;np;a np

(7)

where PDF(q) is the probability density function for a specific set of parameters and F1 is the inverse of the F-statistic with p and np degrees of freedom. When used with a Monte Carlo sampling propagation of errors approach, the p-dimensional parameter space can be sampled by using the PDF obtained from Equation (7). Specifically, we obtain a four-dimensional PDF for the parameters A,B,Tc, and rc found in Equation (5). Then, we randomly sample from this parameter space using the PDF to weight the probability of each parameter set, q. For each parameter set we calculate rL using Equation (5). We can then determine Pc from Equation (4) for the corresponding set of rL, Tc, and rc. By sampling millions of different parameter sets we generate a histogram of Pc values that can be integrated at the desired confidence level. In short, Equations (5) and (7) allow a rigorous estimate of the uncertainty in a predicted value of Pc from Equation (4). The reader is referred to the previous work for more details on this method [5]. 3. Case study: long-chain-length behavior of n-alkanes

The statistical uncertainty of Pc predicted from Equation (4) depends on the uncertainty in rL, Tc, and rc. Here we will discuss how to propagate the uncertainty in rL, Tc, and rc to Pc. Predicting the uncertainty in rL is commonly done by performing replicate simulations at the same temperature. For example, in a previous study we developed an error model for predicting the standard deviation in rL as a function of TR for a 200 molecule system by performing 20 replicate simulations at each temperature [5]. However, the replicate simulation technique cannot be used to independently calculate the uncertainties of rL, Tc, and rc because they are correlated. In other words, any deviation in rL leads to a deviation in Tc and rc because Tc and rc are regressed to rr and rs that depend upon rL (see Equations (2) and (3)). To account for this interdependence the law of rectilinear diameters and the density scaling law may be combined to obtain




ðÞ

C B C BSðqÞ  S b q np ; p; n  pC PDFðqÞ ¼ F 1 B  C B p A @ S b q

(4)

2. Estimating uncertainty in Pc

ðÞ

1

0

ðÞ

T

1Tci

where Rg is the universal gas constant, Mw is the molecular weight, Ti is the temperature corresponding to rL, and g is an empirical parameter estimated to be around 27. Although Vetere’s method has existed for nearly two decades, to our knowledge it had never been utilized with data from molecular simulations [8]. The Rackett approach to predict Pc is not the most useful of methodologies when Tc and rc are obtained from experiment. The reason is that it is much easier to accurately measure Tc and Pc than rc. For example, Nikitin et al. did not report rc values for any of the large n-alkanes they studied. The situation is quite different for GEMC simulations where obtaining Tc and rc from the density scaling law and the law of rectilinear diameters is straightforward. In addition, rL has a much lower inherent uncertainty than Pv when obtained from GEMC simulations. Therefore, although the Rackett method gained little attention for experimental critical measurements it is an ideal approach for molecular simulation when the selected force field can reliably estimate rL, Tc, and rc.

SðqÞ  S b q

where S(q) is the weighted sum squared error of a specific set of parameters, Sðb q Þ is the minimized weighted sum squared error of the optimal set of parameters, p is the number of parameters, n is the number of data points, and Fp,np,a is the F-statistic at the a confidence level with p and np degrees of freedom. An alternative form of this equation is useful when performing a Monte Carlo sampling propagation of errors analysis for Pc and is

 (6)

We have chosen the n-alkane family for a case study of the Vetere approach. There are two objectives for this section. First, to illustrate the effectiveness of the Vetere method for predicting Pc from GEMC simulations. Second, to demonstrate the need for an optimal experimental design to reduce the uncertainty in Pc. For the former, the analysis will include using existing simulation data but applying the Vetere method to predict Pc rather than the Antoine equation. This comparison provides a rigorous test of the abilities of the method by eliminating biases that could be introduced by using simulation codes or methods different than those of the original authors. In this case study we compare the Pc results obtained from four different force fields (TraPPE [15], NERD [16], Exp-6 [17], and Mie 16-6 [18]). Nath et al. used GEMC simulations to predict Tc and rc for n-alkanes up to C48 using the NERD force field [16]; however, they did not report Pc or even Pv values. To analyze this force field, we used their orthobaric density results to calculate Pc using Vetere’s method. We refer to these Pc values as Nath-Vetere. Likewise, we predicted Pc with Equation (4) for the TraPPE model by using the GEMC data reported by two different articles from Siepmann’s group, Martin et al. for C2-C12 and Zhuravlev et al. for C30 [15,21]. For simplicity, we refer to these values collectively as SiepmannVetere. The Siepmann-Vetere results provide a useful benchmark when assessing the performance of Vetere’s method because Martin et al. and Zhuravlev et al. also reported Pc values that were obtained from the traditional method using Equation (1). Fig. 1 depicts the dependence of Pc upon carbon number (CN) for n-alkanes. Panel (a) shows the data for chain-lengths up to C20 and Panel (b) from C15 to C60 (the longest chain-length measured experimentally by Nikitin et al. [2]). Included in Fig. 1 are experimental measurements from three different sources [2e4], three

R.A. Messerly et al. / Fluid Phase Equilibria 425 (2016) 432e442

different prediction models based on the experimental data [2,22,23], molecular simulation results from three different sources [15,17,18,21], and the aforementioned Nath-Vetere and SiepmannVetere results. There are several key conclusions from this figure. First, the Nath-Vetere results in Panel (a) appear to be equally reliable as those reported by Errington et al. and Potoff et al. This is significant because Errington et al. used a much more computationally demanding finite-size scaling Grand Canonical Monte Carlo

435

(GCMC) approach and an Exp-6 model (which reproduces Pv more accurately than the LJ 12-6 model) [17]. Potoff et al. also utilized a more accurate model (the Mie 16-6) with GCMC simulations although they did not correct for finite-size effects [18]. Furthermore, we see a slight improvement when comparing SiepmannVetere to the Pc values reported by Siepmann’s group, which utilized Pv GEMC data and Equation (1). The results in Panel (a) demonstrate that simply changing the data analysis method can produce more accurate results. However, the data in Panel (b) demonstrate that the degree of improvement is not large enough to elucidate the correct trend of Pc for larger nalkanes. For example, the Nath-Vetere result for C36 appears to support Nikitin’s old data whereas the Siepmann-Vetere value at C30 agrees more strongly with the recent Nikitin data. Furthermore, notice that the uncertainties in Pc are too large to offer insight into which prediction model (that of Nikitin et al. [2], Lemmon et al. [22], or Nannoolal et al. [23]) correctly describes the long-chainlength behavior of the n-alkanes. For example, for C24, the NathVetere results seem to support either the Nikitin or Nannoolal model while those of Errington et al. more closely follow the Lemmon model. The situation is even worse for C48 where the results from Errington et al. do not agree well with any prediction model and those of Nath-Vetere only moderately agree with Nannoolal’s model. Notably, the uncertainty for C30 reported by Zhuravlev et al. spans nearly all three prediction models while the Siepmann-Vetere uncertainty is considerably smaller. However, minimizing the uncertainty further is necessary to give statistically-meaningful answers regarding the true behavior of Pc for larger n-alkanes. Doing so requires that the uncertainty in rL, rc, and Tc be reduced e something that changing the analysis method from the Antoine Equation to the Rackett Equation cannot address. In Section 4 we propose an experimental design that will minimize the uncertainty in Pc by minimizing the uncertainty in A, B, rc and Tc. 4. Optimal experimental design to minimize the uncertainty in Tc, rc, and Pc 4.1. General theory of experimental designs In general, an optimal experimental design determines the experimental (or, in this case, the simulation) conditions that provide the most precise estimates of the desired property. Although many different experimental designs exist, here we focus on designs that reduce the uncertainty in one or all parameters obtained from regression. This type of experimental design is ideal for the task at hand because our primary goal is to minimize the uncertainties in Tc, rc, and Pc obtained from Equations (2)e(4). The uncertainties for Tc and rc can be minimized by using a parameterspecific design for Equations (2) and (3) because they appear explicitly as regression parameters. The parameter-specific designs for Tc and rc will be different and, therefore, we refer to these two designs as Tc-optimal and rc-optimal designs. However, since Pc depends upon rL, rc, and Tc it also depends upon parameters A and B through Equation (5). Therefore, to minimize the uncertainty in Pc a so-called D-optimal design is necessary to simultaneously reduce the uncertainty in A, B, rc, and Tc.

Fig. 1. The dependence of Pc on CN. Panels (a) and (b) focus on the lower and higher CN ranges, respectively. Three conflicting experimental data sets are included. Simulation results for the TraPPE, Mie 16-6, and Exp-6 models are as reported by Siepmann’s group (Martin et al. and Zhuravlev et al.), Potoff et al. and Errington et al. respectively. The Nath-Vetere and Siepmann-Vetere results are also shown. Panel (b) includes three different prediction models for the long-chain-length behavior. Error bars correspond to 95% confidence intervals.

4.1.1. D-optimal design A D-optimal design is a type of experimental design that minimizes the parameter imprecision by determining the optimal values for the independent variables [24]. The general approach for a D-optimal design is to find the independent variable values that maximize the determinant of the information matrix (which will be explained later). The complete theoretical derivation of the D-

436

R.A. Messerly et al. / Fluid Phase Equilibria 425 (2016) 432e442

optimal design is beyond the scope of this paper, and readers are referred to the literature for the complete analysis [24]. For the present purposes it is sufficient to begin by defining the covariance matrix of the parameters as

 1 VðbÞ ¼ F T F s2

(8)

where V is the covariance matrix, b is a vector of parameters (in our case A, B, rc, and Tc), F is the “F-matrix”, s2 is the inherent variance in the data, and FTF is the information matrix. The F-matrix is obtained as follows. For a function, f, that depends upon the vector of parameters, b, and the independent variable, x, the elements for the Fmatrix are defined as

Fi;j ≡

vf   vbj xi

(9)

where i and j correspond to the respective rows and columns in the   F-matrix and  denotes evaluating the expression at x ¼ xi. xi Since a D-optimal design aims to reduce the uncertainty in all the parameters it must decrease each element of the covariance matrix. In the assumption of constant s2, the only way to minimize V(b) is to minimize (FTF)1. Since the inverse of a matrix is inversely proportional to the determinant  T of  that matrix, minimization of F F . Hence the “D” refers to the (FTF)1 requires  T maximizing  determinant F F  that is key to this analysis. In the case of a non-constant variance, the s2 term in Equation (8) must be included in the information matrix. This is done by including a weighting matrix (W) such that the “weighted information matrix” for non-constant variance can be expressed as

F T W T WF

(10)

where W is a diagonal matrix with Wi,i ¼ si. Since WF effectively acts as a new F-matrix, we define this as the “weighted F-matrix”. The covariance matrix, V(b), for non-constant variance is obtained from the inverse of Equation the D-optimal  (10). Therefore,  design is obtained by maximizing F T W T WF . As s is included in the weighted F-matrix, a model for s is necessary to develop an experimental design with a non-constant variance. Normally, this requirement might hinder a D-optimal design because an error model is not usually available a priori. However, we recently reported an error model for GEMC simulations of n-alkanes that is accurate for a large range of compound sizes [5]. 4.1.2. Parameter-specific design An experimental design that is similar to a D-optimal design is a parameter-specific design, that is, an experimental design that minimizes the uncertainty in a specific parameter rather than the overall uncertainty in all of the parameters. This can be accomplished by minimizing the variance of the element of the covariance matrix, V(b), that corresponds to that parameter. In other words, to minimize the uncertainty in the parameter that corresponds to the jth column of F we minimize the jth main diagonal element of (FTF)1 or of (FTWTWF)1. 4.2. GEMC experimental design In the specific case of performing an experimental design for GEMC simulations, the two equations used in regression are the law of rectilinear diameters and the density scaling law. Between these two coupled equations there are four parameters (A, B, rc, and Tc) and so it is common to think that four different simulation temperatures are required. However, since each simulation renders two

data points (rr and rs) two temperatures are sufficient. Therefore, our D-optimal design will provide the two temperature values that  will maximize F T W T WF . We will also develop an experimental design to minimize the uncertainty in either Tc or rc (since these are the only two parameters with physical meaning in Equations (2) and (3)) by minimizing the respective V(b) element (see Equation (8)). It is important to reiterate that the Pc-optimal design is the same as the D-optimal design since minimizing the uncertainty in Pc obtained from Equation (4) requires minimizing the uncertainty in rL, rc, and Tc. If rL is obtained from Equation (5) then rL, and subsequently Pc, also depends upon the uncertainties in A and B. Therefore, the Pc-optimal design should be the same as the Doptimal design that minimizes the uncertainty in A, B, rc, and Tc. 4.2.1. GEMC temperature constraints It is common in an experimental design that the optimal value for an independent variable is unattainable due to physical or practical limits. This is the case for GEMC since there are limitations that affect the range of feasible temperatures for molecular simulations. The lowest feasible temperature is limited by the difficulty of a successful particle insertion step into the liquid phase. A successful insertion move is less probable when the liquid density and/ or the molecular size increases. If too few insertion moves are accepted the two simulation boxes representing the liquid and vapor phases will not be in chemical equilibrium [6]. For our purposes, the lowest possible TR was assumed to be 0.7. On the other extreme, the highest possible TR was assumed to be 0.95. This upper constraint is imposed for two reasons. First, as the critical point is approached the two simulation boxes often alternate between the vapor and liquid phases [25], so that very long simulation times and more advanced analysis methods are required to obtain reliable estimates of rL and rv [7]. Second, near the critical point a crossover from Ising-like to mean field behavior causes the optimal value of b in Equation (3) to change [6]. Fortunately, this is not observed for GEMC simulations of small n-alkanes below 0.95 Tc [15]. We will show in Section 4.2.2 that when the variance is assumed to be constant with respect to temperature, the D-optimal design recommends performing simulations at these two extreme values to minimize the uncertainty in the prediction of Tc, rc, and Pc. However, in Section 4.2.3 we will also show that when the variance is not considered constant, as is the case with GEMC simulations, it is undesirable to perform simulations at the upper limit because the inherent error in the GEMC results increases exponentially with respect to reduced temperature. Therefore, the somewhat arbitrary choice of 0.95 as the upper bound does not affect any of the conclusions of this work. By contrast, the value of the lower bound does affect the experimental design. There are several cases where the TR lower bound of 0.7 may not be practical and a higher value will be necessary. For example, particle insertion is more difficult with non-linear molecules, such as branched or cyclic compounds. In Section 6.1 we demonstrate how to account for variability in this lower bound value. 4.2.2. Constant variance A good starting point for a discussion of the proposed experimental design is the case of constant variance. In other words, if the uncertainty in the GEMC data were constant with respect to temperature what would the D-optimal design look like? The answer to this question is best understood by linearizing both the law of rectilinear diameters and the density scaling law. The law of rectilinear diameters can be rewritten as a straight line according to

rr ¼ ðrc þ ATc Þ  AT ¼ b1 þ b2 T

(11)

where b1 ¼ rc þ ATc and b2 ¼ A. Similarly, the linearized form of

R.A. Messerly et al. / Fluid Phase Equilibria 425 (2016) 432e442

437

the density scaling law gives 1

1

1

rsb ¼ Bb Tc  Bb T ¼ b1 þ b2 T 1

(12) 1

where b1 ¼ Bb Tc and b2 ¼ Bb . The D-optimal design for a straightline model suggests simulating at the lower and upper constraints (see Supplementary Material in Appendix A). Therefore, it can be expected that the D-optimal design with constant variance supports performing simulations at reduced temperatures of 0.7 and 0.95. This claim is now substantiated. The first step in a D-optimal design is to construct the F-matrix. In this case, we must construct the F-matrix for the non-linearized forms of Equations (2) and (3). When two equations are coupled together, as are Equations (2) and (3), the definition for the F-matrix is not as straightforward as for a single equation. The F-matrix for the law of rectilinear diameters and the density scaling law coupled together is

3 vrr  vrr  vrr  vrr  6 vrc T1 vA T1 vTc T1 vB T1 7 7 6 7 6    7 6 vrr  v r v r v r r r r     7 6 6 vr T2 vA T2 vTc T2 vB T2 7 7 6 c F≡6 7    7 6 vrs  v r v r v r s s s 7  6 6 vr T1 vA T1 vTc T1 vB T1 7 7 6 c 7 6    5 4 vrs  v r v r v r s s s  vrc T2 vA T2 vTc T2 vB T2 2 3 1 Tc  T1 A 0 6 1 Tc  T2 7 A 0 7 ¼6 40 0 bBðTc  T1 Þb1 ðTc  T1 Þb 5 0 0 bBðTc  T2 Þb1 ðTc  T2 Þb 2

  Fig. 2. A contour plot of F T F  for the coupled law of rectilinear diameters and density scaling law with constant variance (see Equation (14)). The best experimental design condition is TR1 ¼ 0:7 and TR2 ¼ 0:95.

2

(13)

where T1 and T2 are the lower and upper temperatures, respectively. From the F-matrix we can express the determinant of the information matrix as

jF T F j ¼





bBTc2b TR1  TR2

2  b1  b1 1  TR1 1  TR2

2

(14)

where TR1 and TR2 are the reduced temperature values that correspond to the lower and upper temperatures, respectively. From this expression it may not be obvious which values of TR1 and TR2  F T F . For this reason, in Fig. 2 we present a contour plot maximize  T  of F F  with respect  to TR1 and TR2 . The blue region corresponds to the maximum in F T F  and thus the optimal set of reduced temperature values. Notice that these plots are symmetric about the line TR1 ¼ TR2 because the assignment of TR2 > TR1 is arbitrary. As seen in Fig. 2, if the variance in the simulation results is constant with respect to temperature the resulting D-optimal design suggests performing simulations at the two extreme TR values, 0.7 and 0.95. Notice that this is already very different from most GEMC studies since the common assumption is that the low temperature data are not as useful when extrapolating to the critical point. Though the derivation is not provided, a Tc-optimal design or rc-optimal design also supports using 0.7 and 0.95 for the TR values when a constant variance is assumed.

4.2.3. Non-constant variance Previous work has shown that the variance in GEMC simulation data is not constant with respect to temperature [5]. To account for the non-constant variance the information matrix must include the weighted F-matrix (see Equation (10)). For Equations (2) and (3) the weighted F-matrix becomes

1 6 sr ðT1 Þ 6 6 6 1 6 6 sr ðT2 Þ 6 WF ¼ 6 6 6 6 0 6 6 6 4 0

Tc  T1 sr ðT1 Þ

A sr ðT1 Þ

Tc  T2 sr ðT2 Þ

sr ðT2 Þ

A

0

bBðTc  T1 Þb1 ss ðT1 Þ

0

bBðTc  T2 Þb1 ss ðT2 Þ

3 0

7 7 7 7 7 0 7 7 7 b7 ðTc  T1 Þ 7 7 ss ðT1 Þ 7 7 7 ðTc  T2 Þb 5 ss ðT2 Þ

(15)

where sr(Ti) and ss(Ti) are the standard deviations evaluated at Ti for rr and rs, respectively. For n-alkanes, sr(Ti) and ss(Ti) are adequately described by

s ¼ b0 þ b1 eb2 TR

(16)

where bi are fitting parameters specific to sr and ss [5]. We will see that this exponential dependence in sr and ss with respect to TR is significant for the experimental design. The D-optimal design for Equations (2) and (3) is obtained by   finding the T1 and T2 values that maximize F T W T WF  where WF is found from Equation (15). The parameter-specific design for rc and 1

Tc are the T1 and T2 values that minimize ðF T W T WFÞ1;1 and 1

ðF T W T WFÞ3;3 , respectively. Since analytical expressions for   T T F W WF  and ðF T W T WFÞ1 are cumbersome and provide little   1 insight, we have provided contour plots of F T W T WF , 1 , ðF T W T WFÞ1;1

and

1 1 ðF T W T WFÞ3;3

in Figs. 3e5, respectively. We have plotted

1 1 ðF T W T WFÞj;j

so that, in each case, the optimal design is found by maximizing the appropriate quantity. Therefore, the D-optimal, rc-optimal, and Tcoptimal designs correspond to the blue region in Figs. 3e5, respectively. There are two main implications from these figures. First, we see that all three designs support having TR1 at the lower constraint of 0.7. Second, in all three designs the value of TR2 is far removed from the upper constraint of 0.95. Specifically, the D-optimal, rc-optimal, and Tc-optimal designs result in a TR2 value of 0.85, 0.84, and 0.89, respectively.

438

R.A. Messerly et al. / Fluid Phase Equilibria 425 (2016) 432e442

  Fig. 3. A contour plot of F T W T WF  where WF is defined by Equation (15). The Doptimal experimental design condition corresponds to the maximum at TR1 ¼ 0:7 and TR2 ¼ 0:85.

Fig. 5. A contour plot of

1 1 ðF T W T WFÞ3;3

where WF is defined by Equation (15). The Tc-

optimal experimental design condition corresponds to the maximum at TR1 ¼ 0:7 and TR2 ¼ 0:89.

rc-optimal, and Tc-optimal designs) with two designs that are more representative of those found in the literature [16,26]. The errors in the three critical constants for each of the five methods are presented in Fig. 6. The “percent uncertainty” in Fig. 6 corresponds to the two-sided 95% confidence interval divided by the best estimate of the critical property multiplied by 100%. The uncertainties displayed in Fig. 6 were obtained by performing a Monte Carlo sampling propagation of errors analysis for rc, Tc, and Pc. For each design, Pc is obtained from Equation (4). The standard deviations needed for the Monte Carlo analysis came from the error models developed in our previous study to predict sr and ss. As such, the quantitative values in Fig. 6 are only valid when the simulation conditions are the same as those used to develop the error models, namely, a GEMC simulation with 200 molecules, roughly 20% in the

Fig. 4. A contour plot of

1 1 ðF T W T WFÞ1;1

where WF is defined by Equation (15). The rc-

optimal experimental design condition corresponds to the maximum at TR1 ¼ 0:7 and TR2 ¼ 0:84.

These three designs are very different from conventional wisdom that places less emphasis on the low temperature regime. Indeed, these results seem counter intuitive since one might think that the closer the simulations are to the critical point the more beneficial they will be due to reducing the degree of extrapolation needed. However, because the error increases exponentially with respect to TR for GEMC orthobaric densities, simulations near the upper limit actually become less beneficial. We will substantiate these conclusions in the following section.

5. Comparison of experimental design with literature In this section, we compare the uncertainties in Tc, rc, and Pc using the three experimental designs explained above (D-optimal,

Fig. 6. A comparison between the experimental designs developed in this work (1e3) and two common designs found in the literature (4e5). Refer to the text for the temperatures used for designs 1 (D-optimal), 2 (rc-optimal), 3 (Tc-optimal), 4 (pre-2013 literature), and 5 (post-2013 literature). Pc was obtained using Equation (4) in each case. The percent uncertainty for rc, Tc, and Pc are presented at the 95% confidence level.

R.A. Messerly et al. / Fluid Phase Equilibria 425 (2016) 432e442

vapor phase, and the same move probabilities (see Ref. [5] for more details). That being said, it is expected that the qualitative nature of Fig. 6 is maintained for different system specifications. Each “literature” design consists of four equally spaced temperatures between a maximum temperature of TR ¼ 0.95 and a minimum temperature. The two “literature” cases differ in the value used for the minimum temperature e TR ¼ 0.8 or 0.9. The lower TR of 0.8 is more indicative of what the common practice was prior to 2013 while the value of 0.9 follows the more recent practice of only including simulations with TR  0:9 in the regression to the law of rectilinear diameters and the density scaling law. The pre2013 literature design uses only a single simulation at each temperature, whereas the post-2013 literature design uses eight replicates at each temperature [16,26]. By contrast, the three experimental designs developed in this work perform two replicates at two temperatures. Therefore, the uncertainties displayed in Fig. 6 for the D-optimal, rc-optimal, and Tc-optimal designs use the same number of simulations as the pre-2013 literature and only 18 as many simulations as the post-2013 literature. If more replicate simulations are performed at each temperature the uncertainties will be less than those reported in Fig. 6. Several significant conclusions may be drawn from Fig. 6. First, the D-optimal design performed the best for Pc (producing the lowest uncertainty) and was nearly the best for rc and Tc. The rc and Tc performance is surprising as Designs 2 and 3 specifically minimize the errors in these two parameters. Second, Fig. 6 shows that the uncertainty in Tc is nearly the same for all five of the experimental designs while the uncertainty in rc and Pc are very dependent on the type of experimental design. Third, these results are consistent with what Panagiotopoulos originally reported when he stated that “critical densities obtained from the rectilinear diameter rule are subject to significantly larger uncertainties than the critical temperatures” [27]. Finally, it is worth restating that the post-2013 literature uses eight times as many simulations as the D-optimal, rc-optimal, and Tc-optimal experimental designs but the uncertainty in rc and Pc are still much greater. This demonstrates that to lower the uncertainty from the law of rectilinear diameters, it is better to use lower temperatures. In conclusion, the D-optimal experimental design greatly improves the uncertainty in Pc and rc without sacrificing precision in Tc and demands fewer computational resources than the current practice. Therefore, our recommendation is to perform several simulations at reduced temperatures of 0.7 and 0.85, the D-optimal conditions. Although in some cases Tc might be the only parameter of interest, the improvement in Tc when using the Tc-optimal design is minimal while the depreciation in rc and Pc is significant. The rcoptimal design is very similar to the D-optimal design and, therefore, it is a reasonable alternative. In any case, it is more beneficial to have replicates at two temperatures rather than the common practice of spreading out the replicates among several equally dispersed temperature values. Furthermore, simulations near the critical point require longer computer time to obtain reliable averages due to the high degree of fluctuations. The D-optimal design suggests that these simulations are not necessary for precise estimates of the critical point constants.

439

that of the lower constraint but the value of TR2 must be obtained by   maximizing F T W T WF , T T1 1 , or T T1 1 . To demonstrate ðF W WFÞ1;1

ðF W WFÞ3;3

the effect of a different value for the lower constraint, in Fig. 7 we have plotted the dependence of TR2 on TR1 for the D-optimal, rcoptimal, and Tc-optimal designs. The most important point from this figure is that the optimal upper temperature value is typically much lower than the upper constraint of 0.95 Tc regardless of the value of TR1 . The significance of this will be shown in Section 6.2 concerning finite-size effects. In brief, since finite-size effects increase near the critical point, selecting the lowest practical value for TR1 will reduce TR2 and thereby help to mitigate finite-size effects. With a different value for the lower constraint, and thus for TR1 , the uncertainty in Tc, Pc, and rc will change. In Fig. 8 we plot the uncertainty in Tc, Pc, and rc with respect to TR1 . These uncertainties were obtained in the same was as those in Fig. 6, namely from a Monte Carlo sampling propagation of errors. Since we propose use of the D-optimal design, the uncertainties shown in Fig. 8 correspond to the D-optimal design. In other words, for each value of TR1 the value for TR2 was determined from the D-optimal design curve in Fig. 7. Notice that the uncertainties in Tc, Pc, and rc increase with increasing values of TR1 . The increase in uncertainty is most pronounced for Pc while for Tc the error remains relatively small. 6.2. Implications on finite-size effects Finite-size effects arise as the ratio of molecules on the periodic boundary of a box to total molecules increases [6]. This effect is prevalent in simulations of larger compounds (CN > 26) because as the size of the compound increases the number of molecules required to alleviate finite-size effects increases while the amount of molecules that can be simulated in a feasible length of time decreases considerably. Also, it is a well established fact that finitesize effects become more significant as T approaches Tc [6,16,27]. Fortunately, the D-optimal design suggests that it is unnecessary to perform simulations so close to the critical point.

6. Considerations 6.1. Possible modifications to experimental designs In Section 4.2.1 we explained why we assumed that the lower temperature constraint was TR ¼ 0.7. However, in some cases the value of this constraint may change which would not only cause TR1 to change but also TR2 . The optimal value for TR1 will always equal

Fig. 7. The dependence of TR2 on TR1 for the three different experimental designs. The value for the lower reduced temperature is determined by the GEMC limitation for particle insertion of high density liquids.

440

R.A. Messerly et al. / Fluid Phase Equilibria 425 (2016) 432e442

Fig. 8. The effect of TR1 upon the uncertainty in rc, Tc, and Pc when using the D-optimal design (see Fig. 7 for the values used for TR2 ). Pc was obtained using the Vetere approach (Equation (4)). The percent uncertainty for rc, Tc, and Pc are presented at the 95% confidence level.

The initial objective for performing an experimental design was simply to reduce the statistical uncertainty of the parameters in Equations (2) and (3) and thereby reduce the statistical uncertainties in Tc, rc, and Pc. In other words, typically an experimental design only reduces the uncertainty associated with regressing the parameters. However, as explained previously, the experimental designs presented can also be used to reduce finite-size effects. Although it is unclear whether the D-optimal design significantly reduces finite-size effects, it certainly does not contribute to them. Indeed, the amount to which finite-size effects are decreased still depends upon the size of the compound and number of molecules used in simulation. The error model mentioned previously was developed using a relatively small number of molecules (200) so that it should be accurate for the size of simulations most common for large molecules. Therefore, the results presented previously in Figs. 6e8 are quantitatively valid for the size of systems typically used for larger molecules. However, future work is necessary to verify that finitesize effects are reasonably small for long chain-length compounds when performing GEMC simulations with 200 molecules at the Doptimal design conditions.

validation simulations at temperatures other than those from the D-optimal design. For example, this validation can be achieved easily by performing a single simulation at TR ¼ 0.75,0.8,0.9 and 1 0.95. Specifically, if a plot of rr and rbs with respect to T yields a straight line then Equations (2) and (3) are valid over the temperature range in question. We have found this to be the case for TR as low as 0.7 for n-alkanes. Another issue that has not been discussed in detail until now is the fact that TR is not known a priori because TR depends upon the value of Tc. This means that determining the temperature for which TR ¼ 0.7 and 0.85 would be an iterative process where a new estimate for Tc would be obtained after each iteration. At first glance this appears to be problematic; however, this is of little practical importance as Tc can be estimated reasonably well with several prediction models [28]. So the question is simply how sensitive is the experimental design to variations in TR? Figs. 3e5 show how much the uncertainties depend on locating the exact optimal values. When the contours are closer together the uncertainty increases more rapidly with respect to TR. This can be visualized more easily by slicing through a contour at a fixed value of TR1 . Fig. 9 is provided to demonstrate how much the 95% confidence intervals in Tc, rc, and Pc increase with respect to TR2 (TR1 is assumed to be 0.7). The results in Fig. 9 were obtained in the same way as those in Fig. 6. This figure demonstrates that locating the exact value for TR2 is not crucial because the wells for rc, Tc, and Pc are fairly wide. When we compare Fig. 9 with Fig. 8, we see that having TR1 as low as feasibly possible is more important than having TR2 at the exact design conditions. Another consideration is that when the Rackett approach was originally proposed by Vetere he stated that it was inadequate for compounds with CN > 12. However, Vetere’s article was published in 1988, two years before Teja’s experimental critical point data became available for n-alkanes ranging from C12-C18 [4]. When we repeated Vetere’s analysis with Teja’s data the results did not show any depreciation with increasing chain length. However, the fact remains that the Rackett exponent, g, is simply an empirical parameter. Vetere demonstrated that g is not a universal constant

6.3. Possible limitations Although we have discussed some significant advantages for this methodology, in this section we discuss some possible limitations. For example, the proposed D-optimal design does not consider whether the law of rectilinear diameters and the density scaling law are valid at low temperatures. In fact, one of the primary reasons why the practice found in the most recent literature excludes low temperature data is to avoid potentially poor performance of the law of rectilinear diameters [7]. Therefore, it is important to demonstrate that Equations (2) and (3) are reliable over the entire temperature range where simulations are performed. For this reason, we recommend performing a few

Fig. 9. The effect of TR2 upon the uncertainty in rc, Tc, and Pc (TR1 ¼ 0:7). Pc was obtained using Equation (4). The percent uncertainty for rc, Tc, and Pc are presented at the 95% confidence level.

R.A. Messerly et al. / Fluid Phase Equilibria 425 (2016) 432e442

for all compounds and, in fact, it is not even constant with respect to temperature for a given compound [29]. A sensitivity analysis of Equation (4) predicts that Pc will depend less on g when 
1  TTci z1. Therefore, based solely on the Rackett equation, Ti in Equation (4) should be as low as possible to mitigate the effect of g on the estimate for Pc. Fortunately, our experimental design is consistent with this concept of using temperatures far removed from Tc. In addition, the uncertainty in Pc is lowest when Ti is much lower than Tc because the uncertainty in rL decreases exponentially with respect to temperature [5]. Furthermore, the algorithm for estimating the uncertainty in Pc outlined in Section 2 and the results presented in Figs. 6, 8 and 9 correspond only to the statistical uncertainty. In other words, we have neglected possible systematic deviations due to g and b deviating from 27 and 0.326, respectively. Fortunately, our experience is that possible systematic deviations are on the order of 1%, which is usually well within the 95% statistical uncertainties for larger compounds. This was validated by varying g between 0.275 and 0.295 (the range Rackett originally reported [20]) and b between 0.32 and 0.33. Finally, as we mentioned in Section 4.2, the error model presented in Equation (16) is pivotal to the experimental designs. Since this error model was developed assuming a constant set of Monte Carlo move probabilities it is reasonable that different sets of moves would yield a different dependence of s on TR. Indeed, some studies have attempted to determine the dependence of s on the Monte Carlo move probabilities [30]. Therefore, an alternative approach to what we have presented for reducing the uncertainties in Tc, rc, and Pc is to optimize the Monte Carlo move probabilities at different temperatures. However, in Section 4.2.2, we observed that even if the error were constant with respect to temperature the experimental designs would still support performing simulations at the lower bound. Future attention should be given to optimizing the Monte Carlo move probabilities in both the high and low temperature regimes.

441

relatively small for the experimental designs found in the literature, the uncertainty in rc can be quite large. For this reason, we have developed an experimental design that minimizes the uncertainty in Tc, rc, and Pc. This experimental design has several significant implications that run contrary to the common practice found in the simulation literature. First, it is essential to perform simulations at the lowest temperature feasible. In the literature it is quite common to undervalue the low temperature regime and over emphasize the high temperature regime. Second, only two temperatures are necessary for accurate predictions of the critical constants. Third, it is better to perform several replicate simulations at the two optimal temperatures rather than running simulations at several temperatures. Finally, it is not necessary to perform simulations at the upper temperature limit because the inherent uncertainty associated with the simulation results in this region increases exponentially with respect to temperature. This last conclusion is particularly significant because the required simulation length increases substantially near the critical point. Our study has demonstrated that this additional computational cost is actually not as beneficial for predicting the critical point as most have assumed. Furthermore, finite-size effects are reduced by performing simulations farther away from the critical point. These results are very important for simulations of large compounds where the number of molecules used in simulation is often small enough that finite-size effects are expected to be prevalent. Therefore, our D-optimal design should reduce the need for finitesize corrections and permit accurate estimates of the critical point for large compounds with small simulation sizes. The validation of this claim is left as future work. Acknowledgements The authors are grateful to the Design Institute for Physical Properties (DIPPR 801) for funding. Appendix A. Supplementary data

7. Outline of methodology

Supplementary data related to this article can be found at http:// dx.doi.org/10.1016/j.fluid.2016.06.041.

In this section we provide an outline of the methodology to assist the reader:

References

1. Estimate Tc from previous simulation results or a prediction model. 2. Determine the lowest practical simulation temperature (T1) with T ¼ 0.7  Tc as an initial estimate. 3. Obtain TR2 from the D-optimal design curve in Fig. 7 and calculate T2 ¼ TR2  Tc. 4. Perform several replicate GEMC simulations at T1 and T2 to obtain rL and rv. 5. Regress Equations (2) and (3) to the GEMC data to obtain A, B, rc, and Tc. 6. Calculate rL from Equation (5) with Ti ¼ T1. 7. Use Equation (4) to calculate Pc from rc and Tc obtained in Step 5 and rL obtained in Step 6. 8. Follow the uncertainty algorithm outlined in Section 2.

8. Conclusions In this work we have demonstrated that an alternative method, namely the Vetere approach, has many advantages when predicting Pc from GEMC for larger compounds. However, this method is very sensitive to the uncertainty in rc. Although the uncertainty in Tc is

[1] E.A. Guggenheim, The principle of corresponding states, J. Chem. Phys. 13 (7) (1945) 253e261, http://dx.doi.org/10.1063/1.1724033. [2] E.D. Nikitin, A.P. Popov, Critical temperatures and pressures of c-40, c-44, and c-60 normal alkanes measured by the pulse-heating technique, Fluid Phase Equilib. 379 (2014) 191e195, http://dx.doi.org/10.1016/j.fluid.2014.07.017. [3] E.D. Nikitin, P.A. Pavlov, A.P. Popov, Vapour-liquid critical temperatures and pressures of normal alkanes with from 19 to 36 carbon atoms, naphthalene and m-terphenyl determined by the pulse-heating technique, Fluid Phase Equilib. 141 (1e2) (1997) 155e164, http://dx.doi.org/10.1016/S03783812(97)00202-1. [4] A.S. Teja, R.J. Lee, D. Rosenthal, M. Anselme, Correlation of the critical properties of alkanes and alkanols, Fluid Phase Equilib. 56 (1990) 153e169, http:// dx.doi.org/10.1016/0378-3812(90)85100-O. [5] R.A. Messerly, R.L. Rowley, T.A. Knotts, W.V. Wilding, An improved statistical analysis for predicting the critical temperature and critical density with gibbs ensemble monte carlo simulation, J. Chem. Phys. 143 (10) (2015) 104101, http://dx.doi.org/10.1063/1.4928865. [6] D. Frenkel, B. Smit, Understanding Molecular Simulation : from Algorithms to Applications, in: Computational science series, Academic Press, San Diego, 2002. [7] M. Dinpajooh, P. Bai, D.A. Allan, J.I. Siepmann, Accurate and precise determination of critical properties from gibbs ensemble monte carlo simulations, J. Chem. Phys. 143 (11) (2015) 114113, http://dx.doi.org/10.1063/1.4930848. [8] A. Vetere, Estimation of critical pressures by the rackett equation, Chem. Eng. Sci. 44 (4) (1989) 791e795, http://dx.doi.org/10.1016/0009-2509(89)85252-2. [9] J.J. de Pablo, Simulation of phase-equilibria for chain molecules, Fluid Phase Equilib. 104 (1995) 195e206, http://dx.doi.org/10.1016/0378-3812(94) 02649-L. [10] A.C. Biswas, The law of rectilinear diameter for the liquid-gas phase transition, Pramana 1 (2) (1973) 109e111.

442

R.A. Messerly et al. / Fluid Phase Equilibria 425 (2016) 432e442

[11] J.S. Rowlinson, Liquids and Liquid Mixtures, in: Butterworths Monographs in Chemistry, Butterworth Scientific, London ; Boston, 1982. [12] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Clarendon Press ; Oxford University Press, Oxford England New York, 1987. [13] J. Chang, S.I. Sandler, Interatomic lennard-jones potentials of linear and branched alkanes calibrated by gibbs ensemble simulations for vapor-liquid equilibria, J. Chem. Phys. 121 (15) (2004) 7474e7483, http://dx.doi.org/ 10.1063/1.1792572. [14] W.L. Jorgensen, J.D. Madura, C.J. Swenson, Optimized intermolecular potential functions for liquid hydrocarbons, J. Am. Chem. Soc. 106 (22) (1984) 6638e6646, http://dx.doi.org/10.1021/Ja00334a030. [15] M.G. Martin, J.I. Siepmann, Transferable potentials for phase equilibria. 1. united-atom description of n-alkanes, J. Phys. Chem. B 102 (14) (1998) 2569e2577 doi:10.1021/Jp972543þ. [16] S.K. Nath, F.A. Escobedo, J.J. de Pablo, On the simulation of vapor-liquid equilibria for alkanes, J. Chem. Phys. 108 (23) (1998) 9905e9911, http:// dx.doi.org/10.1063/1.476429. [17] J.R. Errington, A.Z. Panagiotopoulos, A new intermolecular potential model for the n-alkane homologous series, J. Phys. Chem. B 103 (30) (1999) 6314e6322, http://dx.doi.org/10.1021/Jp990988n. [18] J.J. Potoff, D.A. Bernard-Brunel, Mie potentials for phase equilibria calculations: applications to alkanes and perfluoroalkanes, J. Phys. Chem. B 113 (44) (2009) 14725e14731, http://dx.doi.org/10.1021/jp9072137. [19] A. Vetere, The riedel equation, Ind. Eng. Chem. Res. 30 (11) (1991) 2487e2492, http://dx.doi.org/10.1021/Ie00059a020. [20] H.G. Rackett, Equation of state for saturated liquids, J. Chem. Eng. Data 15 (4) (1970) 514e517, http://dx.doi.org/10.1021/Je60047a012. [21] N.D. Zhuravlev, M.G. Martin, J.I. Siepmann, Vapor-liquid phase equilibria of triacontane isomers: deviations from the principle of corresponding states, Fluid Phase Equilib. 202 (2) (2002) 307e324, http://dx.doi.org/10.1016/ S0378-3812(02)00137-1.

[22] E.W. Lemmon, A.R.H. Goodwin, Critical properties and vapor pressure equation for alkanes cnh2nþ2: normal alkanes with n ¡¼ 36 and isomers for n¼4 through n¼9, J. Phys. Chem. Ref. Data 29 (1) (2000) 1e39, http://dx.doi.org/ 10.1063/1.556054. [23] Y. Nannoolal, J. Rarey, D. Ramjugernath, Estimation of pure component properties part 2. estimation of critical property data by group contribution, Fluid Phase Equilib. 252 (1e2) (2007) 1e27, http://dx.doi.org/10.1016/ j.fluid.2006.11.014. [24] J. Kiefer, Optimum experimental-designs, J. R. Stat. Soc. Ser. B Stat. Methodol. 21 (2) (1959) 272e319. [25] B. Smit, S. Karaborni, J.I. Siepmann, Computer-simulations of vapor-liquid phase-equilibria of n-alkanes, J. Chem. Phys. 102 (5) (1995) 2126e2140, http://dx.doi.org/10.1063/1.469563. [26] B. Eggimann, P. Bai, A. Bliss, Q. Chen, T. Chen, A. Corest-Morales, E. Fetisov, E. Haldoupis, D. Harwood, R. Lindsey, T. Arachchi, M. Shah, H. Stern, K. Struk, J. Sung, A. Sunnarborg, B. Xue, J. I. Siepmann, T-ua no. 2 ethane. URL http:// www.chem.umn.edu/groups/siepmann/trappe/. [27] A.Z. Panagiotopoulos, Molecular simulation of phase coexistence - finite-size effects and determination of critical parameters for 2-dimensional and 3dimensional lennard-jones fluids, Int. J. Thermophys. 15 (6) (1994) 1057e1072, http://dx.doi.org/10.1007/Bf01458815. [28] K.G. Joback, R.C. Reid, Estimation of pure-component properties from groupcontributions, Chem. Eng. Commun. 57 (1e6) (1987) 233e243, http:// dx.doi.org/10.1080/00986448708960487. [29] A. Vetere, Again the rackett equation, Chem. Eng. J. Biochem. Eng. J. 49 (1) (1992) 27e33, http://dx.doi.org/10.1016/0300-9467(92)85021-Z. [30] A.D.C. Morales, I.G. Economou, C.J. Peters, J.I. Siepmann, Influence of simulation protocols on the efficiency of gibbs ensemble monte carlo simulations (vol. 39, pg 1135, 2013), Mol. Simul. 39 (14e15) (2013), http://dx.doi.org/ 10.1080/08927022.2013.853910, 1293e1293.