Adapting SAFT- perturbation theory to site-based molecular dynamics simulation. I. Homogeneous fluids Ahmadreza F. Ghobadi and J. Richard Elliott Citation: The Journal of Chemical Physics 139, 234104 (2013); doi: 10.1063/1.4838457 View online: http://dx.doi.org/10.1063/1.4838457 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/139/23?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Adapting SAFT- perturbation theory to site-based molecular dynamics simulation. II. Confined fluids and vaporliquid interfaces J. Chem. Phys. 141, 024708 (2014); 10.1063/1.4886398 Thermodynamic of fluids from a general equation of state: The molecular discrete perturbation theory J. Chem. Phys. 140, 234504 (2014); 10.1063/1.4882897 Perturbation theory for multipolar discrete fluids J. Chem. Phys. 135, 134511 (2011); 10.1063/1.3646733 Application of a renormalization-group treatment to the statistical associating fluid theory for potentials of variable range (SAFT-VR) J. Chem. Phys. 134, 154102 (2011); 10.1063/1.3570614 Free energy calculation using molecular dynamics simulation combined with the three dimensional reference interaction site model theory. I. Free energy perturbation and thermodynamic integration along a coupling parameter J. Chem. Phys. 133, 044114 (2010); 10.1063/1.3462276

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.220.71.24 On: Mon, 11 Aug 2014 05:20:30

THE JOURNAL OF CHEMICAL PHYSICS 139, 234104 (2013)

Adapting SAFT-γ perturbation theory to site-based molecular dynamics simulation. I. Homogeneous fluids Ahmadreza F. Ghobadi and J. Richard Elliotta) Department of Chemical and Biomolecular Engineering, The University of Akron, Akron, Ohio 44325, USA

(Received 7 August 2013; accepted 19 November 2013; published online 16 December 2013) In this work, we aim to develop a version of the Statistical Associating Fluid Theory (SAFT)-γ equation of state (EOS) that is compatible with united-atom force fields, rather than experimental data. We rely on the accuracy of the force fields to provide the relation to experimental data. Although, our objective is a transferable theory of interfacial properties for soft and fused heteronuclear chains, we first clarify the details of the SAFT-γ approach in terms of site-based simulations for homogeneous fluids. We show that a direct comparison of Helmholtz free energy to molecular simulation, in the framework of a third order Weeks-Chandler-Andersen perturbation theory, leads to an EOS that takes force field parameters as input and reproduces simulation results for Vapor-Liquid Equilibria (VLE) calculations. For example, saturated liquid density and vapor pressure of n-alkanes ranging from methane to dodecane deviate from those of the Transferable Potential for Phase Equilibria (TraPPE) force field by about 0.8% and 4%, respectively. Similar agreement between simulation and theory is obtained for critical properties and second virial coefficient. The EOS also reproduces simulation data of mixtures with about 5% deviation in bubble point pressure. Extension to inhomogeneous systems and united-atom site types beyond those used in description of n-alkanes will be addressed in succeeding papers. © 2013 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4838457] I. INTRODUCTION

Molecular simulation has become an indispensable tool in developing theories to describe fluid properties. Molecular simulation has now advanced to the stage that a fair degree of confidence can be placed in its inferences when compared to experimental data, at least qualitatively, and sometimes quantitatively.1 Hence, it is justifiable to consider molecular simulation as a basis for testing theories in situations where experimental measurements are difficult or inaccessible. Ideally, atomistic detail applied in simulation should be sufficient to characterize the properties of interest and theories should be refined until they are consistent with the details of the simulations. Currently, the leading approach in the prediction of fluid properties is Statistical Associating Fluid Theory (SAFT).2–5 Equations of state (EOSs) based on SAFT, which benefit from Wertheim’s first order Thermodynamic Perturbation Theory (TPT1),6–9 reproduce molecular simulation results and provide an analytical equation that facilitates interpolation and extrapolation of fluid properties on much shorter time scales than full potential molecular simulation. In its original form, and most subsequent forms,10–13 SAFT models assumed that polyatomic species were composed of tangent sphere segments. Although this form was shown to agree reasonably well with simulations of tangent hard sphere chains, significant discrepancies could arise when compared to simulations of more realistic molecular models such as fused soft-sphere chains. One particular implementation of SAFT that ostensibly includes atomistic detail in fused a) Author to whom correspondence should be addressed. Electronic mail:

[email protected]

0021-9606/2013/139(23)/234104/24/$30.00

sphere descriptions of molecules is the SAFT-γ EOS.14, 15 The SAFT-γ model assigns an interaction site to each sitetype in a united-atom (UA) characterization of the molecule. This level of description is consistent with molecular simulations based on transferable force fields. However, characteristic repulsive and attractive adjustable parameters of the SAFT-γ EOS have been fitted to experimental bulk fluid properties in previous work. Therefore, agreement at the atomistic level with molecular simulation is not guaranteed. In fact, we demonstrate below that there can be significant deviations between simulation and theory when reference and perturbation contributions of the EOS are lumped together and fitted to experimental data. In this context, we seek to develop a form of SAFT-γ EOS such that molecular simulations – based on transferable, united atom force fields – can be checked for consistency with EOS predictions. While our objective is a theory of interfacial properties, we find it necessary to first elucidate the details of the SAFT-γ approach in terms of site-based simulations for homogeneous fluids. Once these details are well-defined at the bulk level, we plan to adapt classical density functional theory (DFT)16–19 to address the interfacial properties. DFT involves estimating the Helmholtz energy at intermediate densities between the vapor and the liquid, in other words, at thermodynamically unstable points inside the binodal. This step requires the EOS to be accurate in order for extrapolation to provide meaningful results. We aim to establish that accuracy based on homogeneous (bulk) fluids, then to test the extrapolations based on molecular simulations of inhomogeneous fluids. We apply thermodynamic perturbation theory to investigate the properties of soft repulsive reference fluids for

139, 234104-1

© 2013 AIP Publishing LLC

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.220.71.24 On: Mon, 11 Aug 2014 05:20:30

234104-2

A. F. Ghobadi and J. R. Elliott

J. Chem. Phys. 139, 234104 (2013)

FIG. 1. The definition of reference (solid line) and perturbation potential (dashed line) for SW (a), BH (b), and WCA (c) perturbation theories.

realistic and fused models of polyatomic molecules. In the process, we adopt the Weeks-Chandler-Andersen (WCA)20 description for separating the reference and perturbation contributions of the potential model. The SAFT-γ methodology plays an important role in developing a master curve for each reference fluid that can be applied to multiple temperatures and multiple values of the well-depth parameter (ε). Once the reference fluid is characterized, we turn our attention to the perturbation contributions. In particular, we compute perturbation terms for multiple chain lengths rather than simply extrapolating corresponding contributions of spherical molecules, finding that such extrapolations are not entirely reliable. Consistent with the recent works of Avendano and coworkers,21, 22 we find that inclusion of the third order perturbation term substantially improves accuracy in the critical region, noting that this accuracy is important for interfacial properties. By following a systematic approach in tuning adjustable parameters, the resulting EOS takes force field parameters (σ and ε) as input variables and reproduces the full potential simulation results. We demonstrate the feasibility of this approach by accurately generating Vapor-Liquid Equilibria (VLE) results for several different potential models while performing reference fluid simulations for only one of those models, The remainder of this paper is organized as follows: In Sec. II, background on the perturbation perspective is presented as applied in liquid theories and SAFT approaches. A general hybrid procedure of molecular simulation and theory is also described in this section. Simulation details are given in Sec. III. Characterization of reference and perturbation contributions of the proposed EOS along with the underlying theory is discussed in Sec. IV. The evaluation of EOS parameters is described in Sec. V. Finally, the predictions of the proposed methodology for thermophysical properties of pure homogenous fluids and also mixtures are presented in Sec. VI.

erence (u0 ) and perturbation (u ) moieties: u(r) = u0 (r) + u (r),

(1)

where r is radial distance between two interaction sites. Figure 1 illustrates the common separation schemes for SW and LJ potentials. In this work, we abbreviate the separation of interaction potential as BH TPT when the BarkerHenderson split is applied and WCA TPT when the WCA split is applied. We do not apply the WCA criterion for effective hard sphere diameter20 and intend no confusion with that theory when we refer to WCA TPT. Implementation of a density dependent hard sphere diameter would substantially complicate the equation of state. Perturbation contributions, corresponding to the perturbation potential of Eq. (1), are expressed as the fluctuation of energy around the reference systems using the Zwanzig high temperature expansion:26–28 Apert A3 A1 A2 = + 2 + 3 + ··· , N kB T T T T 1 U  0 , N kB  −1  (U  )2 0 − U  20 , A2 = 2 2!N kB

A1 =

A3 =

(2)

 1  (U  )3 0 − 3(U  )2 0 U  0 + 2U  30 . 3!N kB3

In Eq. (2),  0 denotes an average over the configurations of the reference system, kB is Boltzmann’s constant, N is the total number of molecules, T is temperature, Apert is the perturbation Helmholtz energy and subscripts 1−3 specify first, second, and third order perturbation contributions. U is total potential of the system such that U =

N 

u (rij ).

(3)

iσαβ × uαβ (|r1 − r2 |)dr1 dr2 ,

(6)

where g(r, ρ) ¯ is the RDF evaluated at distance r and an average density ρ¯ because the exact value of RDF is not known for inhomogeneous systems.31–33 With the choice of g(r, ρ) ¯ = 1, the Carnahan-Starling34 EOS for HS, and a WCA scheme for u , Eqs. (5) and (6) form the basis of the SAFT EOS for interfaces known as iSAFT. The iSAFT EOS has been successfully applied to study self-assembly, structure-property relationships, microstructure and inhomogeneous properties of homonuclear polymers,35–37 heteronuclear chains resembling surfactants,38–40 tethered polymers,41 polymer-colloid mixtures,42 diblock copolymers,43 triblock rod-coil amphiphiles,44 associating polyatomic molecules,45 lipid bilayers,46 and patchy colloids.47 The approach has also been extended beyond TPT1 to treat a variety of molecular topologies including ring, branched, and rigid structures.40, 48 In spite of the broad application of the iSAFT EOS, it only provides a coarse-grained characterization of molecules

in the form of hypothetical tangent sphere chains. In order to obtain an accurate representation of bulk fluid properties for real molecules with a soft and fused molecular topology, this level of coarse-graining usually leads to a non-integer number of segments (m). In the framework of DFT, it is not possible to solve for individual segment density profiles with a non-integer m. Consequently, to obtain the microstructure of real molecules, integer numbers of segments must be enforced during adjustable parameter evaluation. That sacrifices accuracy of the EOS in the bulk limit when the tangent sphere basis is recovered.49 Furthermore, the underlying assumption of g(r, ρ) ¯ = 1 in the iSAFT approach is a substantial oversimplification for complex molecular fluids that restricts the accuracy of iSAFT to narrow temperature ranges. Using a simple mean field approximation and ignoring higher order energy fluctuations would also have negative impacts on predictions of critical and derivative properties. In the last two decades, there has been a systematic effort to improve the description of the SAFT model for attractive and repulsive interactions.10, 13 The original SAFT EOS was extended to HS mixtures50 and soft-spheres51 in the framework of BH TPT25 with the dispersion contribution modified to a high temperature expansion obtained from simulation data of LJ4 and SW53, 54 fluids. Johnson et al.55 combined SAFT theory with the previously developed LJ EOS for spheres56 and presented a LJ EOS for chains. Gubbins and coworkers57, 58 incorporated dipole-dipole interactions into Johnson’s EOS for chains resulting in improved VLE description for small polar molecules such as methanol and ethanol. Johnson’s EOS was extended to heteronuclear chains and mixtures by Blas and Vega59 in a variation known as soft-SAFT. Until that stage, all SAFT EOSs, including Kolafa’s EOS for the LJ fluid60 and Banaszak’s EOS for SW61 and LJ62 freely jointed chains, were either based on MF approximation or were developed with no explicit form of chain RDF for attractive interactions. In other words, a polynomial of temperature and density was used for the last term of Eq. (5), which cannot be reformulated in the form of Eq. (6) within the DFT framework. The RDF for segment-segment interactions of chain molecules is lengthy and, from a practical perspective, burdensome to be included in the EOS. Instead, as proposed by Gross and Sadowski,63 one can fit an integral of the form:  I = g(m; r, ρ)u (r)r 2 dr, m m 1  g(rαβ , ρ), g(m; r, ρ) = 2 m α=1 β=1

(7)

with a polynomial of m and ρ, if the perturbation potential (u ) and g(rαβ , ρ) is known. In Eq. (7), g(m; r, ρ) is the average segment-segment interaction for chains formed of m segments at number density ρ. Gross and Sadowski64 used Chiew’s expression for the chain RDF65 along with the LJ potential for u and fitted the integral of Eq. (7) to a polynomial with 18 constants. The corresponding integral for the second order perturbation contribution in the framework of the Local Compressibility Approximation (LCA) led to another 18-constant polynomial of m and ρ. They further opti-

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.220.71.24 On: Mon, 11 Aug 2014 05:20:30

234104-4

A. F. Ghobadi and J. R. Elliott

J. Chem. Phys. 139, 234104 (2013)

mized these 36 universal constants together with componentspecific parameters (m, σ , and ε) to experimental data of pure n-alkanes. The resulting perturbation contributions and an effective hard sphere diameter (dEHS ) that was obtained from a shoulder potential formed the basis for perturbed-chain SAFT EOS. In PC-SAFT formalism, due to the optimization process of the proposed universal constants, one cannot retrieve perturbation contributions by directly inserting Chiew’s RDF in Eq. (6). Without further modifications, this would result in a non-zero surface tension for bulk phases when the EOS is combined with DFT.33 As a solution, the difference between bulk and DFT expressions of dispersion term is locally augmented to the dispersion functional (refer to Eq. (27) of Ref. 33). Furthermore, having a coarse-grained RDF in the form of an average segment-segment quantity hinders the application of PC-SAFT to study the microstructure of inhomogeneous systems. Another possible alternative to include structural information in the DFT approach is factoring the RDF out of integral of Eq. (6) using the Mean Value Theorem (MVT) as 1  g(σαβ , ρ¯ eff ) 2 α=1 β=1  ρα (r1 )ρβ (r2 )uαβ (|r1 − r2 |)dr1 dr2 , × m

Aatt [{ρα }] =

m

|r1 −r2 |>σαβ

(8) where g(σ αβ , ρ) ¯ is r-independent RDF at HS contact evaluated for an effective density. Effective density can be a polynomial function of density with coefficients that are regressed to simulation data so that Eq. (8) reproduces Eq. (6) for all densities. This idea was implemented in the development of the SAFT EOS for chain molecules with attractive potential of variable range66 (SAFT-VR). Although the SAFT-VR perspective can be adapted to a variety of interaction potentials,67, 68 it is mainly developed for SW chains and it has been successfully applied to phase equilibria calculations of pure compounds and mixtures. As a result of using MVT for inclusion of structural effects, SAFT-VR is conveniently combined with DFT to study interfacial tension and density profile of associating and non-associating SW chains and acceptable agreement with molecular simulation results is obtained.31, 69, 70 But, the predictions of the SAFT-VR EOS for surface tension, especially for associating compounds such as alcohols and water,71 are only acceptable when experimental surface tension data are included in the fitting procedure. In those cases, the trade-off is a decrease in the accuracy of VLE prediction. In addition, for compounds where calculated critical properties deviate significantly from experimental values, a shift corresponding to the difference between predicted and experimental critical properties is necessary to achieve reliable surface tension prediction.72 Finally, it should be mentioned that SAFT-VR is a segment-based EOS and, therefore, an investigation into the microstructure of the system at the atomistic level is not possible without a meaningful relation between the coarse-grained and atomistic description of molecules.

In the DFT framework, a detailed investigation of the microstructure of inhomogeneous systems at the atomistic level is feasible only if the Helmholtz energy (Eq. (5)) is formulated for molecules formed of realistic sites rather than segments. There are three major options in the SAFT family of EOSs that provide this type of information: Group Contribution (GC) SAFT,73–75 GC-SAFT-VR,76, 77 and SAFT-γ .14, 15 The advantage of the SAFT-γ approach over other group contribution methods is that the underlying theory of chain molecules is developed for fused heteronuclear chains that are formed with united atom site types. As proposed by Jackson and coworkers,14, 15, 78 a “shape factor” is defined to characterize the proportion of the spherical united atom group that contributes to the macroscopic properties of the chain. In the SAFT-γ EOS, the shape factor, size, and energy parameters of each site type are fitted to the experimental data of vapor pressure and liquid density of non-associating pure compounds. For associating molecules, two more adjustable parameters are defined along with the binary interaction parameters. Having a site-based EOS holds value only if groupcontribution approaches are targeted, or if it is employed to extract properties at the atomistic level for the system of interest. An example of such a property is the density profile of a particular site-type in an inhomogeneous system. For example, it would be interesting to describe in theory how the density profile of an OH site type differs from that of a CH3 site type at the vapor liquid interface of n-hexanol, then show that simulations yield the same description quantitatively. To reliably obtain such information, acceptable agreement between molecular simulation and theory for reference and perturbation contributions of the EOS is essential. We find that the SAFT-γ formulation provides a level of detail that can be directly validated with molecular simulation. B. Simulation/SAFT/TPT hybrid procedure

Perturbation theory, as described by Eqs. (1)–(4), provides a basis to estimate attractive interactions using statistical averaging over the configurations of the reference system. The approach can be adopted by SAFT methodology to evaluate Aatt in Eq. (5). To do so, one needs to simulate the reference potential at several densities to cover the entire phase diagram as shown by Barker and Henderson.79 The simulation time is significantly less than that of the full potential simulations because the cut-off distance is very short for the reference potential. If the reference system is athermal, as in the case of the SW separation scheme (Figure 1(a)), the simulations can be carried out at a single temperature. Otherwise, the temperature dependency of the generated configurations should be examined. Furthermore, for soft potentials, an effective HS diameter must be introduced to facilitate the use of well-developed HS theory in describing the properties of the reference system.80 Then, perturbation contributions can be computed using the configurations generated by the reference potential according to Eq. (2). The above procedure was employed by Elliott and coworkers81–84 to develop a transferable force field at the atomistic level. They implemented a modified SW TPT to fused HS chains as the reference term and a step function

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.220.71.24 On: Mon, 11 Aug 2014 05:20:30

234104-5

A. F. Ghobadi and J. R. Elliott

with four wells as the perturbative potential. The athermal hard chains were simulated at a single temperature and several densities covering the entire density range using Discontinuous Molecular Dynamics85 (DMD). First and second order perturbation contributions were fitted to the exact simulation data using simple polynomials of density. The perturbation potential was further optimized to improve the predictions of vapor pressure and saturated liquid density. The methodology, named as DMD/TPT, was successfully applied to more than 500 compounds to predict VLE data and critical properties86 of pure fluids as well as mixtures of associating and nonassociating molecules. But, in contrast to the developed force field, the resulting EOS, named as SPEAD, was formulated with no explicit form of RDF in the definition of the A1 term. Therefore, its application to inhomogeneous systems would suffer from the same problems as many other SAFT EOSs. The simulation/perturbation hybrid procedure can be used to systematically examine the accuracy of perturbation contributions in SAFT models.86 Due to lumping of all Helmholtz energy contributions together and fitting directly to experimental data, the consistency of individual perturbation terms with molecular simulation is not always guaranteed. Any inconsistency between the terms in the EOS and their atomistic equivalents undermines our ability to achieve a consistent and verifiable theory for both homogeneous and inhomogeneous fluids. For example, the Macroscopic Compressibility Approximation (MCA) and LCA are common approximations in many SAFT implementations to estimate the A2 term. But these approximations have been demonstrated to provide poor accuracy, even qualitatively, when compared to molecular simulation.87, 88 Our goal is to minimize these kinds of discrepancies, even as we extend to higher order TPT and inhomogeneous fluids. The consistency between SAFT methodology and molecular simulation was improved, to some extent, by Avendano et al.21 where a coarse-grained version of SAFT-γ EOS was introduced for the single-site representation of CO2 , SF6 , and CF4 interacting with a Mie potential. In the SAFT-γ Mie EOS, an empirical correction factor is added to the LCA estimate of the A2 contribution in order to reproduce the fluctuation term of the Mie potential. Also, a third order term

J. Chem. Phys. 139, 234104 (2013)

was introduced in terms of an empirical exponential function. The exponents of the Mie potential (λr , λa ), which represent the steepness and range of attractive interactions, along with σ and ε were fitted to experimental critical temperature, surface tension, and VLE data. The optimized parameters were further tested in molecular simulation to demonstrate the level of agreement between calculated and simulated VLE and second derivative properties such as heat capacity and isothermal compressibility. The satisfactory agreement achieved from this comparison was the motivation to introduce the SAFT-γ Mie EOS as a tool for force-field optimization at a coarse-grained level.21 A meaningful bridge between theory and molecular simulation was also established to improve the global description of thermodynamic properties for chains interacting with the Mie potential,52 to speed up VLE calculations89 and force field optimization90, 91 with other EOSs such as PC-SAFT. Figure 2 compares SAFT-γ Mie and PC-SAFT perturbation contributions with simulation results of spheres in the framework of BH TPT.92 We chose PC-SAFT as a prototype of commonly used SAFT models. The simulation details are presented in Sec. III. It should be noted that we were able to reproduce, within statistical uncertainty, the simulation results of Smith et al.87 with both Monte Carlo (MC) and Molecular Dynamic (MD) simulations at T∗ = TkB /ε = 1. However, the A2 term of Cardenas et al.93 was not consistent with any of other simulations, illustrating the care that is required to maintain consistency at all length scales. Figure 2(a) illustrates satisfactory agreement for the A1 term as obtained by the SAFT-γ Mie EOS, while PC-SAFT underestimates the first order perturbation contribution especially on the liquid side. The overall agreement between simulation and SAFT-γ theory for the A2 term is acceptable, as shown in Figure 2(b), although there are deviations at high densities. The A2 term of PC-SAFT, however, mischaracterizes the qualitative trend of the simulation. As pointed out earlier, this is directly related to application of the MCA. The A3 term elevates questions of consistency to a higher order. As illustrated in Figure 2(c), SAFT-γ Mie only agrees with the simulation data in the gas phase limit. A3 in the SAFT-γ Mie EOS is estimated by fitting to experimental data

FIG. 2. First (a), second (b), and third (c) order perturbation contributions of monomer at T∗ = 1 from BH TPT. MC simulation of Smith et al.87 (squares), MC simulation of Cardenas et al.93 (diamonds), MC simulation of this work (circles), MD simulation of this work (triangles), SAFT-γ Mie EOS (solid line), and PC-SAFT EOS (dashed line). If not shown, error bars are smaller than the symbol size.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.220.71.24 On: Mon, 11 Aug 2014 05:20:30

234104-6

A. F. Ghobadi and J. R. Elliott

and, as a result, represents an effective contribution that may incorporate higher order perturbation terms. So, it is referred to as A3 but, in fact, the correspondence between simulation and theory is uncertain. At moderate densities, the SAFT-γ Mie prediction underestimates simulation data and more importantly, it fails to represent a crucial feature of A3 , which is being “positive” for the dense liquid. A positive A3 was observed from both MC and MD simulations when applying the BH split of the potential model. The positive A3 term may seem surprising in that it indicates a positive contribution to the EOS from what should be an attractive term. Nevertheless, it must be considered in light of the BH split and the A2 term. Note that the magnitude of the A2 term in the BH split increases monotonically with density. If the liquid structure is dominated by repulsive forces, as asserted by perturbation theory of liquids, then the first order term should become more prominent and the higher order terms (A2 and A3 ) should diminish. Therefore, the A2 and A3 terms must be considered in conjunction when a BH split is applied, such that the positive A3 and increasingly negative A2 cancel each other. We interpret this observation to indicate less favorable convergence of BH TPT. It should be noted that our simulation data for perturbation contributions of BH TPT is in excellent agreement with the recently published work of Lafitte et al.52 over the density range for which they display results. They also showed that the qualitative trend of perturbation terms (Figure 2) holds for various repulsive exponents of the Mie potential. Recently, the SAFT-γ Mie methodology was extended to chains of tangent soft spheres with 3 and 6 freely jointed segments as a coarse-grained representation of n-C10 H22 and n-C20 H42 ,22 respectively. It was shown that a third order BH TPT adapted to the Mie potential with adjustable exponents could provide accurate VLE, surface tension, and derivative properties. It is of interest to evaluate the consistency between perturbation contributions of the SAFT-γ Mie EOS for chains with those of simulation. To have a meaningful comparison, we performed the simulation for a three site tangent-sphere chain as a coarse-grained description for nC10 H22 .22 Throughout this work, since MD and MC simulations produced the same results for perturbation contri-

J. Chem. Phys. 139, 234104 (2013)

butions of monomers, MD simulation was carried out for chains to facilitate efficient application of parallel processing. Figure 3 shows that the level of agreement between SAFT-γ Mie and molecular simulation for perturbation contributions is qualitative, but not quantitative.92 Comparing perturbation contributions of monomers and chains, as depicted in Figures 2 and 3, demonstrates that both the magnitude and shape of perturbation terms change with number of segments. Changes in the shape indicate that the RDF of a three-segment chain fluid is different from that of three non-bonded monomers at the same density. Therefore, this result raises a question about the underlying assumption in many SAFT models that Amonomer Ach Ares =m + . N kB T N kB T N kB T

(9)

It seems more reliable if, rather than monomers, chains be recognized as the reference system in the SAFT approach. This concept is implemented in the PC-SAFT EOS, for example, and is believed to be the underlying reason for the better agreement of the A1 term of PC-SAFT with simulation data (Figure 3(a)). It should be noted that while perturbation contributions of PC-SAFT are obtained in the BH TPT framework, the effective hard sphere diameter is formulated for a shoulder potential. This inconsistency causes ambiguity when repulsive and attractive contributions of the PC-SAFT EOS are compared to molecular simulation. Extending the analysis to even higher order, Pavlyukhin has studied perturbation coefficients for SW TPT27 to fourth order in temperature. His results, which are consistent to second order with those of Elliott and coworkers,94, 95 suggest that the perturbation series converges slowly enough that truncation of the series must be carefully considered. The results for A2 of SW potentials exhibit relatively large values even at high density. On the other hand, the A2 of a step potential94 exhibits a steep decay on the liquid side, even when the reference system is identical to that of SW TPT. This means the discrepancy cannot be attributed to the reference potential, so it must be an impact of the SW perturbation potential. As a possible explanation, we note that the SW potential is flat across its entire range, which excessively accentuates the role of fluctuations near the outer limit of the range. Step potentials

FIG. 3. First (a), second (b), and third (c) order perturbation contributions of 3-segment chain obtained from MD simulation at NVT ensemble with T = 300 K and N = 400 using BH TPT. MD simulation of this work (triangles), SAFT-γ Mie EOS (solid line), and PC-SAFT EOS (dashed line). If not shown, error bars are smaller than the symbol size.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.220.71.24 On: Mon, 11 Aug 2014 05:20:30

234104-7

A. F. Ghobadi and J. R. Elliott

approach zero at long range (similar to the LJ model) and this seems to cause the A2 and A3 terms to approach zero more rapidly at high densities. These observations are consistent with those of Zhou and Solana88 suggesting that the convergence of perturbation theories depends somewhat on the nature of both the reference and perturbation potentials. The trends of perturbation contributions presented in this work for BH TPT also suggest that the perturbation series converges slowly. In this case, the reason originates in the definition of the reference potential. In BH TPT, the attractive interactions are improved but the reference potential is truncated so that it does not include the entire repulsive portion of the potential (refer to Figure 1). On the other hand, the WCA splitting of the potential leads to a complete separation between attractive and repulsive interactions. Because more information is encompassed by the reference term of WCA TPT, the corresponding perturbation series converges faster.28 The faster convergence with the WCA split is illustrated for spheres in Figure 4.92 In WCA TPT, second and higher order perturbation terms tend to zero at high density, consistent with the perspective of perturbation theory.79 The drawback of WCA TPT is its indirect adaptation to the existing theory of HS fluids. This caveat is addressed in Sec. IV A. In summary, the preceding review of SAFT approaches reveals that, when there is a robust connection between simulation and theory, as in the case of the iSAFT EOS, the accuracy at the bulk level is limited. On the other hand, when the bulk fluid properties are well-predicted, as in the case of SAFT-γ and PC-SAFT EOSs, a detailed correspondence, at the atomistic level, between molecular simulation and the EOS is not achieved. Note that we distinguish here between the atomistic level (e.g., a three site model for CO2 ) and a coarse grained model (e.g., a one site model for CO2 ). In this work, we propose a methodology to fill this gap between theory and molecular simulation. To do so, a site-based SAFT-γ EOS for fused heteronuclear chains interacting with the LJ potential is developed. A rapidly convergent WCA TPT at third order accounts for attractive interactions. Perturbation contributions (A1 , A2 , and A3 ) are fitted to the fluctuations of the MD simulations for the WCA reference system by introducing an explicit form for the effective RDF of chains. Consequently, a direct correspondence between parameters of the EOS and united-atom description of real molecules for LJ-based force fields is established. The resulting site-based

J. Chem. Phys. 139, 234104 (2013)

EOS provides a convenient platform for both bulk and inhomogeneous fluids near and far from the critical region. III. SIMULATION DETAILS

MC simulations of monomers were carried out in the NVT ensemble for the repulsive reference potential of BH theory (Figure 1(b)) using an in-house developed computer code, which is validated in Ref. 96. In each run, 500 particles were placed in a cubic simulation box with standard Periodic Boundary Conditions (PBC) and an f.c.c. lattice for initial configuration. In BH TPT, the cut-off distance for reference potential simulation is equal to σ and no long-range corrections should be applied during the simulation. In each cycle, one particle was chosen at random and a linear displacement of the particle was implemented using the Metropolis algorithm with an acceptance ratio of 50%. After the initial equilibration period of 5 × 107 cycles, configurations were stored every 5000 cycles. After the simulation, these configurations were employed in conjunction with the perturbation potential shown in Figure 1(b) and a cut-off distance equal to 14 Å to calculate perturbation contributions. This cut-off distance is chosen to be consistent with the spherical cut-off distance of common transferable force fields such as Transferable Potential for Phase Equilibria (TraPPE)97 and Nath-Escobedo-de Pablo revised force field (NERD)98 for bulk fluid. During post processing, standard long-range corrections were applied to A1 while the corresponding correction to higher order terms was practically zero. It should be noted that the number of particles in each run must be large enough to ensure that the cut-off distance is smaller than half of the simulation box. As pointed out by Pavlyukhin,27 A1 and A2 terms reach their equilibrium values very quickly while a large number of configurations are needed to achieve reliable statistics for A3 and higher order perturbation contributions. In this work, 4 × 108 –1 × 109 cycles were used for the production period depending on the system. For WCA TPT, the MC simulation details were the same but the cut-off distance of the reference potential was equal to σ × 21/6 . MD simulations of monomers as well as chains were carried out with LAMMPS99 molecular dynamic package. The Velocity-Verlet time integration method, a time step of 1 fs, a cubic simulation box, and standard PBC were used in all simulations. Initial configurations were prepared by

FIG. 4. First (a), second (b), and third (c) order perturbation contributions of monomer from BH and WCA perturbation theories obtained by MC simulation at NVT ensemble with N = 500 and T∗ = 1. If not shown, error bars are smaller than the symbol size.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.220.71.24 On: Mon, 11 Aug 2014 05:20:30

234104-8

A. F. Ghobadi and J. R. Elliott

J. Chem. Phys. 139, 234104 (2013)

PackMol.100 For BH TPT, simulations were performed in the NVT ensemble with Nose-Hoover101 thermostating and a damping coefficient equal to 200 fs. The initial configuration was equilibrated for 10 ns and system configurations were stored to a file every 1 ps for 80 ns to make sure reliable statistics are obtained for the A3 term. For WCA TPT, the initial configuration was melted by a 2 ns run in the NVT ensemble followed by an 8 ns run in NVE ensemble. The Nose-Hoover101 method was also used for WCA TPT runs. For the production period, NVE simulations were performed for 40–80 ns depending on the system and density of interest. It should be noted that the above simulation times are required to obtain a statistically reliable A3 term. A1 and A2 contributions require about 5 ns to reach their equilibrium values. Using the NVE ensemble for the production period is to ensure that the thermostating method does not affect evaluation of the perturbation contributions. Typically, in WCA TPT perturbation contributions converge faster comparing to BH TPT and fewer configurations are required in the computation of perturbation contributions. Molecular topology and force field parameters are based on TraPPE.97 In the TraPPE force field, bond length is kept constant while in the LAMMPS MD package, bond length cannot be fixed for molecules with more than 3 united atom sites. It has been shown102 that TraPPE results can be reproduced if (Kbond ≈ 1.5 × Kangle ) where Kbond and Kangle are bond stretching and angle bending force constants, respectively. As pointed out in Refs. 27 and 93, the number of particles does not affect computed perturbation contributions if enough configurations are stored. In this work, MD simulation of monomers was carried out with N = 500 and for chains there were enough molecules in the simulation box to have at least 1000 united atom sites. Details of post-processing to calculate perturbation contributions are the same as described above for MC simulation. Samples of the scripts for LAMMPS simulation and the C++ code for post processing to obtain A1 , A2 , and A3 are given in the supplementary material.92 The simulation data for compressibility factor and Helmholtz free energy of all systems studied in this work are provided as the supplementary material.

a = a id + a0 + a pert + a polar ,

To formulate the SAFT-γ EOS for LJ chains consistent with UA force fields such as TraPPE,97 NERD,98 and Optimized Potentials for Liquid Simulations-United Atom (OPLS-UA),103 heteronuclear soft chains formed with fused united-atom sites interacting with the WCA potential are considered as the reference system: ⎧ 12 6 σ σ 1 ⎨ 4ε + ε r < (2 /6 )σ − r r . (10) u0 (r) = ⎩ 1 0 r ≥ (2 /6 )σ

(11)

(12)

where a = A/NkB T, superscript id and polar represent ideal chain and partial charge contributions to Helmholtz energy, respectively. The polar term will be the subject of succeeding papers and is set to zero in this work. The Helmholtz energy of the reference system is introduced as the sum of WCA spheres and chain contributions: a0 = a WCA + a ch .

(13)

Attractive interactions are formulated by a third-order perturbation theory as a pert =

A1 A3 A2 + 2 + 3, T T T

(14)

with Ai defined in Eq. (2). In Sec. IV A and IV B, the functional form and underlying theory of the reference and perturbation contributions are discussed.

A. WCA chain reference system

To construct a chain EOS in the SAFT perspective, one needs an accurate representation of the monomer contribution. In the perturbation theory of soft potentials, properties are described in terms of an effective HS diameter (dEHS ). Among all the proposed equations, those that consider a temperature and density dependent dEHS seem to be more precise in describing WCA fluids over the broader temperature and pressure range.104 On the other hand, the derivative terms for pressure and compressibility become extremely complicated, even for pure fluids. Rather than introducing a complicated dEHS (T, ρ), Heyes and Okumura105 proposed that a slight modification in the Carnahan-Starling34 (CS) EOS for HS can avoid the necessity of a density dependent dEHS . They showed that a generalized CS EOS with Boltzmann’s106 suggestion for effective diameter was accurate for WCA spheres: dE 21/6 . = √ σ (1 + T ∗ )1/6

IV. CHARACTERIZING REFERENCE AND PERTURBATION CONTRIBUTIONS

The perturbative potential in WCA TPT reads as ⎧ 1 ⎪ r < (2 /6 )σ ⎨ −ε . u (r) = σ 12 σ 6 1 ⎪ ⎩ 4ε r ≥ (2 /6 )σ − r r

In our SAFT-γ WCA method, the total Helmholtz energy of the homogeneous fluid is

(15)

Equation (15) can reproduce simulation data for structural properties and compressibility factor of WCA spheres for a wide temperature and density range. They avoid naming dE as effective HS diameter since it is not used with a HS EOS. In this work, we drop the superscript E for simplicity. Nevertheless, the preliminary results of this work show that the application of Eq. (15) to form chain molecules in the SAFT approach leads to underestimation of second virial coefficients at high temperatures and overestimation of saturated liquid density for long chains. On the other hand, the implementation of d proposed by Elliott and Daubert107 for the WCA potential provides better agreement with simulation data and virial coefficients of WCA fluids: d 1 = , (16) 1 ∗ ∗ ∗ 2 rmin (n (b1 n − b2 )(T ) − (n∗ − 1)T ∗ + 1) 2n∗ +1

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.220.71.24 On: Mon, 11 Aug 2014 05:20:30

234104-9

A. F. Ghobadi and J. R. Elliott

rmin = σ

n∗ 6

J. Chem. Phys. 139, 234104 (2013)

1/(n∗ −6) (17)

.

In Eq. (16), b1 = 0.0093 and b2 = 0.0592 and T∗ = TkB /ε and ε is taken directly from force field of interest. n∗ controls the steepness of the interaction potential and for the Lennard-Jones fluid is equal to 12. In this work, n∗ is considered as an adjustable parameter for each united-atom site type. This can be rationalized by the observation that different site types present different degrees of softness owing to repulsive overlaps of fused spheres. It also facilitates mapping the properties of soft fused-chains with a theory that is formulated for tangent hard spheres. In the same fashion as Heyes and Okumura,105 we introduce a Generalized Boublik-MansooriCarnahan-Starling-Leland108, 109 (G-MCSL) EOS for WCA sphere mixtures as Z WCA =

n3 n1 n2 γ1 + 1 − n3 n0 (1 − n3 )2 +

n32 γ2 n3 n32 γ3 − , (18) 3 12π n0 (1 − n3 ) 36π n0 (1 − n3 )3

(nω , ω = 0−3) are Scaled Particle Theory (SPT) variables: nω = ρ

NC 

xi M(ω) i ,

M(0) i =

NS 

a WCA = f3 =



NS 

2 νk,i sk dkk ,

a ch,SAFT−γ = −

i=1

xi

 NS 

(19)

NS π 3 = νk,i sk dkk . 6 k

M(3) i

k

In Eq. (19), ρ is chain number density, NC is number of compounds in the mixture, x is mole fraction, d is given by Eq. (16), NS is number of united-atom site types forming a chain, ν k,i is number of site type k in the chain of type i, and s is the shape factor that is introduced in SAFT-γ EOS.14 It should be noted that n3 is equal to packing fraction (η) for bulk fluid and s is assumed to be temperature independent. For homonuclear chains, G-MCSL EOS with (γw = 1, w = 1−3) will reduce to CS EOS: Z HS =

4η − 2η2 . (1 − η)3

(20)

The Helmholtz energy corresponding to the G-MCSL EOS is obtained by following the procedure of Roth et al.:30

n1 n2 γ1 + f3 n32 , (1 − n3 )   (3γ2 − 2γ3 ) n3 − 1.5 (γ2 − γ3 ) 2n3 − n23 + γ3 (1 − n3 )2 ln (1 − n3 )

1 ρ



−n0 ln (1 − n3 ) +

36π n23 (1 − n3 )2

In the SAFT-γ approach, the Helmholtz energy due to formation of heteronuclear fused chains (ach ) is given by NC 

M(1) i =

k

M(2) i

1 νk,i sk dkk , 2 k NS

νk,i sk ,

WCA

is the compressibility factor of WCA spheres and where Z (γw , w = 1−3) are adjustable constants to fit simulation data of WCA spheres. The fitting procedure is described in Sec. V.

ω = 0 − 3,

i

(22)

k=1

where the RDF at contact (gii ) is a function of effective molecular parameters.14 In the framework of Wertheim’s TPT1, the above equation only holds if all sites have the same contribution to chain formation (homonuclear). Furthermore, the term in parentheses can be easily zero for some values of shape factor. For example, if for a two-site molecule (nitrogen, ethane, etc.) shape factor happens to be 0.5, ach would be zero. When each site type has a different contribution to the chain formation, Wertheim’s TPT1 for a heteronuclear chain would be ⎛ ⎞ {k  } NC NS    −1 ln gkk ⎠ , xi ⎝ sk (23) a ch = 2  k i=1 k=1 where {k } denotes the set of bonded partners of site k. After adapting an effective set of molecular parameters similar

.

to those introduced in the original SAFT-γ EOS,14 we can rewrite Eq. (23) as

 νk,i sk − 1 ln gii ,

(21)

a ch = −

NC 

 xi

i=1

 NS 1 νk,i sk N Bk ln giiWCA , 2 k=1

(24)

where NBk stands for number of bonding sites for site-type k and giiWCA is the RDF at contact for WCA fluid. The chain contribution defined in Eq. (24) is valid for any arbitrary set of site parameters and is only zero when NBk is zero for all k. The expression of giiWCA that corresponds to the G-MCSL EOS is   gijWCA d¯ij ; n2 ; n3 =

(γ1 + γ2 ) n2 1 + 1 − n3 4 (1 − n3 )2 +

(3γ2 − γ3 ) n22 36 (1 − n3 )3

¯ ¯ dii djj ¯ dii + d¯jj

¯ ¯ 2 dii djj , (25) ¯ dii + d¯jj

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.220.71.24 On: Mon, 11 Aug 2014 05:20:30

234104-10

A. F. Ghobadi and J. R. Elliott

J. Chem. Phys. 139, 234104 (2013)

with the effective size parameter of component i (d¯ii ) defined as NS 

d¯ii3 =

3 νk,i sk dkk

k=1 NS 

(26)

. νk,i sk

k=1

From Eq. (25) with (γw = 1, w = 1−3), the RDF at contact of the CS EOS for homonuclear chains is obtained 1 − η/2 g HS (σ ; η) = . (27) (1 − η)3 In summary, the a0 term of our SAFT-γ WCA EOS, has three constants (γw , w = 1−3) and three adjustable parameters per site type (n∗ , σ , s). It should be noted that for a third order perturbation theory, a second order expression for RDF, rather than that of Eq. (25), ought to be used. In this work, since the entire reference term is eventually fitted to the exact simulation data of WCA chains, using a more complicated functional form for RDF only leads to a different set of adjustable parameters without any noticeable change in the ability of Eqs. (18)–(26) to reproduce simulation results. Furthermore, the zeroth order approximation for RDF provides

WCA g0,kl

a clear separation between repulsive and attractive interactions which facilitate the procedure of parameter evaluation for SAFT-γ WCA EOS as discussed in Sec. V. B. Perturbation contributions

As described in Sec. II A, in SAFT-γ and SAFT-VR methodology, MVT is used to estimate the first order perturbation contribution. In this work, A1 resembles that of the SAFT-γ EOS14 but it is reformulated for WCA TPT: A1 = −ρ

NC  NC  i=1 j =1

NS NS  

  WCA1 WCA νk,i νl,j sk sl αkl g0,kl neff1 3,kl .

k=1 l=1

(28) However, one should note that in SAFT-γ and SAFTγ Mie equations, perturbation terms are evaluated for monomers and then are scaled by the number of segments in the chain. In this work, similar to the PC-SAFT approach, perturbation terms are directly evaluated for chains to be consistent with the definition of the reference term. The underlying reasons behind this choice of the reference system are WCA eff1 (n3,kl ) is RDF at condiscussed in Sec. II B. In Eq. (28), g0,kl tact for homonuclear chains evaluated at an effective density:

 2  2   eff1 2    eff1  4 1 − neff1 + (3γ2 − γ3 ) neff1 + 3 γ1 + γ2 neff1 3,kl 3,kl − n3,kl 3,kl n3,kl = , 3  4 1 − neff1 3,kl neff1 3,kl =

eff1 neff1 3,k / (σk εk ) + n3,l / (σl εl )

1 / (σk εk ) + 1 / (σl εl )

11 12 13 2 neff1 3,w  = aw  + aw  n3 + aw  n3 .

In Eq. (31), (aw1w , w = 1−3) are adjustable parameters of site type w  for A1 term. It should be noted that as opposed to SAFT-γ EOS, the effective packing fraction is not constrained to be zero in the ideal gas limit. Although, aw11 appears to be zero for monomers, only a non-zero value fits the simulation WCA1 is the interaction data of chain molecules. In Eq. (28), αkl volume for the perturbative potential defined in (11) and is given by WCA1 αkl

xi xj

 2 16 σkl 1 =− −εkl 4π r 2 dr 2kB λ(T )  ∞ σkl 12 σkl 6 2 4π r dr . (32) + 1 4εkl − r r 2 6 σkl

Note that the lower domain of the first integral of the right hand side is not zero even though the perturbative potential is defined to be non-zero from 0 to ∞. The reason is that for WCA fluids, the RDF vanishes for radial distances less than λ(T). Simulation data for the RDF show that when T → 0,

,

(29)

(30) (31)

λ(T) → 21/6 σ and as T → ∞, λ(T) → 0 which resembles the trend of d versus T for WCA fluid. Unfortunately, setting λ(T) = d(T) causes |A1 | to increase with increasing temperature, which contradicts the trend that is observed for |A1 | from simulation. The inconsistency is probably due to using the MVT to approximate Eq. (6) in the form of Eq. (8). In this work, λ(T) was set to σ kl for all temperatures, leading to the following expression for interaction volume of WCA fluid: 2 16 √ εkl 3 WCA1 − αkl =−π 2 σ . (33) 3 9 kB kl Finally, the corresponding combining rules for ε and σ completes the functional form of the first order term: √ √ (34) εkl = εk εl , σkl = σk σl . Note that the combining rules for A1 (Eqs. (30) and (34)) differs from those of original SAFT-γ EOS. This set of combining rules eliminates the necessity of cross group energy parameters, ensures the transferability of adjustable parameters

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.220.71.24 On: Mon, 11 Aug 2014 05:20:30

234104-11

A. F. Ghobadi and J. R. Elliott

J. Chem. Phys. 139, 234104 (2013)

defined for each united atom site type, and leads to superior results for mixtures without introducing any binary interaction parameter. For the A2 term, as it was discussed earlier, neither LDA nor MCA are capable of reproducing the fluctuation term of WCA chains. In this work, we apply MVT, as implemented for A1 , combined with MCA, to obtain NC NC   xi xj A2 = −ρK WCA i=1 j =1

×

NS  NS 

  WCA2 WCA eff2 νk,i νl,j sk sl αkl g0,kl n3,kl ,

(35)

k=1 l=1 WCA eff2 where g0,kl (n3,kl ) is given by Eq. (29) with an effective packing fraction (neff2 3,kl ) that is optimized to reproduce the A2 computed by simulation. The combining rules for neff2 3,kl resemble those of A1 and are not repeated here. Note that the effecWCA eff2 (n3,kl ) introtive density that is used in definition of g0,kl duces three more adjustable parameters per site type (aw2w , w = 1−3) for the A2 term. The isothermal compressibility (KWCA ) corresponding to G-MCSL EOS (Eq. (18)) reads as (1−n3 )4  , K WCA = 1+(6γ1 −2)n3 +(1+9γ2 −6γ1 )n23 +γ3 n43 −4n33

(36) which reduces to KHS of CS EOS if (γw = 1, w = 1−3). The interaction energy for the second order term is computed by  2 16 σkl 1 WCA2 2 = 2 εkl 4π r 2 dr αkl 4kB λ(T )  ∞ σkl 12 σkl 6 2 2 2 − 4π r dr . + 1 16εkl r r 2 6 σkl (37) Setting λ(T) = σ kl results in  √  4 2−1 εkl 2 3 WCA2 =π σkl . αkl 3 kB

(38)

WCA1 WCA2 Note that both αkl and αkl are positive numbers. For the A3 term, the simulation data of WCA chains of length m show that ε A3 = m A2 (n3 ) , (39) kB

which is simply proportional to the A2 term multiplied by a density-dependent Gaussian-type distribution . After applying the concept of shape factor and appropriate mixing rules, A3 can be reformulated as follows: NC NC   A3 = −ρK WCA xi xj

where (aw3w , w = 1−3) are adjustable parameters of site type w  for the A3 term. The corresponding equations to calculate chemical potential are presented in the Appendix. V. PARAMETER EVALUATION PROCEDURE

In summary, the SAFT-γ WCA EOS has 12 adjustable parameters to be evaluated for each united-atom site type. The parameters of the repulsive reference system and those of the perturbation terms are fitted independently to the corresponding terms computed from molecular simulation of the WCA reference system. In this work, these parameters are evaluated for monomers as well as CH3 and CH2 united-atom site types using the TraPPE97 force field for n-alkanes. Extension to other site types will be pursued in the future. The simulations were carried out at 3 temperatures and 16 densities to cover the entire phase diagram from ideal gas to dense liquid. For each compound, the three temperatures were chosen to cover the subcritical, supercritical, and critical point region of the phase diagram. A. Monomers

For monomers, σ and ε are taken directly from the force field of interest, s = 1.0, n∗ = 12, and NB = 0. Therefore, the chain term is zero and the reference term will have only three adjustable constants, (γw , w = 1−3), to be fitted with simulation data. These constants are fitted to the compressibility factor of the WCA fluid obtained from MD simulations at 3 different reduced temperatures (T ∗ = 0.67, 1.35, and 2.0) and 16 densities. The resulting values are γ 1 = 0.94, γ 2 = 1.45, and γ 3 = 4.5 and are site type independent. Each simulation run of the WCA reference system in the NVT ensemble results in a set of perturbation contributions (Ai , i = 1−3) for that particular temperature and density. The simulation run close to the critical temperature is carried out for ∼40 ns to achieve reliable statistics for A3 term. To save computational time and since A3 is proportional to A2 , we assume that the temperature dependency of A3 can be described if the A2 term is well-parameterized for a wide temperature range. Therefore, A3 is only computed at a single temperature. The simulations for the other two temperatures are performed for relatively shorter times (∼5 ns) because A1 and A2 reach their equilibrium values very quickly. Then, adjustable parameters of perturbation contributions (aij , i, j = 1 − 3) are regressed to reproduce the corresponding perturbation terms (Ai , i = 1−3) from simulation. Figure 5 shows that perturbation parameters that are obtained from simulations with TraPPE for CH4 (as a monomer) are transferable to a wide range of force field parameters.92

i=1 j =1

×

NS  NS 

  2 WCA2 WCA νk,i νl,j sk2 sl αkl g0,kl neff2 3,kl kl ,

(40) (41)

For site types forming chains, ε is taken directly from that of force field. However, the preliminary results of this work show that a satisfactory fit is only obtained if the following relation is assumed:

(42)

σ = σ σ F F ,

k=1 l=1

k / (σk εk ) + l / (σl εl )

kl = , 1 / (σk εk ) + 1 / (σl εl )  2   εkl 31 aw exp −aw32 n3 − aw33 ,

w  = kB

B. Chains

(43)

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.220.71.24 On: Mon, 11 Aug 2014 05:20:30

234104-12

A. F. Ghobadi and J. R. Elliott

J. Chem. Phys. 139, 234104 (2013)

FIG. 5. First (a), second (b), and third (c) order perturbation contributions of monomers with different values of ε and σ as given in the legend of (a). All lines are the results of SAFT-γ WCA EOS and symbols are simulation data.

where σ FF is the size parameter of an arbitrary force field (FF) and σ is a “shift factor” that relates the united-atom size parameter used in the SAFT-γ EOS (σ ) to that of the force field. The reason that σ is not equal to 1 has its origin in the SAFT approach where chains are formed with tangent HS segments. For fused chain molecules, the cavity correlation function is required for distances less than the size of segments.110 For those distances, the cavity correlation function is not known at the current state of theory. The concept of shape factor within the framework of the SAFT-γ EOS provides a convenient mathematical procedure to adapt Wertheim’s TPT1 for fused chains. But, it does not solve the problem from a physical perspective. Consequently, the parameter σ is defined to circumvent the deficiency of SAFT theory for fused chains. Such shift factors are also needed for tangentially formed chains22 to compensate for underlying simplifying assumptions in evaluating correlation functions.111 Note that in the SAFT-γ WCA EOS presented here, neither a shift factor nor any cross interaction parameter is needed for ε (i.e., ε = 1). For CH3 and CH2 , NB is equal to 1 and 2, respectively. The rest of the parameters for reference term, (s, n∗ , and σ ), are fitted to compressibility factor of molecules interacting with WCA potential. It should be noted that the best fit is obtained when CH3 of ethane is distinguished from that of other n-alkanes. So, adjustable parameters of CH3 in ethane are fitted to ethane and those of [CH3 and CH2 ] sites for the rest of alkanes are fitted simultaneously to simulation data of nbutane and n-dodecane. Similar to the case of monomers, the adjustable parameters of the perturbation terms are regressed to reproduce simulation data obtained from Eq. (2). The resulting parameters are summarized in Tables I and II for reference and perturbation contributions, respectively. The large shift factor for CH2 should be understood in conjunction with its small shape factor. Based on the parameters of Table I, the effective volume of CH2 is 54% of a complete sphere while the corresponding number for CH3 is 77%. Note that CH2 is connected to two other segments while CH3 is an end segment. Because of the segment overlap in the chain backbone, it is expected to observe that the CH3 site behaves more like a complete sphere than CH2 . For CH2 , the relatively large value of n∗ = 24 reflects that it has one less H atom than

CH3 and fused-sphere repulsive overlaps on both sides. So, it is effectively “harder” than the CH3 site type. The agreement between the SAFT-γ WCA EOS and simulation data for the WCA reference compressibility factor and perturbation contributions of n-butane and n-dodecane is illustrated in Figures 6 and 7.92 As seen from Figure 7(a), A1 of n-dodecane goes to zero at the ideal gas limit. The ideal chain contribution (intramolecular interactions) for molecules with more than 4 united-atom sites is subtracted from simulation data to ensure the limiting behavior as density tends to zero. Intramolecular interactions are density independent and, therefore, have no impact on VLE calculations. C. Extension to force fields beyond TraPPE

Ideally, a generalized EOS would take non-bonded parameters (ε and σ ) of any force field and reproduce the VLE simulation data for that particular model. For monomers, it is easy to show that the EOS is generalized in this sense. In this section, we present a simple analysis to investigate the correlation, if there is any, between the adjustable parameters of the EOS and those of the force field. Unfortunately, our preliminary results show that generalization is not always possible for chains. To begin, assume that for any arbitrary adjustable parameter P: σ FF εF F , (44) , P = P T raP P E f σ T raP P E εT raP P E where PTraPPE indicates the adjustable parameter that is obtained by fitting to the TraPPE force field, σ TraPPE and εTraPPE are non-bonded parameters of the TraPPE potential, σ FF and TABLE I. Parameters of reference term for SAFT-γ WCA EOS fitted to TraPPE force field. Site type Monomer CH3 (ethane) CH3 CH2



NB

s

n∗

1.0 0.96674 1.01066 1.37077

1 1 1 2

1.0 0.820 0.748 0.210

12 12 12 24

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.220.71.24 On: Mon, 11 Aug 2014 05:20:30

234104-13

A. F. Ghobadi and J. R. Elliott

J. Chem. Phys. 139, 234104 (2013)

TABLE II. Parameters of perturbation contributions for SAFT-γ WCA EOS fitted to TraPPE force field. a11

a12

a13

a21

a22

a23

a31

a32

a33

0.0 0.10 0.142502 0.332588

0.340563 0.498678 0.354847 0.754575

− 0.26100 − 0.39002 − 0.23486 − 0.53453

− 0.45735 2.76712 2.62906 1.27759

0.33188 − 3.78101 − 3.75050 − 0.61385

− 2.42928 19.32828 15.02178 1.426387

0.4 0.7 1.6 4.5

30.0 27.0 3.00 1.95

0.15 0.18 0.12 − 0.50

Site type Monomer CH3 (ethane) CH3 CH2

εFF are non-bonded parameters of an arbitrary force field and f is a correlation function such that f(1,1) = 1. Note that the correlation function is equal to 1 for all adjustable parameters of monomers. We can define dimensionless variables and write P = P T raP P E f (σˆ , εˆ ). σˆ =

σ FF

, T raP P E

εˆ =

εF F

(45)

. T raP P E

σ ε Assuming a separation of variables:

P = P T raP P E X (σˆ ) Y (ˆε) .

(46)

We also assume X and Y are linear functions so that ∂P ∂P , X (σˆ ) = σˆ + 1 − ∂ σˆ ∂ σˆ (47) ∂P ∂P Y (ˆε ) = εˆ + 1 − . ∂ εˆ ∂ εˆ This choice of intercept guarantees X = 1 and Y = 1 when P is not correlated to force field parameters. After putting Eq. (47) in Eq. (46) and rearranging: P = P T raP P E [(δσ (σ F F − σ TraPPE ) + 1) × (δε (εF F − εTraPPE ) + 1)],

(48)

∂P ∂P , δε = F F , F F ∂σ ∂ε To obtain derivatives of all adjustable parameters with respect to σ FF and εFF , several simulations were conducted for ethane, butane, and octane to calculate the corresponding δσ =

Jacobian matrix for CH3 and CH2 site types. Intra-molecular parameters like bond angles and bond length were fixed to those of TraPPE. The results showed that for σ and n∗ the correlation function, f, is equal to 1. Shape factor varies only slightly with the inverse of σ FF as s = s TraPPE [δs (1/σ F F − 1/σ TraPPE ) + 1], ∂P . δs =  ∂ 1 / σ FF

(49)

Because shape factor is roughly proportional to the surface area, volume, and mean curvature integral,112 of the molecule, its value would be affected by particle size as well as temperature for soft particles. Nevertheless, in this work, s is assumed to be temperature independent for simplicity. The adjustable parameters of first and second order perturbation terms were changing with both σ FF and εFF for linear molecules:  TraPPE  δσij (σ F F − σ TraPPE ) + 1 a ij = a ij   ij F F  i = 1 − 2 TraPPE × δε (ε − ε )+1 . (50) j =1−3 These equations hold as long as intra-molecular parameters do not change significantly from those of TraPPE FF. TraPPE are given in Tables I and II, respectively. sTraPPE and a ij Table III lists the parameters of Eqs. (49) and (50) for CH2 and CH3 site types. For the corresponding parameters of A3 the correlation function is equal to 1. It should be mentioned that the applicability of these generalizations to nonlinear molecules has not been tested yet. A major distinction between this work and most alternatives for EOS development is that, rather than fitting the density derivatives of Helmholtz energy to experimental vapor pressure and liquid density, each perturbation term (A1 , A2 , and A3 ) as well as the reference system (A0 ) is individually fitted to the corresponding simulation data in the framework of WCA TPT. The resulting Helmholtz energy, Eq. (12), is only a function of force field parameters (ε and σ ) and serves as a meaningful bridge between molecular simulation and a classical EOS. VI. RESULTS AND DISCUSSION

FIG. 6. Compressibility factor of the reference WCA system for butane (red circles and lines) and dodecane (black squares and lines) at three temperatures. Symbols are simulation data and lines are SAFT-γ WCA EOS results. Error bar on the simulation data is less than the symbol size.

The fitting procedure described in Sec. V guarantees the consistency of reference and perturbation contributions to Helmholtz energy with LJ-type force fields. Nevertheless, all industrial applications require satisfactory agreement for derivatives of Helmholtz energy. Therefore, bulk fluid properties are compared with full potential simulations to evaluate the transferability and viability

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.220.71.24 On: Mon, 11 Aug 2014 05:20:30

234104-14

A. F. Ghobadi and J. R. Elliott

J. Chem. Phys. 139, 234104 (2013)

FIG. 7. First (a), second (b), and third (c) order perturbation contributions of butane (red circles and lines) and dodecane (black squares and lines) at three difference temperatures. Symbols and lines are simulation and SAFT-γ WCA EOS results, respectively. If not shown, error bars are smaller than the symbol size. The grey and blue spheres in the inset of (c) are representatives of CH2 and CH3 united atom sites, respectively. Note the slightly smaller size of CH3 as presented in TraPPE force field. Perturbation contributions are divided by corresponding carbon numbers (CN) for each chain.

of the SAFT-γ WCA EOS for phase equilibria calculations. Figure 8 compares compressibility factor (Z) and second virial coefficient (B2 ) as obtained by SAFT-γ WCA EOS and molecular simulations with TraPPE for parameterization compounds: methane (C1), ethane (C2), butane (C4), and dodecane (C12).92 The agreement between Z from the SAFT-γ WCA EOS and simulation is very satisfactory as expected for WCA TPT.28 The B2 of methane matches that of the LJ fluid for a wide temperature range up to 4000 K, reflecting the accuracy of the effective diameter given in Eq. (16). Similarly, the same level of agreement is achieved for molecules with moderate chain length, as shown in Figure 8(b) for ethane and butane. For dodecane, B2 is underestimated for T > 1000 K where perturbation contributions are relatively small. Considering the degree of agreement for B2 of monomers, the defect, therefore, should be attributed to the chain term presented in Eq. (22). Our preliminary results demonstrate that neither assuming a temperature-dependent shape factor nor increasing the temperature range for evaluation of adjustable parameters improves the B2 predictions at high temperatures. For example, we simulated the reference term of dodecane at T = 2000 K and tuned adjustable parameters of the repulsive chain term (σ , s) only to this temperature. The resulting EOS would still underestimate B2 at high temperatures. Hence, the discrepancy probably stems from the simplifying assumptions in the approximations that are made for the reference fluid cavity correlation function of fused molecules. The viability of such approximations at high temperatures needs more investigation. Vapor-liquid coexistence curves and vapor pressure are given in Figure 9.92 Satisfactory agreement is obtained for

bulk fluid properties when the SAFT-γ WCA EOS and TraPPE FF are compared, even in the critical point vicinity. The ability of the SAFT-γ WCA EOS, as a classical EOS, to provide reasonable prediction of critical properties is due to inclusion of the third-order perturbation term. This was also shown for the SAFT-γ Mie EOS that is based on BH TPT3.21, 22, 116 The third order term, as given in Eq. (2), is a measure of fluctuations in the second order term and, approximately, takes into account the characteristic density fluctuations in the critical region. It should be noted that the results of the SAFT-γ WCA EOS are not expected to match experimental data as the reference and perturbation parameters are fitted to the TraPPE force field. This is more visible for compressed liquid density and vapor pressure at lower temperatures where deviation of the TraPPE model from experimental data is more noticeable. Compressed liquid density for parameterization compounds at several selected isochors is depicted in Figure 10.92 The curvature in the prediction of the SAFT-γ WCA EOS at higher temperatures is an index for supplementing chain molecules with softness. Similar comparisons were carried out for n-alkanes that were excluded from fitting procedure. Figure 11 shows the second virial coefficient and compressed liquid densities for selected validation compounds with different chain length.92 As observed for parameterization compounds, isochoric densities and B2 are overestimated at low temperatures. It is also shown that B2 of n-eicosane is considerably underestimated at moderate to high temperature similar to the case of n-dodecane as discussed earlier. Coexistence densities and vapor pressure for validation compounds are depicted in Figure 12.92 Despite acceptable agreement between the SAFT-γ WCA EOS and the TraPPE

TABLE III. Parameters of Eqs. (49) and (50) for extension of SAFT-γ WCA EOS beyond TraPPE force field. Site type

δ σ 11

δ σ 12

δ σ 13

δ σ 21

δ σ 22

δ σ 23

δs

δ ε 11

δ ε 12

δ ε 13

δ ε 21

δ ε 22

δ ε 23

CH3 (ethane) CH3 CH2

1.4 0.9 2.5

− 0.1 0.4 − 1.9

0.1 0.2 − 2.8

− 0.8 − 2.1 − 0.3

− 4.5 − 5.0 − 0.7

− 3.4 − 3.7 − 0.9

0.5 0.7 0.9

− 0.1 − 0.7 0.17

0.0 − 0.1 − 0.09

− 0.1 − 0. 1 − 0.2

0.1 − 1.0 0.0

0.1 − 2.3 0. 0

− 0.2 − 2.4 − 0.2

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.220.71.24 On: Mon, 11 Aug 2014 05:20:30

234104-15

A. F. Ghobadi and J. R. Elliott

J. Chem. Phys. 139, 234104 (2013)

FIG. 8. Compressibility factor (a) and second virial coefficient (b) of methane (C1), ethane (C2), butane (C4), and dodecane (C12). Filled symbols and solid lines are representatives for full potential TraPPE simulation data and SAFT-γ WCA results, respectively. Simulation data of B2 for C2, C4, and C12 are taken from Schultz and Kofke113 and that of methane is given by Barker et al.114 The simulation results of Z are provided in this work. If not shown, error bars are smaller than size of symbols.

FIG. 9. Vapor pressure (a) and orthobaric densities (b) of methane (C1), ethane (C2), butane (C4), and dodecane (C12). Red unfilled symbols are experimental data from NIST WebBook.115 Black filled symbols are full potential TraPPE simulation data.97 Black solid lines are SAFT-γ WCA results of this work and dashed lines are results of SAFT-γ WCA EOS truncated after the second order term.

FIG. 10. Compressed liquid density of (a) methane (C1, blue lines and symbols), butane (C4, green lines and symbols) and (b) ethane (C2, red lines and symbols), dodecane (C12, black lines and symbols). Unfilled symbols are experimental data from NIST WebBook115 and filled symbols are TraPPE simulation data from this work. Solid and dashed lines are SAFT-γ WCA EOS results.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.220.71.24 On: Mon, 11 Aug 2014 05:20:30

234104-16

A. F. Ghobadi and J. R. Elliott

J. Chem. Phys. 139, 234104 (2013)

FIG. 11. (a) Compressed liquid density of propane (C3, blue symbols and lines), pentane (C5, red symbols and lines), and octane (C8, black symbols and lines). Symbols are the same as Figure 10. (b) Second virial coefficient of propane, pentane, octane, decane (C10), and eicosane (C20). Symbols are the same as Figure 8(b).

force field for molecules with short to intermediate chain length, extension to long chains is inferior. It can be seen that for eicosane, vapor pressure is underestimated comparing to simulation data. Our hypothesis at this time is that poor predictions of the SAFT-γ WCA EOS for long chain molecules are mainly due to ignoring intramolecular interactions in the formulation of the EOS. In TraPPE, and almost all LJ-based force fields, intramolecular interactions do exist between united-atom sites that are four bonds apart. These interactions have a significant impact on molecular topology and thermophysical properties. Additionally, the deficiency can be attributed to the ability of transferable EOSs in extension to long chain molecules.14, 72, 75 The particularly pointed shape of the binodal for n-eicosane may also reflect effects of nonclassical scaling such that the A3 term is not as influential for long chains as for short chains. Table IV summarizes critical properties and normal boiling temperature of selected n-alkanes as predicted by the SAFT-γ WCA EOS and TraPPE force field. Considering the error bars on simulation data, the agreement between simulation and theory is outstanding for a classical EOS without

using any Renormalization Group (RG) corrections.29, 120, 121 Average Absolute Deviation (AAD%) of predicted vapor pressure and saturated liquid density is summarized in Table V. Saturated liquid density and vapor pressure only deviate by 0.8% and 4%, respectively, from those of the TraPPE force field for n-alkanes ranging from methane to dodecane. The SAFT-γ WCA EOS is not developed to only match the TraPPE force field for phase equilibria calculations. In fact, given ε and σ , it can reproduce the results of any force field with a 12-6 LJ type interaction potential as long as the intramolecular potential does not deviate substantially from that of TraPPE. In Figure 13, the vapor pressure and binodals of ethane, butane, and octane are compared for the TraPPE97 (black triangles), NERD98 (red squares), and OPLS-UA103 (blue circles) force fields.92 The corresponding SAFT-γ WCA predictions are shown with solid lines of similar colors. Vapor pressure of n-alkanes follows the sequence: TraPPE > NERD > OPLS-UA. Saturated liquid density far from the critical point obeys the sequence: OPLS-UA > TraPPE > NERD. It is also known that critical pressure and temperature are significantly overestimated in the OPLS-UA

FIG. 12. Vapor pressure (a) and orthobaric densities (b) of propane (C3), pentane (C5), octane (C8), decane (C10), and eicosane (C20). Symbols are the same as Figure 9. TraPPE simulation data of C10 and C20 are taken from Müller and Mejía102 and those of C8 are taken from Errington and Panagiotopoulos.117 The data in the last reference are read from plots using PlotDigitizer software.118 Experimental vapor pressure and liquid density of C20 are taken from DIPPR119 data base.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.220.71.24 On: Mon, 11 Aug 2014 05:20:30

234104-17

A. F. Ghobadi and J. R. Elliott

J. Chem. Phys. 139, 234104 (2013)

TABLE IV. Critical properties and normal boiling temperature (Tb ) from TraPPE force field and SAFT-γ WCA EOS. Subscripts show the statistical uncertainty of the final digit. TraPPE97, 117

SAFT-γ WCA EOS Compound

Tc (K)

Pc (MPa)

ρ c (g/cm3 )

Tb (K)

Tc (K)

Pc (MPa)

ρ c (g/cm3 )

Tb (K)

Methane Ethane Propane Butane Pentane Octane Dodecane

198.8 309.2 369.9 429.9 477.3 579.3 671.4

5.60 5.47 4.61 4.60 4.24 3.09 2.08

0.158 0.205 0.222 0.237 0.242 0.231 0.201

109.9 178.2 222.2 262.1 297.7 387.0 479.6

191.41 3022 3682 4234 4673 5673 6675

4.51 5.14 4.41 4.14 3.71 2.61 2.32

0.1603 0.2063 0.2213 0.2316 0.2384 0.2392 0.2356

111.59 1771 2221 2612 2962 3862 4802

comparing to TraPPE and NERD. The SAFT-γ WCA EOS reflects all of these features while it maintains very satisfactory agreement with available simulation data. In TraPPE and NERD force fields, the C–C bond length is 1.54 Å and C–C–C bond angle is 114◦ . In TraPPE, the bond length is kept constant while in NERD, a harmonic potential with a force constant is governing the vibration of bonds. On the other hand, in OPLS-UA, C–C bond length and C–C–C bond angles are 1.53 Å and 112◦ , respectively. However, none of these minor differences is affecting the predictability of the SAFT-γ WCA EOS. For dihedral angle, the studied force fields share the same potential. Nevertheless, Bernard-Brunel and Potoff122 showed that, at least for n-alkanes, the dihedral angle potential has little effect on VLE properties unless molecular topology is artificially constrained to a specific structure. Therefore, it is likely for the SAFT-γ WCA parameters to be applicable to force fields with a different dihedral angle potential than that used in TraPPE, NERD, and OPLS-UA. Orthobaric densities, vapor pressure, second virial coefficient, compressed liquid density, and critical properties are all related to derivatives of Helmholtz free energy with respect to density. It is also necessary to examine EOS predictions for properties related to derivatives with respect to temperature. To do so, isochoric heat capacity at P = 20 MPa is illustrated in Figure 14 for butane and heptane.92 Heat capacity for the TraPPE force field is calculated using the fluctuation of non-

TABLE V. Average absolute deviation of vapor pressure and saturated liquid density for selected n-alkanes. Average absolute deviation (AAD%)a Compound

ρ sat,liq

Psat

N data

T range (K)

Ref.

Methane Ethane Propane Butane Pentane Octane Dodecane

0.97 0.53 1.17 0.34 1.03 1.0 0.7

6.1 4.3 3.6 4.6 5.2 3.1 4.2

6 6 5 5 5 7 5

111–175 178–275 217–344 262–392 298–439 355–535 450–620

97 97 97 97 97 117 97

Average a

|x

SAFT−γ WCA

0.82 −x

TraPPE

|/x

4.4 TraPPE

× (100/N).

bonded potential energy in the NVT ensemble:124 ig

C V − CV =

 1  2  U NV T − U 2NV T . 2 kB T

(51)

Using PackMol,100 300 n-butane and 144 n-heptane molecules were randomly placed in a cubic simulation box with standard periodic boundary condition. To measure the fluctuations at the pressure of interest, the volume of the box was relaxed for 2 ns at NPT ensemble with a time step of 1 fs. Then, the box was further equilibrated for 8 ns at NVT ensemble followed by 8 more ns as the production period. Block averaging is used with a 2 ns time interval. The Velocity-Verlet time integration method was used along with Nose-Hoover101 thermostating and barostating method with a damping factor of 0.2 ps and 1 ps for temperature and pressure, respectively. A spherical cut off distance of 14 Å with standard long range corrections was applied. Heat capacity can also be evaluated directly through its definition: ∂U ig C V − CV = . (52) ∂T NV Potential energy was calculated at constant volume and number of molecules with a temperature interval of 5 K. Then, the first derivative in respect to temperature was taken numerically. The results of Eq. (52) were in satisfactory agreement, within statistical uncertainty, with those of Eq. (51). As it can be seen from Figure 14(a), the TraPPE force field systematically underestimates the residual heat capacity. The same trend was observed at other pressures and for longer chains. It should be noted that the deficiency in theory, or simulation, is better demonstrated when the residual part of heat capacity is compared to experiment (Figure 14, plot (a) vs. (b)). Scaling the small numbers of residual heat capacity by ideal chain part simply camouflages the problems of theory, or simulation, and should be avoided. In our case, doubling the number of molecules and changing the thermostating method to Berendsen did not improve the prediction of TraPPE for residual heat capacity. On the other hand, the SAFT-γ WCA EOS is in better agreement with experimental data, which is simply a fortunate cancellation of error. Ideally, the EOS has to match the simulation data because all the adjustable parameters are fitted to simulation rather than experiment. To understand the

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.220.71.24 On: Mon, 11 Aug 2014 05:20:30

234104-18

A. F. Ghobadi and J. R. Elliott

J. Chem. Phys. 139, 234104 (2013)

FIG. 13. Vapor pressure (a) and orthobaric densities (b) of ethane (C2), butane (C4), and octane (C8) from OPLS-UA (blue circles), NERD (red squares), and TraPPE (black triangles) full potential simulations. All lines are results of SAFT-γ WCA EOS with the corresponding force field parameters. Simulation data of TraPPE for all compounds and OPLS-UA for C2 and C8 are given by Martin and Siepmann.97 OPLS-UA simulation data of C4 are taken from Lísal et al.123 Orthobaric densities of all compounds for NERD force field are provided by Nath et al.98 Vapor pressure of C2 and C8 are taken from Errington and Panagiotopoulos.117 The data in the last reference are read from plots using PlotDigitizer software.118

FIG. 14. (a) Residual isochoric heat capacity, (b) isochoric heat capacity, and (c) non-bonded internal energy of butane and heptane at P = 20 MPa. The predictions of the SAFT-γ WCA EOS for TraPPE force field are illustrated by solid black line. Filled triangles are simulation data of this work. If not shown, error bars are less than the symbol size. Experimental data (unfilled gray squares) are taken from NIST WebBook.115

FIG. 15. (a) Pressure-composition phase diagram of ethane-heptane mixture given by TraPPE (Martin et al.,126 black filled triangle) and NERD (Nath et al.,127 red filled squares) force fields. The prediction of SAFT-γ WCA EOS is given by solid lines with corresponding colors. (b) Temperature-composition diagram of octane-dodecane mixture predicted by a preliminary version of TraPPE force field (Martin and Siepmann,125 grey filled triangle) and SAFT-γ WCA EOS (grey solid line).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.220.71.24 On: Mon, 11 Aug 2014 05:20:30

234104-19

A. F. Ghobadi and J. R. Elliott

J. Chem. Phys. 139, 234104 (2013)

source of this deficiency, we found it useful to plot internal energy versus 1/Tr where Tr is reduced temperature (T/Tc ). Figure 14(c) shows that the potential energy is in qualitative agreement with simulation data, especially in the temperature range that Helmholtz free energy of butane was fitted to simulation. However, the slope of U versus 1/Tr , which is proportional to the second derivative of Helmholtz energy in respect to temperature, is overestimated. This overestimation suggests that the description of theory for TraPPE chains is too soft although n∗ was as high as 24 for CH2 . One may be able to improve the consistency between theory and simulation by adjusting the temperature dependency of the effective diameter or by introducing a temperature-dependent shape factor. As a final assessment, we show the ability of the SAFTγ WCA EOS in the prediction of VLE properties for binary mixtures. Figure 15(a) shows the P-x phase diagram for ethane-heptane mixture at T = 366 K as predicted by TraPPE and NERD. Figure 15(b) illustrates T-x diagram of octane-dodecane mixture at P = 0.02 MPa as predicted by a preliminary version of TraPPE force field.125 The agreement between simulation and the SAFT-γ WCA EOS is very satisfactory without using any binary interaction parameter. In contrast to the original SAFT-γ EOS, we directly adopt the force field energy parameter (ε) into our formulation of the SAFT-γ WCA EOS with no further modifications for cross group interactions.

VII. CONCLUSION

In this work, we could show that a robust comparison between theory and simulation can lead to an EOS that is consistent with molecular simulation of transferable unitedatom force fields. A direct relation between theory and simulation, as proposed in this work, provides a meaningful bridge between the united-atom description of molecules and softfused heteronuclear chains adapted in the SAFT approach. For example, we showed that the thermodynamics of the reference fluids at several temperatures and densities can be collapsed onto a single density dependent curve for each chain length. Furthermore, six transferable parameters suffice to generate these curves for all chain lengths. We also showed that the perturbation contributions for n-alkanes can also be predicted with a similar procedure. These perturbations include the third order term, which enables notable accuracy in the critical region. We achieve accurate characterization of all perturbation contributions at all temperatures, densities, and chain lengths using an “effective diameter” that is independent of density. This greatly simplifies the EOS and all its derivatives with respect to density when using WCA TPT.

f3 =

The resulting EOS takes force field parameters as input and reproduces simulation data for VLE and thermophysical properties. The EOS inherits many advantages of molecular simulation while it facilitates engineering applications with much shorter time scales than full potential molecular simulation. The direct correspondence between the EOS and simulation can be recognized as a valuable tool to further optimize force field parameters in an efficient way. Furthermore, the terms of the EOS are formulated in a manner that facilitates projections to conditions of density and temperature that are thermodynamically unstable. These conditions are important to the application of classical density functional theory for predicting interfacial properties. We report the investigation of interfacial properties in Paper II of this series.128 APPENDIX: CHEMICAL POTENTIAL COMPUTATION

In this appendix, we present the required equations for chemical potential calculation based on SAFT-γ WCA EOS. We define Helmholtz free energy density as A A = ρ = aρ, (A1) f = V kB T N kB T where a is given by Eq. (12) and ρ is number density. With this definition, the chemical potential (μ) for component i in the mixture would be μi ∂f = μ˜ i = . (A2) kB T ∂ρi Therefore, one can write   WCA μ˜ 1i μ˜ 2i μ˜ 3i pert res ch + 2 + 3 , μ˜ i = μ˜ 0i + μ˜ i = μ˜ i + μ˜ i + T T T where the superscript res denotes the residual chemical potential and subscripts 1–3 represent first to third order perturbation contributions. From Eq. (21) and for the WCA sphere term:   ∂ f1 n0 + f2 n1 n2 + f3 n32 WCA = , (A3) μ˜ i ∂ρi which can be expanded as  ∂n2 ∂n0 ∂n1  + f2 n2 + f2 n1 + 3f3 n22 ∂ρi ∂ρi ∂ρi ∂f1 ∂f2 ∂f3 3 ∂n3 + n0 + n1 n2 + n , (A4) ∂n3 ∂n3 ∂n3 2 ∂ρi

μ˜ WCA = f1 i

f1 = − ln(1 − n3 ), f2 =

γ1 , 1 − n3

  (3γ2 − 2γ3 ) n3 − 1.5 (γ2 − γ3 ) 2n3 − n23 + γ3 (1 − n3 )2 ln (1 − n3 ) 36π n23 (1 − n3 )2

(A5) (A6)

(A7)

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.220.71.24 On: Mon, 11 Aug 2014 05:20:30

234104-20

A. F. Ghobadi and J. R. Elliott

J. Chem. Phys. 139, 234104 (2013)

∂nω = M(ω) i , ∂ρi

And the derivatives with respect to density are ∂f1 1 = , ∂n3 1 − n3

(A8)

∂f1 γ1 = , ∂n3 (1 − n3 )2

(A9)

  f¯3 − 36πf3 2n3 (1 − n3 )2 − 2n23 (1 − n3 ) ∂f3 = , (A10) ∂n3 36π n23 (1 − n3 )2 f¯3 = (3γ2 − 2γ3 ) − γ3 (1 − n3 )(2 ln(1 − n3 ) + 1) − 3(γ2 − γ3 )(1 − n3 )

∂gijWCA ∂ρi

=

(A11)

(A12)

where ω is a generic index and M(ω) is defined in Eq. (19). For the chain contribution, and based on Eq. (24), we have  NS  1 ch μ˜ i = − νk,i sk N Bk ln giiWCA 2 k=1  NS  NC WCA  1 1 ∂gjj − ρj νk,j sk N Bk , (A13) WCA ∂ρ 2 k=1 gjj i j =1 where gWCA is given by Eq. (25) and its first derivative in respect to density is

¯ ¯ dii djj (γ1 + γ2 ) (1 − n3 ) ∂n2 / ∂ρi + 2n2 ∂n3 / ∂ρi ∂n3 1 + 2 ∂ρ 3 ¯ 4 dii + d¯jj (1 − n3 ) (1 − n3 ) i ¯ ¯ 2 dii djj (3γ2 − γ3 ) 2n2 (1 − n3 ) ∂n2 / ∂ρi + 3n22 ∂n3 / ∂ρi + , 4 ¯ 36 dii + d¯jj (1 − n3 )

(A14)

where the effective size parameter of component i (d¯ii ) is defined in Eq. (26). Equations (A3)–(A14) complete the chemical potential of the repulsive reference moiety. The chemical potential for the first order perturbation contribution would be   NC NS NS     eff1  ∂ (ρA1 ) WCA1 WCA =− νk,i νl,j sk sl αkl g0,kl n3,kl 2ρj μ˜ 1i = ∂ρi j =1 k=1 l=1 −

NC  NC 

ρj ρw

j =1 w=1

 eff1  WCA ∂g0,kl n3,kl ∂ρi

NS NS  

WCA1 νk,j νl,w sk sl αkl

k=1 l=1

 eff1  WCA n3,kl ∂g0,kl ∂ρi

,

 eff1  eff1 WCA n3,kl ∂n3,kl ∂n3 ∂g0,kl = . eff1 ∂n3 ∂ρi ∂n3,kl

(A15)

(A16)

Using Eq. (29):  eff1  WCA n3,kl ∂g0,kl ∂neff1 3,kl

     eff1 2 eff1 eff1 c1 − 2c2 neff1 3,kl − 1 − n3,kl + 3 4 + c1 n3,kl − c2 n3,kl = ,  4 4 1 − neff1 3,kl

(A17)

where c1 and c2 are model constants: c1 = 3 (γ1 + γ2 ) − 8,

(A18)

c2 = 3γ1 + γ3 − 4. Finally, according to Eqs. (30) and (31): ∂neff1 3,kl ∂n3 ∂neff1 3,w  ∂n3

=

∂neff1 3,k ∂n3



(σk εk ) +

∂neff1 3,l ∂n3



(σl εl )

,

(A19)

w  = 1 to N S,

(A20)

1 / (σk εk ) + 1 / (σl εl )

= aw13 + 2aw13 n3 ,

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.220.71.24 On: Mon, 11 Aug 2014 05:20:30

234104-21

A. F. Ghobadi and J. R. Elliott

J. Chem. Phys. 139, 234104 (2013)

which completes the definition of chemical potential for the first order perturbation term. For the second order term we can write   NS NC NS     eff2  ∂ (ρA2 ) WCA WCA2 WCA = −K νk,i νl,j sk sl αkl g0,kl n3,kl μ˜ 2i = 2ρj ∂ρi j =1 k=1 l=1 −K

WCA

NC  NC 

ρj ρw

j =1 w=1



NS NS  

WCA2 νk,j νl,w sk sl αkl

 eff2  WCA n3,kl ∂g0,kl ∂ρi

k=1 l=1

NC NC NS NS     ∂K WCA   WCA2 WCA ρi ρj νk,i νl,j sk sl αkl g0,kl neff2 3,kl . ∂ρi i=1 j =1 k=1 l=1

(A21)

WCA eff2 The derivative of g0,kl (n3,kl ) in respect to density resembles that of A1 term (Eqs. (A16)–(A20)). For the isothermal compressibility of WCA fluid, Eq. (36) leads to

and

∂n3 ∂ρi

∂K WCA ∂K WCA ∂n3 = , ∂ρi ∂n3 ∂ρi

(A22)

∂K WCA 4 (1 − n3 )3 + K WCA (∂X / ∂n3 ) , =− ∂n3 X

(A23)

X = 1 + n3 (6γ1 − 2) + n23 (9γ2 − 6γ1 + 1) − 4γ3 n33 + γ3 n43 ,

(A24)

∂X = 6γ1 − 2 + 2n3 (9γ2 − 6γ1 + 1) − 12γ3 n23 + 4γ3 n33 , ∂n3

(A25)

is given by Eq. (A12). Eventually, for the third order term we have   NS NS  NC    eff2  ∂ (ρA3 ) WCA 2 2 WCA2 WCA μ˜ 3i = = −K νk,i νl,j sk sl αkl g0,kl n3,kl kl 2ρj ∂ρi j =1 k=1 l=1 −K

WCA

NC NC   j =1 w=1



ρj ρw

NS NS  



2 WCA2 νk,i νl,w sk2 sl αkl

 eff2  WCA n3,kl ∂g0,kl

k=1 l=1

∂ρi

kl +

WCA g0,kl

 eff2  ∂ kl n3,kl ∂ρi

NC NC NS NS     ∂K WCA   2 WCA2 WCA ρi ρj νk,i νl,j sk2 sl αkl g0,kl neff2 3,kl kl , ∂ρi i=1 j =1 k=1 l=1

(A26)

and for the derivative of Gaussian function in respect to density, Eq. (42) gives   ∂ k l (σk εk ) + ∂ (σl εl ) ∂ kl ∂ρi ∂ρi = , ∂ρi 1 / (σk εk ) + 1 / (σl εl )     2  (3) ∂ w ∂n3 εw 31  ∂ w Mi = = aw −2aw32 n3 − aw33 exp −aw32 n3 − aw33 ∂ρi ∂n3 ∂ρi kB

1

E. D. Maginn and J. R. Elliott, “Historical perspective and current outlook for molecular dynamics as a chemical engineering tool,” Ind. Eng. Chem. Res. 49, 3059–3078 (2010). 2 W. G. Chapman, G. Jackson, and K. E. Gubbins, “Phase equilibria of associating fluids,” Mol. Phys. 65, 1057–1079 (1988). 3 W. G. Chapman, K. E. Gubbins, G. Jackson, M. Radosz, “SAFT: Equation of state solution model for associating fluids,” Fluid Phase Equilib. 52, 31 (1989). 4 W. G. Chapman, K. E. Gubbins, G. Jackson, and M. Radosz, “New reference equation of state for associating liquids,” Ind. Eng. Chem. Res. 29(8), 1709–1721 (1990). 5 G. Jackson, W. Chapman, and K. Gubbins, “Phase equilibria of associating fluids of spherical and chain molecules,” Int. J. Thermophys. 9(5), 769– 779 (1988).



(A27) w  = 1 to N S.

(A28)

6 M.

S. Wertheim, “Fluids with highly directional attractive forces. I. statistical thermodynamics,” J. Stat. Phys. 35, 19 (1984). 7 M. S. Wertheim, “Fluids with highly directional attractive forces. II. thermodynamic perturbation theory and integral equation,” J. Stat. Phys. 35, 35 (1984). 8 M. S. Wertheim, “Fluids with highly directional attractive forces. IV. Equilibrium polymerization,” J. Stat. Phys. 42, 477 (1986). 9 M. S. Wertheim, “Fluids with highly directional attractive forces. III. multiple attraction sites,” J. Stat. Phys. 42, 459 (1986). 10 E. A. Müller and K. E. Gubbins, “Molecular-based equations of state for associating fluids: A review of SAFT and related approaches,” Ind. Eng. Chem. Res. 40, 2193 (2001). 11 P. K. Jog, A. Garcia-Cuellar, and W. G. Chapman, “Extensions and applications of the SAFT equation of state to solvents,

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.220.71.24 On: Mon, 11 Aug 2014 05:20:30

234104-22

A. F. Ghobadi and J. R. Elliott

monomers, and polymers,” Fluid Phase Equilib. 158, 321–326 (1999). 12 P. Paricaud, A. Galindo, and G. Jackson, “Recent advances in the use of the SAFT approach in describing electrolytes, interfaces, liquid crystals and polymers,” Fluid Phase Equilib. 194, 87–96 (2002). 13 I. Polishuk and A. Mulero, “The numerical challenges of SAFT EoS models,” Rev. Chem. Eng. 27(5), 241 (2011). 14 A. Lymperiadis, C. S. Adjiman, A. Galindo, and G. Jackson, “A group contribution method for associating chain molecules based on the statistical associating fluid theory (SAFT-γ ),” J. Chem. Phys. 127, 234903 (2007). 15 A. Lymperiadis, C. S. Adjiman, G. Jackson, and A. Galindo, “A generalisation of the SAFT-group contribution method for groups comprising multiple spherical segments,” Fluid Phase Equilib. 274(1), 85–104 (2008). 16 D. Chandler, J. D. McCoy, and S. J. Singer, “Density functional theory of nonuniform polyatomic systems. I. General formulation,” J. Chem. Phys. 85, 5971 (1986). 17 W. E. McMullen and K. F. Freed, “A density functional theory of polymer phase transitions and interfaces,” J. Chem. Phys. 92, 1413 (1990). 18 C. E. Woodward, “A density functional theory for polymers: Application to hard chain–hard sphere mixtures in slitlike pores,” J. Chem. Phys. 94, 3183 (1991). 19 E. Kierlik and M. Rosinberg, “A perturbation density-functional theory for polyatomic fluids. I. Rigid molecules,” J. Chem. Phys. 97, 9222 (1992). 20 J. D. Weeks, D. Chandler, H. C. Andersen, “Role of repulsive forces in determining the equilibrium structure of simple liquids,” J. Chem. Phys. 54(12), 5237 (1971). 21 C. Avendaño, T. Lafitte, A. Galindo, C. S. Adjiman, G. Jackson, and E. A. Müller, “SAFT-γ force field for the simulation of molecular fluids. 1. A single-site coarse grained model of carbon dioxide,” J. Phys. Chem. B 115(38), 11154–11169 (2011). 22 C. Avendano, T. Lafitte, C. S. Adjiman, A. Galindo, E. A. Müller, and G. Jackson, “SAFT-γ force field for the simulation of molecular fluids. 2. Coarse-grained models of greenhouse gases, refrigerants, and long alkanes,” J. Phys. Chem. B 117, 2717 (2013). 23 L. Verlet, “Computer,” Phys. Rev. 165(1), 201 (1968). 24 J. A. Barker and D. Henderson, “Perturbation theory and equation of state for fluids: The square-well potential,” J. Chem. Phys. 47, 2856 (1967). 25 J. A. Barker and D. Henderson, “Perturbation theory and equation of state for fluids. II. successful theory of liquids,” J. Chem. Phys. 47, 4714 (1967). 26 R. W. Zwanzig, “High-temperature equation of state by a perturbation method. I. Nonpolar gases,” J. Chem. Phys. 22, 1420 (1954). 27 Y. T. Pavlyukhin, “Thermodynamic perturbation theory of simple liquids: Monte Carlo simulation of a hard sphere system and the Helmholtz free energy of SW fluids,” J. Struct. Chem. 53(3), 476–486 (2012). 28 J. P. Hansen and I. R. McDonald, Theory of Simple Liquids (Academic Press, 2006). 29 S. Zhou and J. Solana, “Progress in the perturbation approach in fluid and fluid-related theories,” Chem. Rev. 109(6), 2829–2858 (2009). 30 R. Roth, R. Evans, A. Lang, and G. Kahl, “Fundamental measure theory for hard-sphere mixtures revisited: The White Bear version,” J. Phys.: Condens. Matter 14(46), 12063 (2002). 31 G. J. Gloor, G. Jackson, F. J. Blas, E. M. Del Río, and E. de Miguel, “An accurate density functional theory for the vapor-liquid interface of associating chain molecules based on the statistical associating fluid theory for potentials of variable range,” J. Chem. Phys. 121, 12740 (2004). 32 J. Wu, “Density functional theory for chemical engineering: From capillarity to soft materials,” AIChE J. 52(3), 1169–1193 (2006). 33 J. Gross, “A density functional theory for vapor-liquid interfaces using the PCP-SAFT equation of state,” J. Chem. Phys. 131, 204705 (2009). 34 N. F. Carnahan and K. E. Starling, “Equation of state for nonattracting rigid spheres,” J. Chem. Phys. 51, 635 (1969). 35 S. Tripathi and W. G. Chapman, “Microstructure and thermodynamics of inhomogeneous polymer blends and solutions,” Phys. Rev. Lett. 94(8), 087801 (2005). 36 S. Tripathi and W. G. Chapman, “Microstructure of inhomogeneous polyatomic mixtures from a density functional formalism for atomic mixtures,” J. Chem. Phys. 122, 094506 (2005). 37 A. Dominik, S. Tripathi, and W. G. Chapman, “Bulk and interfacial properties of polymers from interfacial SAFT density functional theory,” Ind. Eng. Chem. Res. 45(20), 6785–6792 (2006). 38 S. Jain, A. Dominik, and W. G. Chapman, “Modified interfacial statistical associating fluid theory: A perturbation density functional theory for inhomogeneous complex fluids,” J. Chem. Phys. 127, 244904 (2007).

J. Chem. Phys. 139, 234104 (2013) 39 C.

P. Emborsky, K. R. Cox, and W. G. Chapman, “Exploring parameter space effects on structure-property relationships of surfactants at liquidliquid interfaces,” J. Chem. Phys 135, 084708 (2011). 40 B. D. Marshall, C. Emborsky, K. Cox, and W. G. Chapman, “Effect of bond rigidity and molecular structure on the self-assembly of amphiphilic molecules using second-order classical density functional theory,” J. Phys. Chem. B 116(9), 2730–2738 (2012). 41 S. Jain, P. Jog, J. Weinhold, R. Srivastava, and W. G. Chapman, “Modified interfacial statistical associating fluid theory: Application to tethered polymer chains,” J. Chem. Phys. 128, 154910 (2008). 42 A. Bymaster, S. Jain, and W. G. Chapman, “Microstructure and depletion forces in polymer-colloid mixtures from an interfacial statistical associating fluid theory,” J. Chem. Phys. 128, 164910 (2008). 43 S. Jain and W. Chapman, “Effect of confinement on the ordering of symmetric diblock copolymers: Application of interfacial statistical associating fluid theory,” Mol. Phys. 107(1), 1–17 (2009). 44 B. D. Marshall, K. R. Cox, and W. G. Chapman, “Supramolecular assembly and surfactant behavior of triblock rod-coil amphiphiles at liquid interfaces using classical density functional theory,” Soft Matter 2012. 45 A. Bymaster and W. G. Chapman, “An i SAFT density functional theory for associating polyatomic molecules,” J. Phys. Chem. B 114(38), 12298– 12307 (2010). 46 L. J. D. Frink, A. L. Frischknecht, M. A. Heroux, M. L. Parks, and A. G. Salinger, “Toward quantitative coarse-grained models of lipids with fluids density functional theory,” J. Chem. Theory Comput. 8(4), 1393–1408 (2012). 47 B. D. Marshall and W. G. Chapman, “A density functional theory for patchy colloids based on Wertheim’s association theory: Beyond the single bonding condition,” J. Chem. Phys. 138, 044901 (2013). 48 B. D. Marshall and W. G. Chapman, “Higher order classical density functional theory for branched chains and rings,” J. Phys. Chem. B 115(50), 15036–15047 (2011). 49 B. D. Marshall, K. R. Cox, and W. G. Chapman, “A classical density functional theory study of the neat N-alkane/water interface,” J. Phys. Chem. C 116, 17641 (2012). 50 G. Jackson and K. E. Gubbins, “Mixtures of associating spherical and chain molecules,” Pure Appl. Chem. 61, 1021 (1989). 51 R. Cotterman, B. Schwarz, and J. Prausnitz, “Molecular thermodynamics for fluids at low and high densities. Part I: Pure fluids containing small or large molecules,” AIChE J. 32(11), 1787–1798 (1986). 52 Lafitte et al., “Accurate statistical associating fluid theory for chain molecules formed from Mie segments,” J. Chem. Phys. 139, 154504 (2013). 53 S. Huang and M. Radosz, “Equation of state for small, large, polydisperse, and associating molecules,” Ind. Eng. Chem. Res. 29, 2284 (1990). 54 S. Huang and M. Radosz, “Equation of state for small, large, polydisperse, and associating molecules: Extension to fluid mixtures,” Ind. Eng. Chem. Res. 30, 1994 (1991). 55 J. K. Johnson, E. A. Müller, and K. E. Gubbins, “Equation of state for Lennard-Jones chains,” J. Phys. Chem. 98(25), 6413–6419 (1994). 56 J. K. Johnson, J. A. Zollweg, and K. E. Gubbins, “The Lennard-Jones equation of state revisited,” Mol. Phys. 78(3), 591–618 (1993). 57 E. A. Müller and K. E. Gubbins, “An equation of state for water from a simplified intermolecular potential,” Ind. Eng. Chem. Res. 34(10), 3662– 3673 (1995). 58 T. Kraska and K. E. Gubbins, “Phase equilibria calculations with a modified SAFT equation of state. 1. Pure alkanes, alkanols, and water,” Ind. Eng. Chem. Res. 35(12), 4727–4737 (1996). 59 F. J. Blas and L. F. Vega, “Thermodynamic behaviour of homonuclear and heteronuclear Lennard-Jones chains with association sites from simulation and theory,” Mol. Phys. 92(1), 135–150 (1997). 60 J. Kolafa and I. Nezbeda, “The Lennard-Jones fluid: An accurate analytic and theoretically-based equation of state,” Fluid Phase Equilib. 100, 1–34 (1994). 61 M. Banaszak, Y. Chiew, and M. Radosz, “Thermodynamic perturbation theory: Sticky chains and square-well chains,” Phys. Rev. E 48, 3760– 3765 (1993). 62 M. Banaszak, Y. Chiew, R. O’Lenick, and M. Radosz, “Thermodynamic perturbation theory: Lennard-Jones chains,” J. Chem. Phys. 100(5), 3803– 3807 (1994). 63 J. Gross and G. Sadowski, “Application of perturbation theory to a hardchain reference fluid: An equation of state for square-well chains,” Fluid Phase Equilib. 168(2), 183–199 (2000).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.220.71.24 On: Mon, 11 Aug 2014 05:20:30

234104-23 64 J.

A. F. Ghobadi and J. R. Elliott

Gross and G. Sadowski, “Perturbed-chain SAFT: An equation of state based on a perturbation theory for chain molecules,” Ind. Eng. Chem. Res. 40, 1244–1260 (2001). 65 Y. C. Chiew, “Percus-Yevick integral equation theory for athermal hardsphere chains. II. Average intermolecular correlation functions,” Mol. Phys. 73(2), 359–73 (1991). 66 A. Gil-Villegas, A. Galindo, P. J. Whitehead, S. J. Mills, G. Jackson, and A. N. Burgess, “Statistical associating fluid theory for chain molecules with attractive potentials of variable range,” J. Chem. Phys. 106, 4168 (1997). 67 L. Davies, A. Gil-Villegas, and G. Jackson, “Describing the properties of chains of segments interacting via soft-core potentials of variable range with the SAFT-VR approach,” Int. J. Thermophys. 19(3), 675–686 (1998). 68 T. Lafitte, D. Bessieres, M. M. Piñeiro, and J.-L. Daridon, “Simultaneous estimation of phase behavior and second-derivative properties using the statistical associating fluid theory with variable range approach,” J. Chem. Phys. 124, 024509 (2006). 69 G. J. Gloor, F. J. Blas, E. M. N. del Río, E. de Miguel, and G. Jackson, “A SAFT–DFT approach for the vapour–liquid interface of associating fluids,” Fluid Phase Equilib. 194, 521–530 (2002). 70 F. Llovell, A. Galindo, F. J. Blas, and G. Jackson, “Classical density functional theory for the prediction of the surface tension and interfacial properties of fluids mixtures of chain molecules based on the statistical associating fluid theory for potentials of variable range,” J. Chem. Phys. 133, 024704 (2010). 71 G. J. Gloor, G. Jackson, F. Blas, E. M. Del Rio, and E. De Miguel, “Prediction of the vapor-liquid interfacial tension of nonassociating and associating fluids with the SAFT-VR density functional theory,” J. Phys. Chem. C 111(43), 15513–15522 (2007). 72 F. Llovell, N. Mac Dowell, F. J. Blas, A. Galindo, and G. Jackson, “Application of the SAFT-VR density functional theory to the prediction of the interfacial properties of mixtures of relevance to reservoir engineering,” Fluid Phase Equilib. 336, 137 (2012). 73 S. Tamouza, J. Passarello, P. Tobaly, and J. de Hemptinne, “Group contribution method with SAFT EOS applied to vapor liquid equilibria of various hydrocarbon series,” Fluid Phase Equilib. 222, 67–76 (2004). 74 T. X. Nguyen Thi, S. Tamouza, P. Tobaly, J.-P. Passarello, and J.-C. de Hemptinne, “Application of group contribution SAFT equation of state (GC-SAFT) to model phase behaviour of light and heavy esters,” Fluid Phase Equilib. 238(2), 254–261 (2005). 75 D. N. Huynh, M. Benamira, J. Passarello, P. Tobaly, and J. de Hemptinne, “Application of GC-SAFT EOS to polycyclic aromatic hydrocarbons,” Fluid Phase Equilib. 254(1), 60–66 (2007). 76 Y. Peng, K. D. Goff, M. dos Ramos, and C. McCabe, “Developing a predictive group-contribution-based SAFT-VR equation of state,” Fluid Phase Equilib. 277(2), 131–144 (2009). 77 Y. Peng, K. D. Goff, M. C. dos Ramos, and C. McCabe, “Predicting the phase behavior of polymer systems with the GC-SAFT-VR approach,” Ind. Eng. Chem. Res. 49(3), 1378–1394 (2010). 78 V. Papaioannou, C. S. Adjiman, G. Jackson, and A. Galindo, “Simultaneous prediction of vapour–liquid and liquid–liquid equilibria (VLE and LLE) of aqueous mixtures with the SAFT-γ group contribution approach,” Fluid Phase Equilib. 306(1), 82–96 (2011). 79 J. A. Barker and D. Henderson, “What is liquid?,” Rev. Mod. Phys. 48, 587–671 (1976). 80 Y. Tang, “Role of the Barker–Henderson diameter in thermodynamics,” J. Chem. Phys. 116, 6694 (2002). 81 J. Cui and J. R. Elliott, J. Chem. Phys. 116, 8625 (2002). 82 J. R. Elliott, “Optimized step potential models for n-alkanes and benzene,” Fluid Phase Equilib. 194(197), 161–168 (2002). 83 J. R. Elliott, Z. N. Gerek, and N. Gray, “Combining molecular dynamics and chemical process simulation: The SPEAD model,” Asia Pac. J. Chem. Eng. 2(4), 257–271 (2007). 84 A. Vahid and J. R. Elliott, “Transferable intermolecular potentials for carboxylic acids and their phase behavior,” AIChE J. 56(2), 485–505 (2010). 85 D. C. Rapaport, The Art of Molecular Dynamics Simulation (Cambridge University Press, 2004). 86 A. F. Ghobadi and J. R. Elliott, “Evaluating perturbation contributions in SAFT models by comparing to molecular simulation of n-alkanes,” Fluid Phase Equilib. 306(1), 57–66 (2011). 87 W. R. Smith, D. Henderson, and J. A. Barker, “Approximate evaluation of the second-order term in the perturbation theory of fluids,” J. Chem. Phys. 53, 508 (1970).

J. Chem. Phys. 139, 234104 (2013) 88 S. Zhou, and J. Solana, “Comprehensive investigation about the second or-

der term of thermodynamic perturbation expansion,” J. Chem. Phys. 131, 134106 (2009). 89 J. M. C. Sanchez, T. Danner, and J. Gross, “Grand canonical Monte Carlo simulations of vapor-liquid equilibria using a bias potential from an analytic equation of state,” J. Chem. Phys. 138, 234106 (2013). 90 P. Paricaud, “A general perturbation approach for equation of state development: Applications to simple fluids, ab initio potentials, and fullerenes,” J. Chem. Phys. 124(15), 154505 (2006). 91 T. van Westen, T. J. Vlugt, and J. Gross, “Determining force field parameters using a physically based equation of state,” J. Phys. Chem. B 115(24), 7872–7880 (2011). 92 See supplementary material at http://dx.doi.org/10.1063/1.4838457 for tabulated data pertaining to points in the figures. 93 F. Betancourt-Cárdenas, L. Galicia-Luna, and S. Sandler, “Equation of state for the Lennard–Jones fluid based on the perturbation theory,” Fluid Phase Equilib. 264(1), 174–183 (2008). 94 J. R. Elliott and N. H. Gray, “Asymptotic trends in thermodynamic perturbation theory,” J. Chem. Phys. 123, 184902 (2005). 95 J. Cui and J. R. Elliott, “Phase envelopes for variable well-width square well chain fluids,” J. Chem. Phys. 114, 7283 (2001). 96 A. F. Ghobadi, V. Taghikhani, and J. R. Elliott, “Investigation on the solubility of SO2 and CO2 in imidazolium-based ionic liquids using NPT Monte Carlo simulation,” J. Phys. Chem. B 115(46), 13599–13607 (2011). 97 M. G. Martin and J. I. Siepmann, “Transferable potentials for phase equilibria. 1. United-atom description of n-alkanes,” J. Phys. Chem. B 102, 2569–2577 (1998). 98 S. K. Nath, F. A. Escobedo, and J. J. de Pablo, “On the simulation of vaporliquid equilibria for alkanes,” J. Chem. Phys 108, 9905–9911 (1998). 99 S. Plimpton, “Fast parallel algorithms for short-range molecular dynamics,” J. Comput. Phys. 117(1), 1–19 (1995). 100 L. Martínez, R. Andrade, E. Birgin, and J. Martínez, “Packmol: A package for building initial configurations for molecular dynamics simulations,” J. Comput. Chem. 30(13), 2157–2164 (2009). 101 D. J. Evans and B. L. Holian, “The Nose–Hoover thermostat,” J. Chem. Phys. 83, 4069 (1985). 102 E. A. Müller and A. S. Mejía, “Comparison of united-atom potentials for the simulation of vapor–liquid equilibria and interfacial properties of longchain n-alkanes up to n-C100,” J. Phys. Chem. B 115(44), 12822–12834 (2011). 103 W. L. Jorgensen, J. D. Madura, and C. J. Swenson, “Optimized intermolecular potential functions for liquid hydrocarbons,” J. Am. Chem. Soc. 106(22), 6638 (1984). 104 Á. Mulero, Theory and Simulation of Hard-Sphere Fluids and Related Systems (Springer, 2008), Vol. 753. 105 D. M. Heyes and H. Okumura, “Equation of state and structural properties of the Weeks-Chandler-Andersen fluid,” J. Chem. Phys. 124, 164507 (2006). 106 L. Boltzmann, Lectures on Gas Theory (University of California Press, 1964). 107 J. R. Elliott and T. E. Daubert, “The temperature dependence of the hard sphere diameter,” Fluid Phase Equilib. 31, 153–160 (1986). 108 T. Boublik, “Hard sphere equation of state,” J. Chem. Phys. 53, 471–472 (1970). 109 G. A. Mansoori, N. F. Carnahan, K. E. Starling, and T. W. Leland, “Equilibrium thermodynamic properties of the mixture of hard spheres,” J. Chem. Phys. 54, 1523–1525 (1971). 110 Y. Zhou, C. K. Hall, and G. Stell, “Thermodynamic perturbation theory for fused hard-sphere and hard-disk chain fluids,” J. Chem. Phys. 103, 2688 (1995). 111 E. A. Müller, K. E. Gubbins, D. M. Tsangaris, and J. J. de Pablo, “Comment on the accuracy of Wertheim’s theory of associating fluids,” J. Chem. Phys. 103, 3868 (1995). 112 M. D. Amos and G. Jackson, “Bonded hard-sphere (BHS) theory for the equation of state of fused hard-sphere polyatomic molecules and their mixtures,” J. Chem. Phys. 96, 4604 (1992). 113 A. J. Schultz and D. A. Kofke, “Virial coefficients of model alkanes,” J. Chem. Phys. 133, 104101 (2010). 114 J. Barker, P. Leonard, and A. Pompe, “Fifth virial coefficients,” J. Chem. Phys 44(11), 4206–4211 (1966). 115 NIST, NIST Chemistry WebBook, 2005, see http://webbook.nist.gov/chemistry/fluid/. 116 T. Lafitte, C. Avendaño, V. Papaioannou, A. Galindo, C. S. Adjiman, G. Jackson, and E. A. Müller, “SAFT-γ force field for the simulation of

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.220.71.24 On: Mon, 11 Aug 2014 05:20:30

234104-24

A. F. Ghobadi and J. R. Elliott

molecular fluids: 3. Coarse-grained models of benzene and hetero-group models of n-decylbenzene,” Mol. Phys. 110(11–12), 1189–1203 (2012). 117 J. R. Errington and A. Z. Panagiotopoulos, “A new intermolecular potential model for the n-alkane homologous series,” J. Phys. Chem. B 103, 6314 (1999). 118 Plot Digitizer is a Java program used to digitize scanned plots of functional data, for more information, see http://plotdigitizer.sourceforge.net/. 119 DIPPR, Design Institute for Physical Property Data, DIADEM 2004, 2004. 120 Y. Meurice, R. Perry, and S.-W. Tsai, “New applications of the renormalization group method in physics: A brief introduction,” Philos. Trans. R. Soc. London, Ser. A 369(1946), 2602–2611 (2011). 121 E. Forte, F. Llovell, L. F. Vega, J. P. M. Trusler, and A. Galindo, “Application of a renormalization-group treatment to the statistical associating fluid theory for potentials of variable range (SAFT-VR),” J. Chem. Phys. 134, 154102 (2011). 122 D. A. Bernard-Brunel and J. J. Potoff, “Effect of torsional potential on the predicted phase behavior of i n/i-alkanes,” Fluid Phase Equilib. 279(2), 100–104 (2009).

J. Chem. Phys. 139, 234104 (2013) 123 M.

Lisal, W. R. Smith, and I. Nezbeda, “Accurate computer simulation of phase equilibrium for complex fluid mixtures. Application to binaries involving isobutene, methanol, methyl tert-butyl ether, and n-butane,” J. Phys. Chem. B 103(47), 10496 (1999). 124 M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Clarendon Press, 1999). 125 M. G. Martin and J. I. Siepmann, “Predicting multicomponent phase equilibria and free energies of transfer for alkanes by molecular simulation,” J. Am. Chem. Soc. 119(38), 8921–8924 (1997). 126 M. G. Martin, B. Chen, and J. I. Siepmann, “Molecular structure and phase diagram of the binary mixture of n-heptane and supercritical ethane: A Gibbs ensemble Monte Carlo study,” J. Phys. Chem. B 104(10), 2415– 2423 (2000). 127 S. K. Nath, F. A. Escobedo, J. J. de Pablo, and I. Patramai, “Simulation of vapor-liquid equilibria for alkane mixtures,” Ind. Eng. Chem. Res. 37(8), 3195–3202 (1998). 128 A. F. Ghobadi and J. R. Elliott, “Adapting SAFT-gamma perturbation theory to site-based molecular dynamics simulation. II. Confined fluids and vapor-liquid interfaces,” J. Chem. Phys. (unpublished).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.220.71.24 On: Mon, 11 Aug 2014 05:20:30

Adapting SAFT-γ perturbation theory to site-based molecular dynamics simulation. I. Homogeneous fluids.

In this work, we aim to develop a version of the Statistical Associating Fluid Theory (SAFT)-γ equation of state (EOS) that is compatible with united-...
3MB Sizes 0 Downloads 0 Views