Extension of the SAFT-VR Mie EoS To Model Homonuclear Rings and Its Parametrization Based on the Principle of Corresponding States.

The statistical associating fluid theory of variable range employing a Mie potential (SAFT-VR-Mie) proposed by Lafitte et al. (J. Chem Phys. 2013, 139, 154504) is one of the latest versions of the SAFT family. This particular version has been shown to have a remarkable capability to connect experimental determinations, theoretical calculations, and molecular simulations results. However, the theoretical development restricts the model to chains of beads connected in a linear fashion. In this work, the capabilities of the SAFT-VR Mie equation of state for modeling phase equilibria are extended for the case of planar ring compounds. This modification proposed replaces the Helmholtz energy of chain formation by an empirical contribution based on a parallelism to the second-order thermodynamic perturbation theory for hard sphere trimers. The proposed expression is given in terms of an extra parameter, χ, that depends on the number of beads, ms, and the geometry of the ring. The model is used to describe the phase equilibrium for planar ring compounds formed of Mie isotropic segments for the cases of ms equals to 3, 4, 5 (two configurations), and 7 (two configurations). The resulting molecular model is further parametrized, invoking a corresponding states principle resulting in sets of parameters that can be used indistinctively in theoretical calculations or in molecular simulations without any further refinements. The extent and performance of the methodology has been exemplified by predicting the phase equilibria and vapor pressure curves for aromatic hydrocarbons (benzene, hexafluorobenzene, toluene), heterocyclic molecules (2,5-dimethylfuran, sulfolane, tetrahydro-2H-pyran, tetrahydrofuran), and polycyclic aromatic hydrocarbons (naphthalene, pyrene, anthracene, pentacene, and coronene). An important aspect of the theory is that the parameters of the model can be used directly in molecular dynamics (MD) simulations to calculate equilibrium phase properties and interfacial tensions with an accuracy that rivals other coarse grained and united atom models, for example, liquid densities, are predicted, with a maximum absolute average deviation of 3% from both the theory and the MD simulations, while the interfacial tension is predicted, with a maximum absolute average of 8%. The extension to mixtures is exemplified by considering a binary system of hexane (chain fluid) and tetrahydro-2H-pyran (ring fluid).


■ INTRODUCTION
The Statistical Associating Fluid Theory (SAFT) equation of state (EoS) is one of the most versatile, advanced, and accurate molecular-based EoS used to predict the thermophysical properties of pure fluids and fluid mixtures. Its formulation is based on Wertheim's first-order thermodynamic perturbation theory (TPT1) 1−4 and expressed as a sum of contributions to the Helmholtz free energy, namely:  (1) where a is the dimensionless Helmholtz energy, defined as a = A/(Nk B T); A being the total Helmholtz energy, N the number of molecules, k B the Boltzmann constant, and T the temperature. a IG is the ideal gas reference value, while the other terms correspond to increased complexity in the model; a MONO adds the monomer (unbounded) contribution for a chain composed of m s tangential segments, a CHAIN accounts for the formation of said chains, while a ASSOC is the Helmholtz energy contribution accounting for intermolecular association. Since its inception within the Gubbins group in the late 80s, several variants of the SAFT EoS have been proposed in order to incrementally improve its accuracy and range of applicability. It is arguably one of the most important advances in the area of correlation and prediction of fluid phase properties of the last decades. For an abridged overview of the SAFT EoS and its most recent versions, the reader is directed to the available reviews in the open literature. 5−9 Without prejudice to other implementations, the SAFT-VR Mie proposed by Lafitte et al. 10 is employed here, as it provides a remarkable route to connect experimental determinations, theoretical calculations, and molecular simulations results. 11 Its success is attributed to two factors, the versatility of the underlying potential and its accuracy in representing said potential. The Mie potential, 12 ϕ Mie ,  (2) where r is the distance between centers of the particles, ε is the minimum energy of interaction, σ is the short distance where the repulsive and attractive contribution of the potential cancel and λ r and λ a are the exponents that characterize ultimate range of the interaction. It represents an "upgrade" of the commonly used Lennard-Jones (λ r = 12; λ a = 6), incorporating a flexibility in the treatment of the attractive range and repulsive nature of the intermolecular potential, whereas the inclusion of a perturbation expansion of the underlying monomer term, a MONO , to third order allows the theory to describe with precision the expected macroscopic behavior of the potential. This latter aspect is particularly important, as it allows the exploitation of the link between the theoretical description by means of an EoS and the direct molecular simulation of the corresponding fluid. The accuracy of this one-to-one correspondence is crucial in order to represent molecules as coarse-grained (CG) beads with molecular parameters provided by the fitting of the SAFT model to experimental data. 13 The procedure essentially defines a top-down coarse grained approach that has been widely successful in modeling both the phase equilibria, interfacial properties, and self-assembly of organic molecules, polymers, surfactants, and mixtures. 14−19 In one of the seminal manuscripts defining the above theory, Laffite et al., 16 recognized the importance of having the appropriate overall shape of the molecular model. TPT1, on which most of the SAFT formulations are based on, adds to the monomer term of the Helmholtz energy an explicit contribution corresponding to the free energy required to form a chain from the corresponding monomer segments. While TPT1 makes no explicit provision for the connectivity of the beads, an analysis of the assumptions made in the theory with regards to the three (and higher body) distribution functions immediately suggests that the theory will be most accurate for linear chains (or completely flexible and extended chains). 20 Molecules formed by beads connecting each other in a nonlinear fashion require a fundamentally different approach, essentially invoking Wertheim's second order perturbation theory (TPT2). 21 If one is to consider phase equilibria applications, whereas the parameters of the EoS are fitted to experimental data, and the resulting analytical expression is used to correlate and extrapolate the bulk thermodynamic properties, the nuances of the shape of the molecule and those of the underlying potential are to a great extent irrelevant. As a consequence of that the various versions of SAFT (and indeed of any other modern EoS) fundamentally provide comparable accuracy among themselves, and all serve as appropriate tools for use in engineering and spreadsheet calculations. 9 A key point raised by Laffite et al., 16 is that the shape of the underlying molecular model, while having modest effect on the bulk thermodynamic properties, has an important role in the accurate description of the structure of the fluid, as probed by molecular simulation. In ref 16, tasked with the challenge of representing a discotic molecule like benzene as a ring of three beads bonded to form an equilateral triangle, it was suggested to replace the a CHAIN term in the SAFT-VR Mie EoS (see eq 1) by an empirical a RING term proposed by Sear and Jackson. 22 The performance of this new model was tested by comparing the phase equilibria (i.e., temperature−density coexistence envelope and Clausius−Clapeyron curve) and interfacial tension results for benzene to alternative models such as a single-segment CG sphere, a geometrically accurate six-segment TraPPE United Atom 23 model and available experimental data. As concluded in ref 16, the three-segment CG ring is an excellent approximation displaying good agreement with experiments for subcritical and critical phase equilibria, interfacial tension, and as an ancillary benefit avoids the premature freezing exhibited by the single bead model. More importantly, however, the ring model displays the correct behavior for the fluid structure, as evidenced by the center of mass pair distribution function. The Laffite et al., 16 approach suffers from two drawbacks: the accuracy with which the Sear and Jackson a RING term represents the potential is limited (hence the need to rescale the EoS parameters when attempting to employ them in simulations) and the absence of information on the actual ring geometry. The latter limitation is particularly important for the case of ring components with more than three segments. In principle, an expression for ring molecules can be obtained by taking the limit of an associating fluid that can either self-associate (as a snake biting its tail) or associate forming rings. This latter approach is discussed recently by Haghmoradi et al. 24 and the reader is referred to said manuscript for a review of the existing alternative approaches.
This work reports a revised expression for the a RING term of the SAFT-VR Mie EoS that allows the parameters to be used within molecular simulations. The proposed model is developed with input from TPT2, adapted and informed by molecular simulations of models of ring molecules. The resulting model is further parametrized, employing a corresponding states correlation, to establish a link between the force field parameters and a selected set of macroscopical properties.
Helmholtz Energy for Rings. The TPT1 theory, upon which SAFT is based, "forbids" the combination of monomers in such a way that a closed ring (or a branched chain) is formed. In spite of this limitation, versions of the SAFT EoS have been used to model rings compounds, by employing the molecular chain length (m s ) as a lumped parameter that absorbs the geometric shape of molecule. An example of this approach was presented by Huang and Radosz, 25 where correlations for the parameters of polynuclear aromatic compounds (PAH) are presented in terms of molecular weights. In this (and similar approaches), the molecular parameters used for these compounds can not be transferred to molecular simulations and, in general, the corresponding values of the chain length m s do not provide any physical link to the actual molecular geometry.
A closer look at the perturbation form of the SAFT EoS reveals that the contribution to the Helmholtz energy of forming a ring structure must be encapsulated in the a CHAIN term of eq 1, which has the general form of where g ref (σ) denotes the contact value of the radial pair distribution function for the reference monomeric fluid. Sear and Jackson 22,26,27 proposed that if within TPT1 C can be loosely interpreted as the contribution to form (m s − 1) covalent-like bonds between m s monomers and takes the form of C = (m s − 1) for a linear chain, then for the case of an m smembered ring of tangential segments, the expression for the Helmholtz energy for rings (a RING ) should be given by eq 3 with C = m s . The background for this ansatz is that in the case of a circular structure an "additional" bond is being made between the first and last member of a chain to form the ring structure. The use of this closure does not always give satisfactory results: the results of Lafitte et al. 16 suggested that the EoS parameters obtained when using this ansatz and fitting experimental data of vapor pressure and liquid density, need to be rescaled in order to be further used in molecular simulations. The fact that the linear version of the theory does not require such scaling implies some degree of inconsistency. Muller and Gubbins 20 explored the extension of TPT1 to second order (TPT2) and obtained rigorous expressions for the EoS for molecular geometries of rings of m s = 3 for hard fluids. These expressions are based on the knowledge of the threebody distribution functions of the hard sphere reference fluid and depend on the density (packing fraction) of the system. One can compare analytically the results of a linear chain of three beads with that of a ring of three beads and extract from the comparison the "equivalent" value of C to be used within TPT1 to accurately represent a ring. Further details are supplied in the Appendix, but it suffices to say that the derivation leads expression for C with the following form where η is the packing fraction, η = m s πρσ 3 /6, ρ is the molecular number density, and χ is a parameter which depends on the geometry and the reference potential. It is noteworthy that for χ = 0, the usual TPT1 expression for linear chains is recovered. While χ = 1.3827 is an exact result for a hard sphere system of triangles (see Appendix), in order to employ a similar approach for other geometrical realizations of planar rings, we postulate the following general form for the contribution to ring formation for the SAFT-VR Mie EoS where χ is defined as a parameter, function of m s and the actual geometrical connection of the ring. For chain fluids of m s > 3, and soft potentials TPT2 is not fully developed and analytical expressions are not available, hence, the procedure outlined above for hard triangles is not applicable. For the same reason, for soft potentials, the packing fraction is ill-defined and reverts essentially to a measure of the system density. We seek to exploit the correspondence between the SAFT-VR Mie EoS Table 1. Ring Molecule Geometries and Values for the Parameter χ (Eq 5) a a All molecular models are planar, built from m s tangent Mie spheres, bonded rigidly at a distance σ from the center of adjacent beads. The lines between the centers of neighbouring beads form equilateral triangles. , fixed geometry and molecular parameters: m s , ε, σ, λ r , λ a ) and use them as pseudoexperimental data to calculate χ from the EoS. We explored the phase equilibrium behavior of planar pure ring fluids for the cases of m s = 3, m s = 4, m s = 5 (two configurations), and m s = 7 (two configurations). While we use a single set of molecular parameters: ε/k B = 250 K, σ = 3.0 Å, λ r = 11.0 and λ a = 6.0 for the simulations, we will later extend these results employing a corresponding states principle. Simulation details and numerical values of the simulations are given as part of the Supporting Information. In Table 1, we summarize the resulting values of χ along with the graphical representation of the ring geometries used. Figures 1−4 illustrate the phase equilibrium (ρ−T projection) and vapor pressure diagram (T−P projection) of these models. In Figures 1−4, we include the calculation from the SAFT-VR Mie EoS, molecular dynamics results for rings, and also include the theory results for a chain (χ = 0) fluid formed by the same number of beads, m s . These figures quantify the importance of considering ring geometries separately; for a fixed value of m s , eq 5 is sensitive to the geometric connectivity of the ring. Figures 3 and 4 (m s = 5 and 7) are particularly revealing, as they show how even among ring molecules with the same value of m s , their particular conformation is relevant. In general, we see that for a fixed value of m s , the ring compounds display higher critical coordinates (i.e., temperature, pressure, and density) than the corresponding linear versions, in accordance with what is expected for common fluids, that is, compare the critical temperature of cyclooctane (647.2 K) to that n-octane (568.8 K). 28 Corresponding States Parametrization. In order to use the proposed models for specific fluids, it is necessary to define the values of the molecular parameters χ, m s , ε, σ, λ r , λ a . This can be done by using the criteria of the critical conditions of a   Table 1; (•) MD results for geometry (c) from Table 1; () SAFT-VR Mie EoS calculations with χ = 4.7188 (geometry (d) from Table 1); ( · ) SAFT-VR Mie EoS calculations with χ = 3.2222 (geometry (c) from Table 1 Langmuir XXXX, XXX, XXX−XXX pure fluid together with the evaluation of the liquid density at a specific condition 29 or alternatively by fitting to experimental data, for example, vapor pressure and liquid density. 10 Such an approach is exemplified in the Supporting Information by fitting the model of m s = 3 (triangles) for three discotic-like molecules: benzene, hexafluorobenzene, and toluene. Alternatively, the parameters can be found by using a corresponding state parametrization. The procedure is described in detail in ref 30 and used therein to provide parameters for thousands of molecular fluids employing the linear chain approximation. 31 It will be briefly detailed here for completeness. From the set of parameters required to describe a pure fluid, χ and m s are predetermined according to the shape and geometric connection of molecule (see Table 1). A further reduction in the number of parameters can be made if one recognizes that, for fluid phase properties, there is an intimate relationship between the exponents λ r and λ a of the Mie potential, that is, their specification is not independent of each other. 32 The observation is based on the premise that the van der Waals constant, α, which captures the first order contribution to the mean field dispersion energy of the Mie potentials, establishes a unique dependence between pairs (λ r , λ a ) of exponents, that is, fluids described by the same value of α exhibit conformal fluid phase diagrams. In that sense, for most common applications and without lack of generality, one can fix the attractive exponent λ a to 6.0. This choice will be made here and in what follows in the manuscript λ = λ r is the only one of the Mie exponents to be fitted. The van der Waals constant is then related directly to the repulsive exponent as For each of the molecular geometries (i.e., for a given set of values of m s and χ from Table 1) it is possible to correlate λ with the slope of the vapor pressure curve. This is done by taking a unique point of the saturation curve at a reduced temperature of T r = T/T c = 0.7. Such point is used to define Pitzer's acentric factor, ω. Tables for the acentric factor of pure substances are commonly available, otherwise the value can be found by its definition. 33 The core of our procedure is to recognize that given the fact that we can use the SAFT-VR Mie EoS to relate macroscopic thermophysical properties with potential parameters, we can explicitly evaluate the relationship between these two conjugate properties (one microscopic and the other macroscopic) and fit the results to a Padéseries From the knowledge of the corresponding λ, the van der Waals constant, α, may be found from eq 7. For each possible value of α, the principle of corresponding states states that the critical properties and phase behavior can be expressed uniquely for each fluid in reduced terms, that is, scaling with respect to the appropriate energy (ε) and length (σ) scales. Thus, we express in reduced units the temperature T* = k B T/ε, as well as the density ρ* = ρσ 3 . This unique relationship may also be expressed in parametric form if one does not desire to resolve the parent EoS. We propose the following expressions for the critical temperature T c *, and the reduced saturated liquid density at T r = 0.7, ρ*| Tr=0.7 Finally, to match the reduced behavior to that of a real fluid, we need then only two experimental points for which the corresponding values for ε and σ for a given fluid may be obtained by simple scaling  In summary, the model relies on three molecular parameters ε, σ, and λ, which require fitting to experimental conditions. This is done by matching three chosen thermodynamic state points: the acentric factor, ω (essentially a point on the vapor pressure curve), the critical temperature, T c , and the saturated liquid density, ρ| Tr=0.7 , at the reduced temperature T r = T/T c = 0.7, guaranteeing a good fit of the range of the potential, the energy and the length scales. In Table 2 we present the values of the parameters a i , b i , ..., f i required for eqs 8−10) for each of the ring geometries discussed in Table 1 along with their range of applicability.
In the following section we illustrate the application to some representative cases, where SAFT-VR Mie is compared to MD simulations and the available reference 28,34,35 or experimental 36 data.
Application to Pure Fluids. In order to test the applicability of the proposed expression for a RING , its parametrization based on corresponding state principia and the transferability of its molecular parameters from SAFT EoS model to MD simulations, we have selected some archetypal cases of ring fluids represented by m s = 3, 4, 5, and 7. Table 3 summarizes the pure ring fluids tested in this work together with the values for the critical temperature (T c ), the acentric factor (ω), the liquid density at T r = 0.7 (ρ Tr=0.7 ) used to obtain  The molecular parameters listed in Table 3 have been simultaneously used in both the SAFT EoS model and MD simulations to carry out phase equilibria and interfacial properties calculations. Selected numerical values obtained from MD simulations presented in Table 3 are included in the Supporting Information. The corresponding statistical deviations of the results presented here and those obtained from other authors have been summarized in Table S19 (Supporting  Information). These values are obtained from the comparison between molecular simulation and the reference or experimental data. The quoted figures also include the reference or experimental data reported in databases such as NIST-REFPROP 28 and DECHEMA, 36 and reported results are from molecular simulations that employ other molecular models (e.g., CG or united atom force fields). The SAFT EoS calculations and MD simulation agree with each other to within a deviation of 1% for phase equilibria predictions. Table S19 also includes the error observed when comparing the interfacial tension between the liquid and vapor bulk phases calculated from MD by employing the virial route against the experimental information. Interfacial tensions are not an input to the fitting of the potential or the EoS and provide a sensitive gauge to the quality and consistency of the density predictions obtained from the force fields.  Langmuir XXXX, XXX, XXX−XXX While the tables are too lengthy to discuss in detail, we selected two representative triangle fluids and compare the results graphically. Figure 5 and 6 show the comparisons for benzene, and tetrahydrofuran. Benzene and tetrahydrofuran are particularly relevant, as both have been previously modeled both as fully atomistic models, united atom representations where the hydrogens are lumped into the descriptions of the heavy atoms and as CG three-bead (equilateral triangle configuration) employing a related top-down approach through the use of fitting via the SAFT-VR Mie. 16,41 The high-fidelity models 23,39−41 are computationally much more expensive and do not provide any further level of accuracy, 23,39−41 exception being the fully atomistic TraPPE-UA-EH 40 model of benzene.
For the case of previous models of ring fluids fitted using the SAFT top-down methodology, 16,41 as these did not employ the appropriate form of the EoS, the molecular parameters obtained by fitting the EoS to experimental data had to be further refined (scaled) in order for the simulations to match experimental data, limiting the predictive nature of the procedure.
For the case of naphtalene and anthracene, the CG for models for rings display low ADDs in density and vapor pressure (data in the Supporting Information) . For the other fluids presented in Tables 3, this work represents predictions of the hypothetical phase equilibria, as the temperatures will most likely be above the decomposition temperatures for these organic fluids. Example is given in Figure 7 with the presumed phase equilibria of coronene.
Besides the phase equilibrium, Table S19 includes a comparison of the accuracy of the calculation of the interfacial tension, obtained through the Square Gradient Theory (details in the Supporting Information), showcasing the accuracy and predictive nature of the model. Figure 8 displays the density profiles across a vapor−liquid interface, ρ(z), profiles computed both from the theory and MD simulations at three temperatures for benzene. The classical hyperbolic profiles described by the theory agree with the MD results, and as expected, the interfacial width increases with temperature.
Applications to Mixtures. The extension of the model for the case of fluid mixtures is straightforward, as the eq 5 reverts to the case of a linear chain when χ = 0. For mixtures (which may include linear chains, rings or both) it is necessary to replace the original contribution of a CHAIN for the following expression, which accounts for both chain and ring fluids: where the subscript i runs over all the individual n c components (not beads) in the mixture. As the other terms in the perturbation expansion (eq 1) remain unchanged, the methodology to solve phase equilibria and obtain thermophysical properties remain unchanged with respect to the original versions of the EoS. In particular, for a mixture, one must define combinations rules to be used both in the simulations and within the theory. We employ the same choice originally suggested by Lafitte et al., 10 There is no expectation that the mixture behavior will match these idealized results, hence, a binary interaction parameter, k ij , is included in the cross energy term. This interaction parameter can be obtained from experimental data.
To illustrate the performance of this extension for fluid mixtures, we consider the binary mixture composed by hexane and tetrahydro-2H-pyran (THP). For this application, hexane is modeled as chain fluid formed by two tangential spheres, m s = 2 with parameters taken from the literature, 43 σ = 4.508 Å, ε/k B = 376.35 K, λ = 19.57, whereas THP is modeled as an equilateral triangle with m s = 3 with parameters taken from Table 3. Figure 9a displays the vapor−liquid phase equilibria (VLE) at 94 kPa, while Figure 9b displays the phase equilibria in a temperature−mass density diagram. The theory is resolved using a value of k ij = −0.022, which is obtained from the fitting of the EoS model to the experimental data. 44 Details of the calculation of the equilibria are given in the Supporting Information. The SAFT model correctly accounts for the zeotropic behavior of this mixture in the whole mole fraction range with a low absolute average deviation in boiling temperature (0.05%) and vapor mole fraction (0.91%.). While arguably such a quality of fit could be obtained by employing other versions of SAFT (and indeed other EoS), Figure 9 also includes the results obtained from performing molecular dynamics simulations using the same molecular parameters employed in the theory. The agreement between the theory and the simulations is unique to this treatment. Simulation details and numerical results are provided in the Supporting Information (see Tables S14 and S15).
A further test is the prediction of the homogeneous liquid density of the mixture at isothermal−isobaric conditions. Figure  10 displays the mass liquid density as a function of the mole fraction of the mixture at 298.15 K and 101.325 kPa (numerical data are presented in Table S16 of the Supporting Information) including also the SAFT calculations with k ij = −0.022 and MD

■ CONCLUSIONS
We present a modification to the SAFT-VR Mie EoS that allows the accurate representation of the properties of planar rings. Specific expressions for use in the case of rings composed of 3, 4, 5 (two configurations) and 7 (two configurations) beads are given in terms of a unique fixed parameter χ. Further expressions are provided that allow the parametrization of the EoS invoking the corresponding states principle. This approach provides molecular parameters from the knowledge of three state points: the critical temperature, acentric factor, and a reference liquid density of the pure fluid. The proposed approach has been tested by predicting the phase equilibria for discotic molecules and mixture of discotic and linear molecules and, when available, compared to experimental data. An important aspect of the proposed model is the direct link that exists between the molecular parameters underlying the Hamiltonian that defines the theory and the molecular simulations that employ the same parameters. This correspondence provides evidence that the approximations made within the theory are robust and accurate. The parallelism between the theory and the simulation models becomes important when dealing with interfacial and transport properties, where the theoretical developments are less refined. As an example, in Table S15 of the Supporting Information we present the results of the interfacial tension for a heptane−THP mixture where, as expected, the tension decreases as temperature and/or mole fraction increases, and while no isobaric experimental data is available to compare to, we have enough confidence in the models to trust these results. 45 In essence, the procedure of fitting molecular EoS parameters to experimental data becomes a shortcut to parametrize force fields for molecular simulations, with the advantage of producing effective pairwise intermolecular potentials which are fitted not to a few properties, but to provide the best compromise fit of the full free energy landscape. Comparison of structural and transport properties, such as radial distribution functions, diffusion, shear viscosity, and thermal conductivity are beyond the scope of this manuscript, although there are indications 16,41 that suggest a similar level of accuracy.

■ APPENDIX
Following the second order thermodynamic theory (TPT2) for hard spheres as proposed initially by Wertheim, 21 one can express the compressibility factor Z = P/ρk B T of a hard sphere chain of m s spherical segments by where Z ref is the reference term for hard spheres, Z 1 and Z 2 correspond to the first and second-order terms in the perturbation series, respectively, P, T, and ρ correspond to the pressure, temperature and molecular number density, respectively. In order to produce a closed expression for the reference and first-order terms of the Z expansion, one can  Table S15 of Supporting Information).   In eq A.5, g HS (σ) is the contact value of the radial distribution function for hard spheres, and g HS (3) , is the triplet correlation function of three touching hard spheres. Muller and Gubbins 47 evaluated this latter quantity and provided 21 a closed form expression for g HS (3) , which for the case of an equilateral triangle, θ = 60°, becomes