Accurate statistical associating fluid theory for chain molecules formed from Mie segments

Accurate statistical associating fluid theory for chain molecules formed from Mie segments Thomas Lafitte,1 Anastasia Apostolakou,1 Carlos Avendaño,1,2 Amparo Galindo,1 Claire S. Adjiman,1 Erich A. Müller,1 and George Jackson1,a) 1Department of Chemical Engineering, Centre for Process Systems Engineering, Imperial College London, South Kensington Campus, London SW7 2AZ, United Kingdom 2School of Chemical Engineering and Analytical Sciences, The University of Manchester, Sackville Street, Manchester M13 9PL, United Kingdom


I. INTRODUCTORY BACKGROUND
The goal of this contribution is the development of an accurate equation of state (EOS) for fluids based on a generic and versatile force field that captures the key features of the interactions for a broad range of real substances. The gradual unmasking of the role played by intermolecular forces in determining the thermophysical properties of fluids and the marked sensitivity of the observed behaviour on the nature and range of the interaction (the softness/hardness of the repulsions or the extent and range of the electrostatic attractions) has a long and illustrious history. 1 In order to place our generic approach in the appropriate context we first summarize the key findings, as the subtleties are not always fully appreciated or have even been forgotten altogether. This is a) Author to whom the correspondence should be addressed. Electronic mail: g.jackson@imperial.ac.uk particularly true in the development of EOSs where the precise form of the underlying force field is often deemed as secondary to the theoretical formulation. A detailed knowledge of the various contributions to the intermolecular potential is central to the accurate description of molecular systems; a robust link between the force field and the thermodynamic properties is provided in our paper. Much of the studies on "corpuscular" interactions in the late seventeenth and eighteenth centuries centered around the ideas of Newton, who had already considered attractive forces between particles that varied as an (undetermined) inverse power of their separation in his Principia Mathematica, 2 but these were largely speculative with minimal empirical validation. The existence and nature of repulsive forces was even more disputed and poorly understood, though Newton and his supporters Rowning and Bošković recognized that repulsions must also be active at short distances to account for the incompressibility of dense fluids and materials; some of the sketches of the interactions between "Boscovichian" particles 3 look strikingly similar to the potential functions that we are familiar with today. Perhaps the exception in terms of firm experimental confirmation was the discovery by Priestley, 4 Cavendish, 5 and Coulomb 6 that charged bodies (no matter what their physical dimension) interact through an inverse-square force law, i.e., a potential function of the gravitational form which involves the reciprocal of the separation, 1/r. Young 7 and Laplace 8 were convinced at the turn of the nineteenth century that one had to take into account both repulsive and attractive interactions between particles to explain the observed cohesion and capillarity in fluids, but it took much longer to establish the precise forms of these interactions.
In trying to interpret the attractive potential between particles in terms of electrostatic forces, Mosotti 9 was one of the first to develop an explicit potential for the interactions in dense electrically neutral fluids as an exponentially damped version of the Coulombic interaction (i.e., u a (r) ∝ exp (−λ/r)/r with λ as a system-dependent potential-range parameter) based on Poisson-Laplace electrostatics; this form of force field has since made an appearance in a number of physical scenarios and is now commonly referred to as the Yukawa 10 potential. 11 By 1857 Clausius 12 fully recognized and explicitly stated that molecules repel each other at short relative separations and attract each other over a moderate range, starting a passionate quest amongst his contemporaries for the specific mathematical form of the overall interaction potential that best described the transport properties of dilute gases. In his seminal paper on the dynamical theory of gases, Maxwell 13 extended his earlier work on perfectly elastic hard spheres (which borrowed heavily from the kinetic theory of Clausius) to particles interacting through a pairwise repulsive force taken to be inversely proportional to a power of the separation between their centres of mass (i.e., the repulsive analogue of the Newtonian attractive form). Maxwell 13 noted that the treatment of molecules as spherical hard cores leads to an inadequate description of the diffusion of gases. By examining his experimental data for the viscosity of air at different temperatures (where he observed a direct proportionality with temperature, which is now known to be erroneous), Maxwell 13 was able to deduce that, according to his theory, the exponent characterizing the repulsive force had to take the value of five to enable one to reproduce the observations, i.e., f r (r) ∝ 1/r 5 ; the repulsive pair interaction potential between these Maxwellian particles would thus be of the form u r (r) ∝ 1/r 4 . Though in modern terms this may appear to be an overly soft interaction, Maxwell pointed out that he is representing the gas as "small bodies or groups of molecules repelling one another"; at the end of our paper we will return to the interesting consequences of this type of coarse graining of the interactions between collections of molecules in terms of single spherical cores. In work along the same vein, Boltzmann 14 came to the contrasting conclusion that a completely different intermolecular force law (including a generic attractive contribution) could be used to reproduce the experimental gas viscosity.
The emergence of these ideas enabled van der Waals 15 to develop his equation of state for fluids which, for the first time, allowed one to represent vapour-liquid equilibria (VLE).
Within the van der Waalsian model, molecules are represented as perfectly hard-spherical cores (which maintain the low compressibility of the dense fluid) with additional longrange attractive interactions (which stabilize the condensed phase); however, the specific form of the attractive forces remained unspecified, the contribution of the interactions to the thermodynamic properties of the fluid being described instead with an integrated energy which conceals the precise algebraic dependence of the interaction on the molecular separation. Van der Waals toyed with the idea of a potential of the Yukawa form but did not pursue it further. In a series of contributions to the Philosophical Magazine spanning the last two decades of the nineteenth century, Sutherland [16][17][18] followed van der Waals' lead and coupled a hard-sphere model with pair attractive forces which he postulated would vary as the inverse fourth power of the separation (i.e., a pair potential of the form u a (r) ∝ 1/r 3 , obtained rather mysteriously as the second term of an expansion of the gravitational force which also depended on the particle's mass). With this inverse-fourthpower law for the attractions, Sutherland [16][17][18] was able in large part to reconcile the existing experimental data for the thermophysical properties such as the equation of state, 16 the vapour-liquid interfacial tension, 17 and the viscosity of gases. 18 A significant leap in unveiling the empirical form of the repulsive interactions between particles was made by Mie 19 in 1903, almost half a century after the pioneering works of Clausius 12 and Maxwell. 13 Like Maxwell, Mie represented the repulsive pair potential between molecules in terms of a generic power law with a term of the form u r (r) ∝ 1/r λ r , where λ r is the exponent characterizing the nature of the repulsions, but he now also included the attractive interactions implicitly with an integrated van der Waals energy. Mie's treatment of the potential is generic in the sense that it can be used to describe explicitly repulsive interactions of varying softness/hardness, and attractive interactions of different range as described through the magnitude of the van der Waals constant. We will make use of this particular versatility of the Mie potential form to represent the interactions between molecular segments later in our paper. By examining the equation of state and thermal properties of metals (namely, mercury, copper, nickel, and zinc), Mie showed that in order to provide a consistent description of the experimental values of the compressibility, a steeper repulsive exponent was required (in his case λ r ∼ 5 instead of the Maxwellian value of λ r = 4), something that Boltzmann had also become increasingly convinced about in his later years. The key finding that thermodynamic derivative properties such as the compressibility are particularly sensitive to the steepness of the repulsive pair potential is a principal theme running through our current work. It is not widely recognized, however, that Grüneisen 20,21 appears to have been the first to write down an explicit expression for the overall interaction potential between pairs of particles (again for metallic systems) as the sum of a generic Maxwellian repulsive contribution u r (r) ∝ 1/r λ r and a generic Newtonian attractive contribution u a (r) ∝ 1/r λ a , where λ a characterizes the variable range of the attractive interaction; this form of interaction is often ascribed to Mie or, even more-commonly referred to as the generalized Lennard-Jones potential. In his work Keesom 22,23 focussed on the form of the attractive potential by examining the second virial coefficient of molecules modeled as rigid hard spheres with Newtonian attractive interactions; this type of force field is now generally referred to as the model of Sutherland in recognition of the Australian's seminal contribution. Keesom 22,23 found that the description of the temperature dependence of the second virial coefficient (another thermodynamic derivative property) of argon was relatively insensitive to the value of the attractive exponent (at least in the range λ a = 3 to 5), with a marginal preference for λ a = 4 as the optimal value.
Lennard-Jones [24][25][26] is now widely credited as having explicitly expressed the full pair potential in the form u(r) = C r /r λ r − C a /r λ a (C r and C a representing the corresponding repulsive and attractive constants) which, as has already been mentioned, is the generic expression that had been reported by Grüneisen twelve years earlier. In his paper of 1924 on the temperature dependence of the viscosity of gases, Lennard-Jones (Jones at the time) 24 assumed an interaction potential with a variable λ r and λ a = 2, making this choice for the attractive exponent to render the triple integration for the collision dynamics tractable; he was able to describe the experimental viscosity of argon gas over a broad temperature range for repulsive exponents spanning λ r = 6 2 3 to 20. The representation of the second virial coefficient of argon was the focus of a subsequent well-known paper, 25 in which Lennard-Jones used the more-generic form of the potential; he concluded that values of λ r = 13 1 3 and λ a = 4 provided the best description of the equation of state at low density, though again no unique set of exponents could really be favoured over another. The particle diameter that he obtained from the viscosity and virial coefficients of argon were found to be equivalent to within 10%. By the time that Lennard-Jones 26 delivered his lecture to the Physical Society in 1931, he had settled on the ubiquitous LJ potential form which bears his name with exponents of λ r = 12 and λ a = 6, the latter value being consistent with the findings of the budding discipline of quantum mechanics.
With the advent of quantum mechanics in the early part of the twentieth century, a more rigorous and detailed description of the electrostatic interactions between particles became possible. The repulsive interactions could now be seen as the result of the Pauli exclusion principle between electronic clouds (or the exchange interaction discovered by Heisenberg and Dirac), though the precise dependence on the interparticle separation was still difficult to describe algebraically; a compact exponential form (exp (−r/σ 0 ), with σ 0 representing a size parameter related to the Bohr radius) can be used for the repulsive potential resulting from the overlap of the electron density of s orbitals but only for the simplest atoms (cf. the papers of Slater et al. 27,28 ). On the other hand, Wang 29 had shown in 1927 that the attractive interactions (dispersion forces) between atoms are the result of fluctuating electric dipoles, which when treated at the quantum-mechanical level give rise to a contribution proportional to 1/r 6 ; this form of attractive potential was later confirmed by London and Eisenschitz, 30 and is now often referred to as the London dispersion interaction in recognition of London's accessible theoretical description. 31 The fluctuations of the electrons in larger molecules also leads to quadrupoles and higher multipoles which give rise to a sum of contributions in powers of λ a = 6, 8, 10, etc. 32,33 It is therefore apparent that in order to retain a simple yet rigorous quantum-mechanical form of the pair interaction between molecules one could combine an exponential repulsive dependence with a sum of attractive terms in 1/r λ a ; in the particular case of an attractive contribution of the London form this corresponds to the Buckingham exp-6 potential. 28,34,35 The reader is directed to the monograph by Stone 36 for an upto-date account of the different types of intermolecular multipolar interactions. Before we conclude our brief overview of the important quantum-mechanical insights into the development of classical force fields we should also mention that the square-well (SW) pair potential (or more precisely its quantum-mechanical analogue the "particle-in-the-box") is also rooted in the birth of the wave mechanics of Schrödinger, Fermi, Dirac, and others. The popularity of the square-well potential in the description of the interactions between classical particles stems from its very simple mathematical form, which allows for exact and tractable statistical mechanical calculations (e.g., the algebraic evaluation of the second and third virial coefficients 37,38 ).
By the 1950s the nature of the intermolecular pair potential was well understood, as was the sensitivity of the thermophysical properties of fluids to the precise form of the force field. 39 Molecular-beam experiments allowed for a direct determination of the pair interaction, and both the exponential form of the repulsions at short distances and the London form of the attractions at moderate distances were confirmed. 40 It was also well recognised that a universal form of pair potential (e.g., with fixed exponents characterizing the repulsive and attractive interactions) would not be generally appropriate to represent the intermolecular forces of different substances. 39,40 As Hirshfelder et al. 39 pointed out in 1954: "The only conclusion that can be drawn [from the body of existing work] is that a potential function more realistic than the Lennard-Jones (12-6) potential must be chosen if consistent information about the intermolecular forces is to be derived from the various physical properties." The fact that the description of the intermolecular interactions in fluids and materials with the Lennard-Jones (12-6) potential is still so prevalent in modern theoretical and computersimulation studies is rather surprising. The reason can in part be ascribed to the larger number of intermolecular parameters that need to be determined for potentials of a more-generic form. This is hampered by the lack of an accurate algebraic theory for the thermodynamic properties of fluids described with more-realistic models, which would facilitate the efficient development of the force fields for broad classes of systems.
In our paper we describe the fluid-phase equilibria and thermodynamic derivative properties of fluids of varying complexity by retaining a generic representation of the softness/hardness of the repulsive interactions and the range of the attractive interactions. The Mie (generalized Lennard-Jones) pair potential affords the appropriate versatility in this regard making it an ideal choice of force field; a power law of variable range can be used to provide a good representation of the more-realistic exponential form of the repulsive interactions, while still retaining mathematical tractability. A prerequisite to assessing the adequacy of the Mie potential in representing the interactions in real substances is an accurate description of the thermodynamics and structure of the fluid.

II. EQUATIONS OF STATE FOR COMPLEX FLUIDS
The development of an accurate molecular equation of state for the thermodynamic properties of fluids of associating chain molecules is an important goal as it provides a framework for the representation of the behaviour of a variety of real fluids and their mixtures with minimal computational cost. Moreover, a rigorous theory of fluid systems is invaluable in enabling an in-depth molecular-level understanding of the complex phenomena often encountered in engineering applications. The statistical associating fluid theory (SAFT), 41,42 which is based on the first-order thermodynamic perturbation theory (TPT1) of Wertheim, 43-48 is among the most successful of the modern algebraic descriptions of associating chain fluids. The literature on the development and use of the SAFT EOS is too extensive to be enumerated here; the reader is directed to excellent reviews which provide a comprehensive description of the theory and its successful application to the phase equilibria of complex fluids ranging from small associating compounds to multi-component systems including polymers, surfactants, micelles, liquid crystals, asphaltenes, and electrolyte solutions. [49][50][51][52][53][54][55] The main feature of the SAFT approach is that the thermodynamic and structural properties of a fluid of associating chain molecules are described in terms of those of a corresponding monomer (reference) system comprising the segments that form the molecule. One must first define explicitly the intermolecular potential of the reference fluid, typically an isotropic potential, and then evaluate the system's free energy and structure via the radial distribution function (RDF). There are now many different embodiments of the generic SAFT approach, differing mainly in the interaction potential used to model the reference fluid. The earlier versions of the SAFT EOS were based on a reference system of hard spheres. In morerecent versions of the theory, reference fluids represented with the square-well, [56][57][58][59][60] 56 and LJ 65 potentials) has been reported, [72][73][74] and it was demonstrated that EOSs based on a hard-core or LJ pair potential do not generally provide one with an adequate representation of thermodynamic second-derivative properties such as the isothermal compressibility or the speed of sound in the compressed liquid phase; the description of these properties is found to be no better than that obtained with EOSs of the conventional van der Waals form such as the Peng-Robinson 75 (PR) or Soave-Redlich-Kwong 76 (SRK) cubic relations. In order to overcome these limitations a generic description of the interactions can be introduced, 72 using the Mie potential to describe the monomer-monomer interactions between the segments making up the molecules. This potential can be expressed as 19-21, 25, 26 where r is the distance between the spherical segments, σ is the position at which the potential is zero (segment diameter), is the potential depth, and λ a and λ r are the attractive and repulsive exponents, respectively; the pre-factor C ensures that the minimum of the potential is at − regardless of the values of the exponents. In early developments within the SAFT-VR framework 56, 65 expressions were derived for the EOS and RDF in the particular case of Lennard-Jones chains (SAFT-VR LJC). To extend further the predictive capabilities of the SAFT-VR treatment for the more-general case of Mie potentials, a description of the RDF at contact of the monomer Mie reference fluid was then proposed; 72 here we will refer to the resulting EOS as SAFT-VR Mie 2006 to avoid any confusion with our current work. The SAFT-VR Mie 2006 EOS has been shown to provide a significantly improved description of the dense liquid phase, especially for the secondderivative properties, while still retaining an accurate representation of the vapour-liquid equilibria. 72,73 As a particular example of the marked improvement that can be achieved in the description of thermodynamic second-derivative properties by using a Mie interaction with a variable repulsive exponent to control the softness/hardness of the potential, the speed of sound (a stereotypical derivative property which is related to the compressibility and heat capacity of the fluid) of selected n-alkanes calculated with SAFT-VR Mie 2006 72 is compared in Table I to the corresponding description with SAFT-VR SW, SAFT-VR LJC, and PC-SAFT, as well as with the PR and SRK cubic EOSs. The previously reported values of the molecular parameters for each approach are adopted in this assessment; a more-complete comparison of the EOSs (also including soft-SAFT) with parameters estimated from thermodynamic derivative properties as well as VLE will be made in Sec. V B. The added flexibility of varying the range of the repulsive interactions is clearly the key to a reliable description of the derivative properties of the dense fluid, in addition to providing an accurate and versatile platform for the representation of the fluid-phase equilibria. The adequacy of the description of the density and thermodynamic derivative properties of model fluids with the SAFT-VR Mie 2006 EOS has been assessed by comparison with molecular-dynamics (MD) simulations of spherical molecules interacting via λ r -6 Mie potentials with λ r = 8, 10, and 12; the theory is found to provide a good description of the simulated thermodynamic properties of Mie fluids, with the greatest deviations seen for the heat capacity. 77 An assessment of the thermodynamic properties of model chain molecules formed from Mie segments with the SAFT-VR Mie EOS developed in 2006 72 has not been made until now, TABLE I. Speed of sound u of selected n-alkanes in the condensed-liquid phase. Comparison of experimental data with the description of equations of state including the Soave-Redlich-Kwong 76 (SRK), Peng-Robinson 75 (PR), PC-SAFT, 71 SAFT-VR SW, 56 SAFT-VR LJC, 65 and the SAFT-VR Mie approach developed in 2006. 72 The average absolute deviation AAD(%) is computed as AAD(%) = 1 and a detailed analysis will be presented later in our current paper. It is important to realise that the quality of the theoretical representation will depend not only on the accuracy of the approximations inherent in the Wertheim TPT1 description of the chain contribution but also on that of the underlying theory for the EOS and contact value of the RDF of the reference monomeric Mie fluid. Such a comparison will therefore represent a stringent test of the approximations used at each level of the development of our new SAFT-VR Mie EOS. An accurate algebraic expression for the EOS of simple monomeric fluids can be developed by using an empirical or semi-empirical scheme based on exact computer-simulation data. This has proven to be a very successful approach in the case of the LJ fluid as exemplified by the work of Johnson et al. 78 or Kolafa and Nezbeda 79 who made use of multiparametric relations to correlate the simulation data over a wide range of densities and temperatures. As far as chains formed from LJ segments are concerned, a marked sensitivity of the Wertheim TPT1 description to the precise representation used for the RDF of the reference momoner LJ fluid has been demonstrated. 80 The EOS of Johnson et al. 78 has been coupled with an accurate empirical description of simulation data for the contact value of the RDF of LJ fluids to develop a high-fidelity EOS for LJ chains; 63 this EOS for the LJC system is the precursor to the subsequent SAFT approaches for associating LJ and LJC fluids 81,82 including the soft-SAFT approach. 66,67 Despite the accuracy of the resulting description with these approaches, it should be pointed out that semiempirical expressions for the EOS and RDF of this type are restricted to a specific intermolecular potential model (the LJ (12-6) potential in this particular case), and the applicability of the methodology is, of course, limited to the temperature and density ranges of the computer-simulation data used in the parameterization.
A more-promising route to the description of the equation of state and structural properties of a generic reference fluid is through a rigorous statistical mechanical theory. Among these methodologies, thermodynamic perturbation theories are more suitable than integral-equation approaches for the development of an algebraic formulation. The perturbation theories of Weeks, Chandler, and Andersen (WCA), 83 and of Barker and Henderson (BH) 84,85 are the most-ubiquitous approaches in this context. It has been shown, 80,86 however, that the WCA recipe does not provide the most-accurate description of the properties of chain molecules within a SAFT-like treatment. In our current paper we will show that the perturbation approach used to develop the SAFT-VR Mie EOS in 2006 72 is also somewhat inaccurate in its description of the vapour-liquid equilibrium of long LJ chains. Here, we propose a new formulation for both the EOS and RDF of the monomeric Mie fluid. The novel approach is based on a rigorous application of the Barker and Henderson 84, 85 hightemperature perturbation expansion of the Helmholtz free energy of the Mie fluid up to third order. In addition, a significantly improved expression for the RDF of Mie monomers at contact is presented, following a self-consistent description of the pressure from the Clausius virial theorem and the density derivative of the free energy. In the special case of LJ fluids, our new generic SAFT-VR Mie EOS will be seen to provide a description of similar fidelity to that of the empirical EOS of Johnson et al. 78 (cf. Sec. V A 1).
In order to assess the adequacy of our new SAFT-VR Mie EOS, the calculated thermodynamic properties are first compared with the corresponding computer-simulation data for fluids of spherical Mie particles, for chains comprising up to 100 LJ segments, and for chains of Mie λ r -6 segments. The association contribution of the SAFT-VR Mie EOS is then also revisited in order to capture the effect of placing the bonding sites at different distances from the centre of the reference core. This turns out to be a useful additional feature of the new theory in the development of more-realistic models of the hydrogen-bonding interaction. A generic analytical expression of the association contribution to the free energy for an arbitrary geometry of the bonding sites is presented. Its accuracy is evaluated by comparison with the available molecular-simulation data for selected thermodynamic properties. As a final demonstration of the enhanced capabilities that our new SAFT-VR Mie EOS has to offer, the theory is used to describe the vapour-liquid equilibria and secondderivative properties for representative real substances where we demonstrate the importance and versatility of using a variable repulsive range with Mie potentials of the λ r -6 form. The accurate theoretical description allows for the size and energy intermolecular parameters and the repulsive exponents to be estimated with confidence from the bulk thermodynamic data of the different systems.

III. BARKER AND HENDERSON PERTURBATION THEORY OF MIE FLUIDS
In this section we derive the Helmholtz free energy and radial distribution function at contact of Mie fluids using the BH perturbation theory for soft-core potentials. Different approximations for the EOS and RDF are available and the corresponding description is compared with that obtained from our new methodology. A comparison is also made with computer-simulation data for a wide range of values of both the repulsive and attractive exponents λ r and λ a of the Mie potential. Based on this analysis, we propose simple and compact expressions for the EOS and RDF by making use of the mean-value theorem (MVT). The detailed implementation of these expressions in the development of the SAFT-VR Mie EOS for associating chains of Mie segments is described in Sec. IV.

New SAFT-VR Mie formulation
The starting point of the BH thermodynamic perturbation theory is to decompose the Mie potential given in Eq. (1) into a sum of a reference u 0 (r) and a perturbation u 1 (r) term: 84 with Barker and Henderson 84 showed that it is then possible to write the residual Helmholtz free energy a res = A/(N s kT) as a res = ∞ n=0 (β) n a n , where a 0 is the Helmholtz free energy of the reference fluid, a 1 , a 2 , . . . , a n are the first-order, second-order, and nth-order perturbation terms, respectively, β = 1/(kT) (with k representing the Boltzmann constant, and T the absolute temperature), and N s is the number of particles (spherical segments). The reference (repulsive) system defined by Eq. (4) is not convenient in practice because its Helmholtz free energy a 0 and RDF g 0 (r) are not generally known. However, this reference system can be mapped onto a system of effective hard spheres, whose diameter d remains to be determined. 85 In order to evaluate this diameter, an expansion of the Helmholtz free energy has to be carried out with respect to two coupling parameters (α and γ ) which vary the steepness and depth of the potential in such a way that α = γ = 1 corresponds to the original Mie potential u Mie (r) and α = γ = 0 to the effective hard-sphere potential u HS (r); the reference potential u 0 (r) is obtained in the limit of γ = 0 and α = 1. Barker and Henderson 84, 85 demonstrated that one can nullify the firstorder term of the Helmholtz free energy in the α-expansion about a hard-sphere system of diameter d when All the thermodynamic functions of the reference system as well as its RDF can therefore be approximated by those of this effective hard-sphere system. For instance, a 0 and g 0 (r) may be obtained from where a HS and g H S d (r) are the free energy and radial distribution function of the hard-sphere fluid of effective diameter d. The integral involved in the evaluation of d can be computed using an appropriate numerical integration method, e.g., one based on Gaussian quadrature. 87 It follows that the Helmholtz free energy of a soft-core potential can be written as a series of the form The reference hard-sphere free energy a HS is calculated by integrating the well-known Carnahan and Starling 88 EOS, a H S = P H S /(ρ s kT ) − 1 /ρ s dρ s , where η = ρ s π d 3 /6 is the packing fraction of the reference hard-sphere system, and ρ s = N s /V is the segment density of the Mie fluid (with V the volume of the system).
In our current work, we propose a truncation of the series at third order instead of the more-widely used second-order expansion. This is motivated by the fact that for large values of the repulsive and/or attractive exponents, the vapour-liquid critical point can be located below T * = kT/ = 1; in such cases truncating the high-temperature expansion at second order generally results in a significant overshoot of the critical temperature. 89 It will be shown in Sec. V A that the incorporation of a third-order term considerably improves the description of the VLE in the vicinity of the critical region.
The first-order perturbation term a 1 (mean-attractive energy) is calculated from The key quantity is the reference hard-sphere RDF g H S d (r), which we obtain by solving the Ornstein-Zernike integral equation using the reference hypernetted chain (RHNC) closure and the Malijevsky and Labik 90 formula for the hard-sphere bridge function. A comparison of the meanattractive energy obtained in this way with the corresponding molecular-simulation data is shown in Fig. 1. At this stage, we emphasize that Eq. (12) is slightly different to the approximations made in the previous work on the SAFT-VR Mie EOS of 2006 72 where Note that though the latter approach is not a formally correct application of the BH theory, it is a useful approximation to derive a simple and compact algebraic expression for a 1 when the SAFT-VR mapping 56 is used. The corresponding description for a 1 using Eq. (14) with the Malijevsky and Labik 90 RHNC recipe for g H S d (r) is also shown in Fig. 1(a). It can be seen that the rigorous application of the BH theory leads to a more-accurate representation of a 1 . In order to develop an EOS for use in engineering applications an algebraic expression for a 1 is highly advantageous; such an expression will be presented in Secs. III A 2 and III A 3. The approximations inherent in Eqs. (13) and (14) have also been employed recently by Betancourt-Cardenas et al. 91 to derive the perturbation terms in the Helmholtz energy expansion of the LJ fluid; this explains in part why these authors find discrepancies between their simulation data for the perturbation terms and the original values of Smith et al. 92 An exact evaluation of the second-order term a 2 (fluctuation contribution) is much more difficult as it requires a knowledge of two-, three-, and four-body correlation functions of the reference system. 92 When the superposition approximation is used to evaluate the three-and four-body functions, the integrals may be evaluated numerically. Such calculations were performed by Smith et al. 92 showing significant improvement compared to the local compressibility approximation (LCA) 85 and macroscopic compressibility approximation (MCA) 85 that are commonly used to represent a 2 . Despite the good description obtained with this more-rigorous approach compared to the simpler LCA and MCA recipes, Smith et al. 92 showed that the agreement with molecular simulation remains only fair at higher densities. Moreover, the integrals involved in the calculation of a 2 using the superposition approximation cannot be calculated analytically. With the aim of deriving an algebraic equation of state, our objective is therefore to find a simpler, but accurate, expression for the fluctuation term a 2 . The approach that we advocate here is based on the work of Zhang 93 who proposed an improved MCA expression by assuming that the numbers of molecules in neighbouring coordination shells are correlated. Within the improved MCA approach, the term a 2 can be expressed as 93 where χ = 8.23η 2 is a simple correction pre-factor, and K HS = kT(∂ρ s /∂P HS ) T is the isothermal compressibility of the hard-sphere reference fluid, which is obtained from the density derivative of the Carnahan and Starling 88 expression for the compressibility factor Z = P HS /(ρ s kT): By introducing this simple correction, Zhang 93 showed that a significant improvement in the representation of a 2 is achieved in the case of the square-well potential of range λ = 1.5, especially for systems at intermediate to high densities. This correction fails however in the case of soft-core potentials. A recent modification was proposed by Paricaud 87 to address this deficiency by including an additional correction contribution proportional to η 5 . A similar approach is used in our work, but employing a more-generic correction factor: where x 0 = σ /d. The functions f i (i = 1, 2, 3) depend on both the attractive and repulsive exponents of the Mie potential through a single dimensionless van der Waals-like attractive constant α defined as (18) The specific functional dependence of f i on α is given in Eq. (20).
It is important to emphasize that while the fluctuation term displays a complex non-monotonic behaviour with density, the third-order term is known to exhibit a relatively simple dependence, at least for the square-well potential. This was shown by Espindola-Heredia et al. 94 who evaluated the perturbation terms up to fourth order by means of Monte Carlo simulation for fluids characterized by square-well potentials with ranges in the interval 1.1 ≤ λ ≤ 3; they then proposed an accurate equation of state using a set of empirical functions that simultaneously reproduced the exact values of the perturbation expansion up to fourth order and the virial coefficients up to third order. In view of the success of their approach, a similar strategy is chosen here for the Mie fluid with the use of the following empirical expression: where the expression is cast in a functional form so that it is independent of temperature. This approximation is chosen in order to reduce the number of adjustable empirical coefficients and is a reasonable assumption as the third-order term is found to be essentially insensitive to the temperature. There is no unique way of determining the functions f i . A rigorous approach would be to first determine a 2 and a 3 as functions of density and temperature for different combinations of λ r and λ a by means of molecular simulation. The functions f i could then be determined by correlating the simulation data as in the work of Espindola-Heredia et al. 94 The drawback of this methodology is the relatively large uncertainty in calculating the third-order term a 3 . Moreover, very small differences in the behaviour of a 3 with density are found to lead to significant differences in the location of the critical point. A different methodology is therefore employed using a combination of a knowledge of the accurate simulated values of the fluctuation term a 2 and selected simulated VLE and critical-point data for different Mie potentials. The main difficulty of such an approach lies in finding a closed functional form for f i which can be used to represent the selected data using the unique constant α. We have found a satisfactory de-scription with the following functional form: (20) Two steps are made in order to determine the coefficients φ i, n . First, we consider only the coefficients φ i,n (i = 1, . . . , 3) that are involved in the representation of a 2 , obtained by simultaneously correlating the exact values of a 2 and the VLE data, as determined by molecular simulation. For the exact values of a 2 , we have performed traditional Metropolis Monte Carlo simulations for reference Mie (λ r -λ a ) fluids with selected values of the repulsive exponent, keeping the attractive exponent fixed at its London value of λ a = 6 (8-6, 12-6, 14-6, 20-6, 30-6) for a fixed reduced temperature of T * = kT/ = 1; the simulation procedure is described in the review by Barker and Henderson. 85 The literature values of the VLE data 95,96 for the saturated liquid densities and vapour pressures for Mie fluids characterized with exponents (8-6), (10-6), , , , , , , , and (24-12) are considered. Once the coefficients φ i, n (i = 1, . . . , 3) are determined, the remaining terms φ i, n (i = 4, . . . , 6) of Eq. (20), corresponding to the third-order perturbation contribution, are estimated from the VLE data, including critical temperatures and critical densities. We should stress that the resulting a 3 term therefore corresponds to an effective perturbation contribution which accounts for higher-order terms that have been neglected in the expansion. The resulting values of the coefficients φ i, n are summarized in Table II.
The reliability of the semi-empirical expressions for a 2 can be assessed by comparing the fluctuation term proposed in Eq. (15), which makes use of the correction function χ defined in Eq. (17), to the traditional MCA approach for the Lennard-Jones fluid (cf. Fig. 1(b)). As can be seen, the corrected MCA expression leads to a significant improvement of the description of a 2 , especially at the higher densities where the correct trend is captured. Note that the traditional MCA is appropriate at low and intermediate densities but unreliable at high densities. In our comparison we also include calculations with the previous form of a 2 employed in the SAFT-VR Mie 2006 EOS, 72 which was based on the approximate perturbation potential described by Eq. (13) and the LCA approach.
As can be seen from Fig. 1(b) the LCA is too crude an approximation to accurately represent the second-order perturbation terms of Mie fluids. In Sec. V A 1 we will show that despite this the LCA expression together with the approximate formulation for a 1 given in Eq. (14) provides a surprisingly good description of the VLE in the specific case of the Lennard-Jones fluid. This is due to the cancellation of errors between the first-and second-order perturbation contributions.
The behaviour of the perturbation terms of the Helmholtz free energy of Mie fluids is further examined in Fig. 2 where the a 1 , a 2 , and a 3 terms calculated with our new methodology are compared with exact simulation data for different values of the repulsive exponent (λ r = 8, 12, 14, 20, and 30) keeping the attractive exponent fixed at λ a = 6. Metropolis Monte Carlo simulations of the corresponding BH reference fluids are carried out here to provide exact values for the three perturbation contributions (see Ref. 85 for details of the methodology). Excellent agreement with the simulation data is apparent for the first-order term for systems characterized by this broad range of repulsive exponents. Small systematic discrepancies can, however, be observed at high densities due to the slight inaccuracies resulting from the calculation of the RDF of the reference hard-sphere fluid with the RHNC integral equation. It is important to recall that the representation of the first-order mean-attractive term a 1 follows the rigorous application of the BH theory and is not optimized to simulation data. This negligible difference between the theory and the simulation data for a 1 does, however, have some repercussions on the evaluation of the fluctuation term. The second-order contribution a 2 is simultaneously correlated to exact values of the fluctuation term and to VLE data. Our methodology therefore results in some discrepancies between the computed values of a 2 and the simulation data at high densities (particularly in the case of the more-repulsive potentials) due to a compensation of errors with the first-order term in this region. Nevertheless, the agreement at low densities is very satisfactory, and the overall behaviour for a 2 is qualitatively correct for systems with varying repulsive exponents. We should also note that the qualitative trend observed for the fluctuation term of the Mie fluids is very similar to the one reported by Espindola-Heredia et al. 94 for the squarewell fluids. The novelty of our approach lies in the evaluation of the third-order term using VLE data, including the critical region, for different Mie fluids. One could argue that this constitutes a mere correlation of empirical functions to the coexistence curves of the various systems, but it is reassuring to observe that the resulting density dependence of a 3 (cf. Fig. 2(c)) follows the same general behaviour as the simulation data for the third-order term obtained by Espindola-Heredia et al. 94 for square-well fluids of variable range. Indeed, it appears that the magnitude of the a 3 term is largest near the critical density. The trends observed in Fig. 2(c) are consistent with the slow convergence of the high-temperature expansion in the near-critical region. Interestingly, the magnitude of a 3 decreases for increasing values of the repulsive exponent (corresponding to a decrease in the range of the interaction) which is somewhat counterintuitive. This behaviour is again in agreement with the findings of Espindola-Heredia et al. 94 We will show in Secs. V A 1 and V A 2 that the inclusion of a third-order term in the high-temperature perturbation expansion allows for a very good description of the critical region without using a specific cross-over treatment. The density ρ * s = N s σ 3 /V dependence of the (a) first a * 1 = a 1 / , (b) second a * 2 = a 2 / 2 , and (c) third a * 3 = a 3 / 3 Barker and Henderson perturbation terms for selected Mie , , , , and (30-6) potentials at a temperature T * = kT/ = 1. The continuous curves represent the calculation using our theoretical expressions, Eqs. (12), (15), and (19). The circles represent the exact values we have determined by Monte Carlo simulation.

Algebraic expressions for the first-and second-order perturbation terms
In order to avoid using numerical quadrature techniques to evaluate the integrals in Eqs. (12) and (15) with the reference RDF calculated with the RHNC integral equation, we propose a formal mapping of the integrals by making use of the mean-value theorem. The procedure is at the heart of the SAFT-VR methodology for general spherically symmetrical potentials with hard-repulsive cores. 56 The approach is somewhat more complicated to implement when a softcore reference fluid is considered as will now be explained in detail.
The first-order perturbation term a 1 in the BH approach for a Mie potential is obtained by re-expressing Eq. (12) as The computation of the first integral I 1A is facilitated by using the algebraic expressions of the first-order term of the hardcore Sutherland potential as follows: where x = r/d is the reduced centre-centre distance between two hard spheres of diameter d.
The expression represents the first-order term a S 1 (η; λ) of the Helmholtz free energy of a system of hard-core Sutherland particles of diameter d characterized by the range parameter λ. This contribution can be evaluated within the SAFT-VR MVT mapping, 56 which allows one to express a S 1 (η; λ) as a simple analytical function without loss of accuracy compared to using numerical quadrature. Details of the procedure are discussed in Sec. III A 3. With the availability of an algebraic form for the Sutherland first-order term a S 1 (η; λ), the integral I 1A can be evaluated analytically as In order to obtain an analytical expression for the second integral in Eq. (21), we approximate g H S d (r) with a linear relation obtained from a first-order Taylor expansion of the contact value of g H S d (r) and its first derivative with respect to the separation r. Nezbeda and Iglesias-Silva 97 have pointed out that such a description is expected to be reasonable for short ranges of integration. The hard-sphere RDF g H S d (r) can therefore be evaluated with the following expression: The second integral I 1B can then be written as with To simplify Eq. (26) we define and By then using the Carnahan and Starling expression for ) and the Percus-Yevick approximation for the derivative at contact, a simple analytical expression for I 1B can be obtained as where As a result the first-order term a 1 for the Mie fluid can be approximated with the following algebraic expression: The same procedure is used to develop an algebraic expression for the second-order term a 2 : The integrals I 2A and I 2B can be written in terms of the first-order perturbation terms a S 1 of the corresponding hard-core Sutherland potentials. It can be easily shown that Finally, we recall that the expression for the thirdorder perturbation contribution a 3 proposed in Eq. (19) does not require a knowledge of the RDF of the reference fluid. Its evaluation is straightforward using Eqs. (19) and (20).

Improved mean-value theorem mapping
for the first-order perturbation term of the Sutherland potential a s 1 As will have become apparent from Sec. III A 2, the development of an analytical expression for the Helmholtz free energy of a fluid of Mie spherical segments within the SAFT-VR approach 56 involves the evaluation of the firstorder terms a S 1 (η; λ) of systems of hard-core Sutherland potentials of range λ (cf. Eq. (23)), where In the SAFT-VR treatment, 56 the RDF for the hard-sphere reference g H S d (xd) is calculated numerically by making use of an integral-equation theory. In our current work the RHNC theory with the accurate hard-sphere bridge function of Malijevsky and Labik 90 is used as it is known to provide a hardsphere RDF which is in very close agreement with computersimulation data. The integral in Eq. (37) is then evaluated as a function of the exponent λ and packing fraction η. Once the range and density dependence of a S 1 (η; λ) are known numerically the MVT can be applied to determine an analytical form for the integral by expressing Eq. (37) as where ξ is a distance satisfying the MVT equality. As shown in the original SAFT-VR treatment 56 it is possible to map g H S d (ξ ) to the contact RDF g H S d (d) evaluated at an effective density η eff . By using the expression for the contact value of the hard-sphere RDF which is consistent with the Carnahan and Starling 88 EOS (cf. Eq. (30)), a compact relation for a S 1 (η; λ) can be obtained as 56 An algebraic expression for the integral in Eq. (37) can be obtained once the dependence of the effective packing fraction η eff on λ and η is established. The original parameterization 56 for η eff (η; λ) was derived for exponents of the Sutherland potential in the range 3 < λ ≤ 12. We now propose a new mapping in order to account for a broader range of exponents, particularly for the more-repulsive systems, i.e., 5 < λ ≤ 100: The resulting values of the first-order Sutherland term a s 1 (η; λ) obtained with Eq. (39) along with the numerical results calculated from Eq. (37) with the hard-sphere RDF from the RHNC integral-equation theory of Malijesvsky and Labik 90 are shown in Fig. 3. It can be seen that the algebraic expression provides an excellent representation of a S 1 (η, λ) over the entire range of λ for which the effective packing fraction has been parameterized.
The density ρ * s = N s σ 3 /V dependence of the first a S * 1 = a S 1 / Barker and Henderson perturbation term of the Sutherland fluid with varying attractive exponent λ. The circles represent the essentially exact values of numerical quadrature with the radial distribution function of the hard-sphere reference system obtained using the RHNC integral equation. The continuous curves correspond to the analytical correlation determined by employing the mean-value theorem, cf. Eqs. (39)-(41).

B. Contact value of the radial distribution function for Mie fluids
The calculation of the RDF of Mie fluids of variable ranges evaluated at contact is revisited and derived through a rigorous application of the perturbation theory of Barker and Henderson. 85 It is important to recall that an accurate representation of the contact value of the RDF of the reference monomer fluid is a prerequisite in a Wertheim TPT1 description of the thermodynamic properties of chain fluids formed from Mie segments.

New SAFT-VR Mie formulation for g Mie (σ )
Just as for the Helmholtz free energy, a high-temperature perturbation expansion can be used to represent the RDF for Mie fluids: Another possibility is to expand the logarithm of the RDF, ln g Mie (r). When the expansion is taken to second order the RDF can be expressed as 85 The latter representation is preferred as it allows for a faster convergence of the series for the contact value g Mie (σ ) than the high-temperature expression of Eq. (42). Another attractive aspect of expanding ln g Mie (r) is that the resulting expression for g Mie (r) will never give rise to negative values. This turns out to be very important in the Wertheim TPT1 implementation, since the logarithm of the RDF at contact is required for the calculation of the contribution to the Helmholtz free energy due to the formation of chains of Mie segments. One should note that a second-order expansion for the structure of the fluid is consistent with a third-order expansion for the free energy as the RDF essentially corresponds to a variation of the latter with respect to density. In order to calculate the zeroth-order term g 0 (r), the approximation of Barker and Henderson is used so that the RDF for the soft-repulsive potential u 0 defined in Eq. (4) is approximated as 85 where g H S d is the RDF of the reference system of hard spheres of diameter d defined in Eq. (7), which accounts for the softness of u 0 . For the contact value at r = σ , the zeroth-order term can therefore be approximated as It is important to recognize that even though σ is the contact distance of the Mie fluid, it is not the contact distance of the hard-sphere reference system; the contact distance of the HS fluid is d which is implicitly temperature dependent (cf. Eq. (7)). As σ > d the calculation of the term g H S d on the right-hand side of Eq. (45) therefore requires a knowledge of the RDF of the HS fluid at reduced distances x 0 = σ /d > 1. Various semi-empirical representations are available for g H S d (r). In our current work we use the expression proposed by Boublík 98 which is valid for 1 < x 0 < √ 2: where the density-dependent coefficients k i are given by 98 and This representation for g 0 (σ ), cf. Eqs. (45) and (46), differs from the approximation made in previous work 65,72 where the zeroth-order RDF was simply evaluated with the Carnahan and Starling 88 expression for the contact value at the packing fraction η, essentially corresponding to the approximation that d = σ . The values of our new zeroth-order term are always smaller than those of the previous formulation but both tend to the correct zero-density limit, lim ρ→0 g H S d (σ ) = 1. In order to determine the contact value of the first-order term g 1 (σ ) we follow the usual SAFT-VR recipe 56 employing a self-consistent representation for the pressure P from the Clausius virial theorem 85 for the virial relation, and for the thermodynamic derivative, where Z HS is the compressibility factor of the HS system of diameter d. When one approximates the reference potential u 0 (r) ≈ u HS (r) as that of hard spheres of diameter d, noticing that the exponential function, is a step function so that its derivative with respect to r is a δ function at r = d, it then follows that Eq. (51) can be rewritten as After replacing g Mie (r) by Eqs. (43) and (44), and expanding the exponential function one can write By comparing Eq. (55) with Eq. (52), compact expressions for the first-order and second-order terms of the RDF at the distance r = d are obtained: and At this stage it is useful to make some observations about these expressions. It would appear that by making use of this self-consistent method, the evaluation of the perturbation terms of the RDF at the distance r = σ is not possible. We note, however, that the approximations g 1 (σ ) ≈ g 1 (d) and g 2 (σ ) ≈ g 2 (d) yield a good representation of the full g Mie (σ ) as long as the zeroth-order contribution g 0 (σ ) ≈ g H S d (σ ) is calculated rigorously from BH theory at the contact distance r = σ with the Boublík expression. 98 The crude approximation that g Mie (σ ) ≈ g Mie (d) does not yield an accurate representation due to the rapid increase in g Mie (r) between d and σ . It is important to emphasize that when the approximated perturbation potential as described in Eq. (13) is used in our current development then one recovers the expression for g 1 (d) used in the previous SAFT-VR Mie 2006 EOS. 72 While the first-order term is easily evaluated with a knowledge of the hard-sphere RDF of the reference HS fluid, the second-order term requires an intricate calculation of g 1 (r) for the Mie potential over all distances. Nevertheless, Eq. (57) can be simplified further by making use of the expression for the configurational energy together with the MCA approximation. After inserting the high-temperature expansion of the RDF in the energy, a correspondence between the first-order term of the RDF and the second-order contribution to the Helmholtz free energy can be obtained as 99 If one then considers the MCA approximation, which relates a 2 to the RDF of the HS fluid, and then equates Eqs. (58) and (59), the corresponding MCA approximation for g 1 (r) for separations r > σ is obtained as Note that the MCA approximation for g 1 (r) at the contact distance σ is not useful since it leads to the incorrect value, g MCA 1 (σ ) = 0. By substituting Eq. (60) into Eq. (57) we obtain the MCA approximation for the second-order RDF of the Mie potential at the distance d which is now written in terms of the known RDF of the HS fluid: In order to validate the approximations made in our new approach, we combine the zeroth-, first-and second-order contributions given by Eqs. (46), (56), and (61), respectively, into the perturbation expression, Eq. (43), to evaluate the contact value of the RDF of the Mie fluid, assuming that g 1 (d) ≈ g 1 (σ ) and g 2 (d) ≈ g 2 (σ ). The resulting values of g LJ (σ ) for the LJ (12-6) fluid are shown in Fig. 4(a) as function of density at a reduced temperature T * = 1, making a comparison with the "exact" contact values of the RDF obtained from correlations to simulation data. 63 To gain insight into the influence of truncating the series expansion, we also display the resulting values obtained with the zeroth-and firstorder perturbation expansions. The introduction of a secondorder term in the perturbation is seen to enhance the accuracy of the RDF at low temperatures (cf. Fig. 4(b)), although the MCA approximation for g 2 still fails to capture the slope of the contact value of the RDF at zero density for T * = 1 (cf. Fig. 4(a)). The positive slope obtained for the temperature dependence of the contact value of LJ fluids is also seen for square-well fluids, as described theoretically by Tavares et al., 100 who demonstrated that for a SW interaction with a range λ ≥ 1.5, a positive slope is obtained for any temperature. It has been argued that it is not possible to capture the subtle temperature dependence of the contact value with high-temperature perturbation expansions. 100 It is, however, reasonable to assume that, if no approximations (such as the MCA) were made in calculating the second-and higher-order terms, then the inadequacies at low density and temperature would not be present.
Here we propose the following simple empirical correction to the second-order term in order to make the contact value of the RDF of the Mie potential accurate at low temperature: where γ c is a correction factor that is a function of both density and temperature, where θ = exp(β ) − 1, and α is the reduced constant defined in Eq. (18). The coefficients φ 7, i (i = 0, 3, and 4) are obtained by correlating the contact values of the RDF of the LJ fluid, as described by the expression of Johnson et al., 63 at low temperatures and densities. The functional form of the factor is chosen such that both Eqs. (61) and (62) provide an equivalent description at high density and high temperature. We find that a temperature correction, characterized by the parameter θ , is necessary to capture the RDF of LJ fluids accurately at very low temperature (cf. Fig. 4(b)). In this sense the expression for g 2 given by Eq. (62) represents an effective term which accounts for the contributions of the higher-order terms in the perturbation expansion. A hyperbolic tangent is used in order to bring the correction γ c to zero for very shortranged potentials. This allows one to capture the decreasing magnitude of the density dependence of the slope of the RDF at zero density for Mie potentials with increasing repulsive exponents. We illustrate this effect in Fig. 4(c) where our corrected second-order RDF is compared with contact values of the RDF of LJ (12-6) and Mie (30-6) fluids obtained by molecular simulation. 63 Overall the agreement is satisfactory in both cases, with a slight over-estimation of the contact values of the RDF at higher densities; nevertheless, the theory enables one to capture accurately the slope of the RDF at high density and its dependence on the repulsive part of the potential. To assess further the accuracy of the corrected second-order contact value of the RDF, our new description is compared in Fig. 5 with the empirical correlation to simulation data for the LJ fluid 63 as well as with the corresponding approximate expressions used in the development of the SAFT-VR Mie 2006 EOS 72 and the SAFT-VR LJC EOS. 65 It can be seen that our corrected second-order RDF is in much better agreement with the simulation data than the other approximations, and is now seen to approach the correct limit at high temperatures; this is a result of using of the correct contact distance σ in the zeroth-order term instead of the approximate diameter of the reference HS system d employed previously. 72 Overall, the reformulation proposed here is accurate over the entire fluid region including low temperatures.
A key aspect of this work is that the new expression presented for the contact value of the RDF is generic and applicable for any value of the repulsive and attractive exponents of the Mie potential (so long as λ a > 4), as opposed to the correlation of Johnson et al. 63 which is appropriate only for the LJ fluid.

Algebraic expression for g Mie (σ )
It is now straightforward to obtain an algebraic expression for the first-order contribution to the RDF at contact as the integral in Eq. (56) can easily be expressed as function of the first-order terms of the Sutherland fluid a S 1 (λ). The treatment outlined in Sec. III B 1 is followed whereby The second-order term requires a somewhat lengthier derivation. The details are omitted and we give only the final expression, where The expression for the contact value of the RDF of the Mie fluid is then obtained by substituting Eq. (46), (64), and (65) into Eq. (43). One should note, however, that this algebraic expression is valid only for Mie potentials with values of the repulsive and attractive exponents over the interval which was used in the parameterization of a S 1 (η; λ) [5 ≤ λ ≤ 100].

IV. NEW SAFT-VR MIE EQUATION OF STATE
With the Helmholtz free energy and the RDF at contact of the Mie fluids in hand, the traditional SAFT 41, 42 formalism for the equation of state of fluids of associating chain molecules formed from Mie segments can be implemented in a straightforward manner. We summarize here the expressions for pure substances. The extension to mixtures with the appropriate mixing rules is described in the Appendix.

A. Ideal contribution
The ideal-gas term is given by 85 where 3 is taken to represent the de Broglie volume which incorporates all of the kinetic (translational, rotational, and vibrational) contributions to the partition function.

B. Monomer contribution for Mie fluids
The contribution to the Helmholtz free energy for the reference monomeric system (corresponding to the atomized chains) is where a M = A M /(N s kT) is the corresponding residual Helmholtz free energy per monomer. As outlined in Sec. III A we have developed a perturbation series expansion in the inverse of the temperature β = 1/(kT) up to third order: where a HS is given by Eq. (11), a 1 by Eq. (34), a 2 by Eq. (36), and a 3 by Eq. (19).

C. Chain contribution for the formation of chains of Mie segments
Assuming that the Mie monomers are tangentially bonded at r = σ to form chains of m s segments, then the chain contribution to the Helmholtz free energy can be expressed in the usual Wertheim TPT1 form as 47, 101 where the RDF at contact is obtained from with the zeroth-order term g H S d (σ ) given by Eq. (46), the firstorder term g 1 (σ ) by Eq. (64), and the second-order term g 2 (σ ) is obtained from Eqs. (65) and (66).

D. Association contribution for bonding between Mie segments
The contribution to the Helmholtz free energy due to association between chains of Mie segments with s bonding site types, n a sites of type a, a = 1, . . . , s, is obtained from the usual Wertheim TPT1 expression: [101][102][103] a ASSOC = s a=1 n a ln X a − 1 2 X a + 1 2 .
Here X a is the fraction of molecules not bonded at sites of type a, which is obtained from the solution of the following mass-action relation: The function ab , which characterizes the strength and bonding volume of the association interaction, is given by where f ab (r ab ) = exp (−βφ ab (r ab )) − 1 is the Mayer-f function of the association potential described here with shortranged square-well potentials: The vector r ab denotes the centre-centre distance between the association sites a and b, and r c ab is the range of the association interaction; the square-well bonding sites are placed at a distance r d ab from the centres of the Mie segments. By introducing the square-well site-site interaction into Eq. (75) and performing the angle average of the Mayer-f function one can express ab as where F ab = exp(β H B ab ) − 1 and I ab is a dimensionless association integral defined as The determination of this integral therefore requires a knowledge of the RDF of the reference Mie fluid over a range of separations around the contact distance σ , so that the expression for the RDF at contact derived in Sec. III B is not sufficient. Even in the specific case of the Lennard-Jones potential, an accurate analytical expression for this function is challenging to obtain. 99,104,105 A compromise involves the use of empirical or semi-empirical correlations of the corresponding computersimulation data for fluids with association sites. The bonding integral has been evaluated numerically for LJ fluids with single SW bonding sites for a specific site position and range of interaction (bonding geometry) by Müller and Gubbins 106 using accurate tabulated values for the LJ RDF at different particle separations, the latter obtained from molecular-dynamics simulation. The results of the quadrature were then correlated to a 25-parameter empirical relation in terms of reduced density and temperature. This is, however, only appropriate for the LJ system with the specific bonding geometry, i.e., for the fixed value of range r c ab and position r d ab of the association sites relative to the centre of the LJ segment employed in the integration. Another sound treatment of the association in such models has been reported by Tang and Lu 107 who proposed a new method to solve the Ornstein-Zernike integral equation in order to obtain the RDF analytically around r = σ . This methodology has proven to be as accurate as the empirical correlation of Ref. 106 but suffers from two limitations: the approach cannot be applied directly in the case of the generic Mie potential; the procedure leads to a complicated expression for the RDF which does not allow for the integral in Eq. (78) to be evaluated analytically. A simpler route is to use a perturbation theory such as WCA 83 which leads to a simple expression for the LJ RDF while retaining a satisfactory description of the association integral I ab for a broad range of thermodynamic conditions. This type of analytical treatment can be improved by adding the correct low-density limit (LWCA) in an ad hoc manner. 108 Though the LWCA approximation is less accurate than the integral-equation method of Tang and Lu, 107 it can easily be applied to any spherically symmetrical segment and for any geometry of the association site.
As a complement to the full numerical solution of the integral in Eq. (78) with an accurate representation of the RDF of the reference Mie fluid (obtained from integral-equation theories or molecular simulation), we also propose a simple approximate expression for the algebraic determination of I ab . Our approach is based on the zeroth-order high-temperature expansion of the RDF, whereby g Mie (r) is related to the cavity function y H S d (r) of a hard-sphere fluid of diameter d through with the reference potential u 0 (r) given by Eq. (4). By substituting this expression into Eq. (78) It is important to note that an implicit dependence on temperature is retained in this expression for I ab through the temperature dependent hard-sphere diameter d, cf. Eq. (7). One can perform the integration to give where the dimensionless bonding volume is 102 The theoretical treatment of the association contribution described in this section complies with the usual Wertheim TPT1 scheme [43][44][45][46][47][48] requiring that: bonding at a given site is independent of bonding at another; multiple bonding at a given site is forbidden; and only dimeric, chain, and network aggregates are formed (with no ring-like or double bonded structures).

A. Validation of SAFT-VR Mie equation of state by comparison with molecular simulation
The preceding exposition has centred on the development of an improved theory for both the EOS and the contact RDF of monomer Mie fluids as well as the chain contribution and the association integral, collectively referred to as SAFT-VR Mie. In this section we assess and discuss the accuracy of the new theory in describing the thermodynamic properties of systems of chain molecules formed from Mie segments. At this stage the focus is on pure fluids by making use of the expressions developed in Sec. IV. A detailed analysis for mixtures is left for future work.

Vapour-liquid equilibria of the Lennard-Jones (12-6) fluid
The vapour-liquid coexistence and vapour pressure curves of the ubiquitous Lennard-Jones fluid determined by solving for the equality of pressure, temperature and chemical potential in the two phases with our reformulated SAFT-VR Mie EOS are shown in Fig. 6 and compared with the essentially exact data obtained by computer simulation. The description obtained with our third-order perturbation expansion of the Helmholtz free energy is also compared with the corresponding results with the expansion truncated at first and second order. As can be seen from the figure, the coexistence densities and vapour pressures obtained with the new SAFT-VR Mie EOS are in excellent agreement with molecular simulation. It is worth noting that the SAFT-VR Mie description is almost of identical quality to that obtained with the EOS of Johnson et al. 63 which is based on a correlation of simulation data. One of the key findings of our current work is the striking impact of the third-order perturbation term on the accurate representation of the near-critical region of the VLE, allowing for a significant decrease in the characteristic overshoot of the critical temperature which is obtained with the lower-order expansions. This feature may at first appear rather surprising as it is widely acknowledged that any analytical EOS (including those obtained with a traditional Barker and Henderson 84, 85 perturbation theory) will exhibit a classical quadratic convergence of the difference in the liquid and vapour coexistence densities near the critical temperature, in disagreement with the experimental near-cubic critical behaviour: the exponent characterizing the classical van der Waals critical behaviour is 1/2 while the universal critical exponent is ∼ 0.325. [109][110][111][112] Despite this inherent deficiency of classical perturbation expansions, we note here that our third-order expansion yields a quadratic convergence only very close to the critical temperature thus allowing for a very good description over the rest of the near-critical region without the need of a cross-over treatment. We also emphasize that experimentally the cross-over between the classic and universal behaviour is seen only within a tenth of a kelvin of the critical point.  113 Okumura and Yonezawa, 95 Orea et al., 115 and Potoff and Bernard-Brunel 96 (circles).

Vapour-liquid equilibria and isotherms of Mie (λ r -6) and (2λ a -λ a ) fluids
In Fig. 7(a) the simulated VLE for Mie (λ r -6) fluids obtained by Okumura and Yonezawa 95 for λ r = 8, 10, 15, and 20, by Lofti et al. 113 for λ r = 12, and by Potoff and Bernard-Brunel 96 for λ r = 36 are compared with the description with our new SAFT-VR Mie EOS and with that of the previous version of the theory from 2006. 72 In the first comparison, only the repulsive exponent of the potential is varied in order to assess the capability of the theory in capturing the nonconformal behaviour that is due to the steepness of the potential. As shown in Fig. 7(a), the SAFT-VR Mie EOS of 2006 provides an excellent description of the coexistence curve of the Mie (10-6) fluid, a fair description for the Mie (12-6) and  fluids, but poor agreement for the softer potentials such as the Mie (8-6) fluid, especially for the liquid branch. It is also important to point out that this previous version of the theory 72 cannot be applied to potentials with exponents that are higher than λ r = 12; the reader is referred to Sec. III A 3 for a detailed explanation of the source of this limitation. Our new SAFT-VR Mie EOS is seen to provide a very satisfactory agreement for the full range of repulsive exponents considered, as well as an excellent description of the VLE near the critical region. The SAFT-VR Mie EOS allows for an accurate description of very steeply repulsive potentials such as the Mie  fluid; one should note that this potential is very short ranged and the VLE envelope appears at very low temperatures. The high-fidelity agreement obtained with the SAFT-VR Mie EOS supports the validity of the form of the third-order expansion proposed in our current work. A slight deterioration in the theory is still, however, apparent for very soft potentials such as the Mie (8-6) fluid. We attribute this discrepancy to the representation of the reference Helmholtz free energy a 0 of the Mie potential which is mapped on to that of the hard-sphere reference a HS through the use of a simple density-independent effective diameter d; this approximation is known to lead to some shortcomings for very soft potentials. 114 We also consider the VLE of the (2λ a -λ a ) Mie fluids. The theoretical description is displayed in Fig. 7(b) along with the Monte Carlo simulation data of Orea et al. 115 Four different Mie potentials are considered (with λ a = 6, 7, 9, and 12). It can be seen that the new SAFT-VR Mie EOS provides a very accurate description of the coexistence densities, and can be applied simultaneously to very high values of both the repulsive and attractive exponents.
The accuracy of the theory in capturing the vapour pressure of the Mie fluids for different forms of the repulsive interactions is also analysed. We compare the available molecularsimulation data of Okumura and Yonezawa 95 (for λ r = 10), Lofti et al. 113 (for λ r = 12), and Potoff and Bernard-Brunel 96 (for λ r = 14 and 36) with the vapour pressure obtained using our new SAFT-VR Mie EOS (and with the previous version of the theory 72 ). Excellent agreement is found in each case with a slight discrepancy for the Mie (10-6) potential, cf. Fig. 8. The new theory can be applied without loss of accuracy up to extremely steep potentials. From this analysis it is clearly apparent that the present form of the theory with the third-order expansion greatly enhances the accuracy of the description of the exact molecular-simulation results compared to the previous version, 72 particularly in the vicinity of the critical point. The third-order term is useful in allowing one to correctly capture the properties of very short-ranged potentials (e.g., with large values of the repulsive exponent) for which the reduced critical temperature is located below T * = 1.
We have also determined pressure-density isotherms using our new theory for Mie (2λ a -λ a ) and (λ r -6) fluids to compare with the Monte Carlo simulation data reported by Orea et al. 115 As can be seen from Fig. 9 the new theory provides an excellent description of the simulation data in all cases, even for "unusual" forms of the potential such as the Mie (32-6) and Mie  fluids. It is therefore apparent that our

Vapour-liquid equilibria of Lennard-Jones (12-6) chain fluids
In this section the VLE of fluids of fully flexible chain molecules formed from tangent Lennard-Jones segments as described using our new SAFT-VR Mie EOS is presented and compared to the description obtained with the previous version of the theory, 72 as well as the SAFT-VR LJC, 65 and the EOSs of Johnson et al. 63 The vapour-liquid coexistence densities and vapour pressures for LJ chains formed from m s = 2, 4, and 8 segments are shown in Fig. 10. The coexistence envelopes are compared with the corresponding molecular-simulation data of Lofti et al. 113 for the monomer (m s = 1) and of Escobedo and de Pablo 119 for the chains, except for the dimer fluid where the more-recent simulation data of Vega et al. 120 are used; the simulated vapour pressures for chains of m s = 4 and 8 LJ segments are taken from the work of Mac-Dowell and Blas. 121 It is interesting to point out the marked discrepancies with the previous SAFT-VR Mie 2006 and SAFT-VR LJC approaches. This is an important comparison as these EOSs differ in the level at which the perturbation   theory is truncated and most critically the expression used for the RDF at contact to describe the chain contribution of the Helmholtz free energy. Such an analysis indicates that the Wertheim TPT1 treatment is particularly sensitive to approximations employed in estimating the pair correlation function for the reference fluid (the monomer Mie fluid in our case). This feature of Wertheim's theory has already been pointed out 80 where it was shown that much of the disagreement observed between a TPT1 representation and molecularsimulation data can be attributed to the approximations used to derive the RDF at contact. By contrast, the EOS 63 of Johnson et al. involves an empirical correlation of the EOS and contact value of the RDF obtained from exact simulation data for the monomer LJ fluid, and as a consequence represents the best agreement that one could obtain using a SAFT-like equation for these systems. Our new SAFT-VR EOS is derived with a chain term using a RDF treated up to second order and is therefore expected to provide a marked improvement in the description of LJ chains. Essentially, no overshoot of the critical temperature of the LJ dimer is observed. The progressive increase in the over-prediction of the critical temperature found for longer chains can be attributed to the use of the linear approximation to evaluate the many-body distribution function at contact in the chain contribution to the Helmholtz free energy. 122

Vapour-liquid equilibria of long Lennard-Jones (12-6) chain fluids
The description with our new SAFT-VR Mie EOS for the coexistence densities and vapour pressures of fully flexible chains of m s = 20, 50, and 100 tangent LJ segments is compared in Fig. 11 with the corresponding VLE obtained using the semi-empirical EOS of Johnson et al. 63 As can be seen, both EOSs provide equivalent satisfactory descriptions of the fluid-phase equilibria of these long chains. Discrepancies between the two approaches are apparent near the critical region, where a slightly larger over-prediction of the critical temperature is observed with the SAFT-VR Mie EOS, together with a deterioration of the vapour pressure with increasing chain length. These deviations are expected for very long chains and are a natural consequence of the linear approximation used to represent the many-body correlations (in terms of the pair distribution function) within the Wertheim TPT1 treatment of the chain contribution. In order to address this deficiency of the TPT1 approximation, Blas and Vega 123 have built on the work of Ghonasgi and Chapman 124 and Chang and Sandler, 125 demonstrating that an improved description of the coexistence densities and vapour pressure of chain fluids of moderate length can be obtained by considering a dimer (instead of a monomer) reference fluid, thereby including more structural information in the chain contribution to the Helmholtz free energy.

Vapour-liquid equilibria of Mie (λ r -6) chain fluids
As no molecular-simulation data have to our knowledge been reported for fluids of chain molecules formed from Mie segments, we have performed direct moleculardynamics simulations of the VLE for two model Mie (λ r -6) flexible chain molecules: three tangent segments (m s = 3) interacting through a Mie (22-6) potential; and six tangent segments (m s = 6) interacting through a Mie (16-6) potential. For details of the simulations the reader is directed to our series of papers 126-128 on the use of the Mie potential in the development of coarse-grained models. A comparison of the description obtained for VLE of these Mie chain fluids with our SAFT-VR EOS and the simulated values is shown in Fig. 12. The good agreement obtained in both cases provides a strong confirmation of the adequacy of the expression derived for the contact value of the RDF of Mie fluids when considering potentials of variable repulsive ranges. No deterioration of the description is apparent for chains with steep repulsive segment interactions such as the Mie (22-6) potential. We emphasize again that this system is a challenge for any high-temperature perturbation theory because as the repulsive exponent λ r is increased (i.e., as the range of the interaction is shortened), the VLE envelope is found at progressively lower reduced temperatures. It is therefore necessary to make use of a very accurate description of both EOS and the RDF of the reference Mie monomer fluid, which our SAFT-VR Mie treatment gratifyingly provides. As will be shown in Sec. V B an interaction potential with moderately large values of the repulsive exponent turns out to be very useful in describing simultaneously the VLE and the second-derivative thermodynamic properties of long n-alkanes.

Vapour-liquid equilibria of associating LJ (12-6) fluids
When a soft-core reference interaction is used to model the particles, association sites can interact even when fully embedded within the molecular core. An original feature of our current approach lies in the approximate algebraic expression for the association integral (cf. Eqs. (77), (83), and (84)) which allows one to take into account different geometries of the sites, as opposed to the fixed position (r d ab = 0.4σ and r c ab = 0.2σ ) of the semi-empirical correlated expressions. 106 Unfortunately, molecular-simulation results for different geometries of the sites are not readily available so we will assess only the association geometry used in the existing correlation. 106 In Fig. 13 we display the vapour-liquid coexistence densities for LJ particles with a single embedded square-well association site for two different values of the site-site association energy ( H B ab = 10 and H B ab = 20 ). For the dimerizing LJ system with an association energy of H B ab = 10 , the agreement in the coexistence envelope between theory and simulation is very good, while small discrepancies are apparent for the strongly associating system with H B ab = 20 . For these associating systems the description with the SAFT-VR Mie EOS turns out to be of similar accuracy to that obtained with other perturbation theories for the reference radial distribution function, such as the low density WCA recipe. 80

B. Application of the SAFT-VR Mie EOS to real substances
In our current formulation of the SAFT-VR Mie EOS, molecules are treated in a homonuclear fashion (i.e., all of the segments in the chain model representing the molecule are taken to be identical). We apply the theory to selected compounds, including members of three homologous series (nalkanes, n-perfluoroalkanes, and n-alkanols), carbon dioxide, benzene, and toluene. The methodology commonly applied to the determination of the molecular parameters of SAFT-like EOS is to estimate them from experimental saturation pressures P sat and liquid densities ρ sat . The task involves a minimization of a relative least-squares objective function consisting of the appropriate sums of individual residuals: where θ is the vector of parameters of the intermolecular model, and N P sat and N ρ sat are the numbers of vapour pressure and saturated liquid density data points, respectively; the superscripts exp and calc refer to experimental data points and calculated values. The saturation values are calculated with the EOS by determining the vapour and liquid densities that satisfy the equality in the pressure P = −(∂A/∂V ) N,T and chemical potential µ = (∂A/∂N) V ,T in both phases at a given temperature T. Equal weights are generally chosen for both properties in this objective function so that w 1 = w 2 = 1. Other properties can be obtained from the standard thermodynamic relations involving the appropriate derivatives of the Helmholtz free energy.
It has, however, been demonstrated 72 that the inclusion of single-phase condensed-liquid density ρ and speed of sound u data often leads to a better balance of the different thermophysical properties. This has proven to be very useful in determining intermolecular potential parameters, and in the particular case of the Mie interaction, physically relevant values of the repulsive exponent λ r are obtained. The objective function is expressed as where N ρ and N u are the number of condensed-liquid density and speed of sound data points, and w 3 and w 4 are the corresponding weights. Several values of the different weight factors have been tested, and a good balance in the accuracies of the description of the different estimated properties is obtained using w 1 = 1, w 2 = 1, w 3 = 1, and w 4 = 1/4. The resulting Mie intermolecular parameters for the substances examined in our current study with the SAFT-VR Mie EOS are reported in Table III. In the case of the n-alkane series, we find rather different values for the repulsive exponent λ r with our new SAFT-VR Mie EOS than were found with the previous version of the theory. 72 The new SAFT-VR Mie models are consistent with the force fields developed in the simulation studies of Potoff and Bernard-Brunel 96 who found that the non-bonded interactions between two CH 2 groups in nalkanes can be accurately modeled with a steeper repulsion than the classical Lennard-Jones (12-6) potential. The SAFT-VR Mie and Potoff and Bernard-Brunel 96 models are not, however, identical: in SAFT-VR Mie the bending and torsional potentials are not included explicitly and the nonbonded interactions are the same between all the spherical segments while a heteronuclear model is used by Potoff and Bernard-Brunel; 96 we also note that the theory is developed for a tangent-chain model, while the simulation force field is for highly overlapping segments. The difference in the repulsive exponent obtained with the SAFT-VR Mie EOS for eicosane (λ r = 22.926) and the value of Potoff and Bernard-Brunel 96 (λ r = 16) can also be attributed to the optimization procedure used by these authors who did not include vapourpressure data in the development of the force field. The TABLE III. The Mie (λ r -λ a ) molecular parameters of some selected compounds estimated from the experimental vapour-liquid equilibria and single-phase data with our new SAFT-VR Mie equation of state. The value of m s denotes the number of Mie segments making up the chains, σ characterizes the diameter of the segment, the strength of the interaction, λ r and λ a the repulsive and attractive exponents; for molecules with association sites, r c ab characterizes the range of the site-site interaction, H B ab the strength of the bonding, and r d ab = 0.4σ is the position of the site in all cases. The n-alcohols are modeled with a two-site association scheme corresponding to one donor and one acceptor site (denoted as 2B in the classification scheme of Huang and Radosz 136  percentage average absolute deviation AAD(%) of the theoretical description from the experimental data for the vapour pressure P sat , saturated-liquid densities ρ sat , enthalpy of vaporization H v , as well as properties in the condensedliquid region, such as the density ρ, speed of sound u, and isobaric heat capacity C p , are reported in Table IV; the reader is referred to Table V for details of the thermodynamic conditions considered for each substance and the sources of the experimental data. [129][130][131] It is worth noting the very small overall deviation obtained for both vapour-liquid equilibria and second-derivative thermodynamic properties. A separate comparison of the corresponding values of the critical temperature, pressure, and density is made in Table VI. A key advantage of the formulation of the theory to third order in the perturbation expansion is the significant reduction of the overshoot of the critical temperature. This is illustrated in Fig. 14 where the experimental vapour pressures, saturated densities, and enthalpies of vaporization of the n-alkanes from methane to n-decane are compared to the description with the SAFT-VR Mie EOS. The new theory is also significantly more accurate at low temperatures than the previous SAFT-VR Mie EOS developed in 2006. 72 We should note that for some low-temperature isotherms, more than three volume roots can sometimes exist at a given pressure as a result of the non-cubic nature of the perturbation expansion in density. These additional roots do not compromise the adequacy of the EOS as they are found inside the phase coexistence regions and correspond to unstable/metastable maxima/local minima in the Gibbs free energy of the system. Care should however be taken to ensure that the numerical solver used to identify volume roots converges to the global solution of the Gibbs free energy of the system (corresponding to the smallest and largest density roots). In our experience with the SAFT-VR Mie EOS, we have not found any instance of unphysical behaviour, such as the two separate regions of vapour-liquid equilibrium found with PC-SAFT. [132][133][134] As a further critical assessment of the new theory we compare the representation of the VLE and thermophysical properties with the descriptions obtained with other commonly employed SAFT EOSs such as SAFT-VR SW, soft-SAFT, and PC-SAFT. The comparison with SAFT-VR SW is relevant since both equations are developed with a BH perturbation theory for the reference monomeric fluid. An improvement in the description can in a large part be attributed to the choice of the potential used to model the reference fluid, i.e., the more-realistic soft-core Mie potential rather than the hard-core square-well potential. The soft-SAFT EOS also affords a useful comparison as the theory is formulated using a TABLE IV. Average absolute deviation AAD(%) of the description obtained with our new SAFT-VR MIE equation of state compared to experimental data for the vapour pressure P sat , the saturated-liquid density ρ sat , the enthalpy of vaporization H v , single-phase condensed-liquid density ρ, speed of sound u, and isobaric heat capacity C p . The definition of the AAD(%) is given in the caption of Table I semi-empirical EOS and RDF of the reference LJ fluid which are correlated to molecular-simulation data. An analysis of the relative performance of the SAFT-VR Mie and soft-SAFT EOSs is of interest since the deviations observed from the experimental data are therefore mainly a result of the choice of intermolecular potential model (a generalized LJ or a fixed LJ (12-6) form) and not to approximations used to evaluate the EOS and RDF of the LJ reference. The AADs(%) of the theoretical values from the experimental data obtained for n-hexane, n-pentadecane, and n-eicosane for fluid-phase equilibria and selected thermodynamic properties are reported in Tables VII and VIII. The PC-SAFT, SAFT-VR SW and SAFT-VR LJC molecular parameters for these compounds are taken from previous work. 56,65,71,72 The soft-SAFT parameters for n-eicosane are not available in the literature and are only reported as correlations in terms of the molar mass; furthermore, in order to present a fair comparison, parameters of the LJ chain segments considered in this approach are re-determined by estimation to experimental vapour pressures and saturated-liquid densities. It is important to recall that one should avoid using speed of sound data in the parameter estimation procedure with "traditional" SAFT approaches that incorporate a fixed model of the repulsive interactions. Indeed, it has been shown 72 that one is unable to represent the speed of sound with these EOSs while retaining an accurate description of the vapour-liquid equilibrium. Our goal here is not to provide a comprehensive comparison for the entire n-alkane series but just to confirm the conclusions of the previous study 72 with the new version of SAFT-VR Mie: the versatility of a potential model with a variable repulsive exponent is required to provide accurate simultaneous descriptions of the VLE and other thermodynamic properties, something which is not possible with a fixed LJ or SW form of interaction.
The second-derivative properties are closely interrelated through standard thermodynamic identities. For example, the isobaric and isochoric heat capacities, C p and C v , respectively, depend on the coefficient of thermal expansion α p = (1/V )(∂V /∂T ) p and isothermal compressibility κ T = −(1/V )(∂V /∂P ) T through the relation C p − C v = V T α 2 /κ T . We have already shown that the isobaric heat capacity, thermal expansivity, and isothermal compressibility are accurately represented with the SAFT-VR Mie EOS. In view of the thermodynamic interrelationships between these properties it follows that the isochoric heat capacity (a property that is often poorly represented with commonly employed EOSs) should also be accurately described with the theory. The temperature and pressure dependence predicted for C v with the SAFT-VR Mie approach for selected hydrocarbons (n-hexane, n-octane, n-decane, and toluene) is  Fig. 15, where the description is seen to be very satisfactory for the vapour and liquid branches, including the discontinuity at the vapour-liquid transition. We should point out that neither the heat capacity nor the thermal expansivity or isothermal compressibility are considered in the parameter estimation procedure so this represents a true indication of the predictive capability of our methodology. We now turn our attention to the n-perfluoroalkane homologous series. These fluorinated molecules are very challenging to describe with traditional intermolecular potential models due to their unique properties related partly to their high polarizability and repulsive nature of the fluorocarbonfluorocarbon interactions. 40,135 The molecular parameters for the first six n-perfluoroalkanes of the series are reported in Table III. In order to take into account the multipolar effects in an effective way for this class of substance, it is reasonable to relax the criterion on the value of the attractive exponent λ a = 6, which can now be considered as an extra adjustable parameter; one should also recall that the segments forming the chains are coarse-grained representations of the molecules comprising a number of carbon and fluorine atoms so a simple London dispersion interaction is not necessarily the most-appropriate form. This is done in the case of perfluoromethane and perfluoroethane; for the sake of convenience the value of λ a = 5.7506 obtained for perfluoroethane is transferred to higher members of the series. It is interest-ing to note that very repulsive Mie potentials (with repulsive exponents in the range λ r ∼ 19 to 42) are required to correctly capture the VLE and the speed of sound of these compounds. This is a particular feature of the n-perfluoroalkanes which are characterized by very low values of the speed of sound in the liquid phase compared to their nonfluorinated n-alkane analogues. The theoretical description of the VLE and heat capacity of perfluoromethane with the PC-SAFT and SAFT-VR Mie EOSs are compared in Fig. 16. In the case of PC-SAFT, the model parameters with the number of segments fixed to m s = 1 (perfluoromethane is a near spherical molecule) have not previously been reported, and we have determined them in the usual way by estimation to the experimental saturated liquid densities and vapour pressures. The resulting values are σ = 4.1712 Å and /k = 194.51 K. It can be seen that one is not able to adequately capture the VLE of perfluoromethane with PC-SAFT (cf. Fig. 16(a)). The AAD(%) from the experimental data for the vapour pressure is 17.10% with PC-SAFT compared to 0.92% with SAFT-VR Mie. Three different isobars (10, 25, and 50 MPa) are also presented for the density and isobaric heat capacity in Figs. 16(a) and 16(c), respectively. A significantly improved description is achieved with the SAFT-VR Mie EOS. We also illustrate the excellent agreement between the experimental data and the SAFT-VR Mie approach for other perfluorinated compounds in Fig. 17. A key feature of our new theory is that  it provides the capability of describing the VLE of pure substances even close to the critical region. Associating substances have also been considered by assessing the adequacy of our approach for the first members of the n-alkanol series. In the current work we investigate a two-site association scheme corresponding to one donor and one acceptor site (denoted as 2B in the classification scheme of Huang and Radosz 136 ). Seven molecular parameters need to be determined for each of these compounds: the number of segments m s , the diameter σ , the dispersion energy , the repulsive exponent λ r , the attractive exponent λ a , and the parameters associated with the site-site square-well potential including the association energy H B ab , the distance to the centre of the Mie segment r d ab , and the range of interaction r c ab . The attractive range is kept at the London value of λ a = 6 in order to reduce the number of adjustable parameters; the position of the site is also fixed to a constant value and set to r d ab /σ = 0.4, i.e., the sites are placed close to the surface of the Mie segment. This is the same site geometry that is used in the soft-SAFT approach and has proven to be very successful in modeling a wide variety of associating substances. The SAFT-VR Mie molecular parameters for methanol, ethanol, propan-1-ol, and n-butan-1-ol are reported in Table III. Once again, it can be seen (cf. Table IV) that one is able to accurately capture both the VLE and the second-derivative   properties over a wide range of thermodynamic conditions with the SAFT-VR Mie EOS. A relatively soft Mie effective potential is found to be required to model these alkanols underlining the usefulness of developing an equation of state for a broad range of repulsive exponents, in accordance with previous findings. 73 A preliminary assessment of the adequacy of the new SAFT-VR Mie EOS in describing the vapour-liquid equilibria of binary mixtures is now presented. Mixtures of ethane + n-decane and carbon dioxide + n-decane are examined here, and the resulting description is compared with the existing studies based on a SAFT-VR SW description. In order to make a fair comparison of the two equations of state we first estimate the optimal unlike interaction dispersion energy parameters from the experimental VLE data. In the case of SAFT-VR SW, a deviation from the geometric mean of k ij = 0.0204 is obtained for ethane + n-decane, and k ij = 0.1203 for carbon dioxide + n-decane. The corresponding unlike energy parameters obtained with SAFT-VR Mie are k ij = −0.0222 for ethane + n-decane and k ij = 0.0500 for carbon dioxide + n-decane.
The resulting theoretical description with the SAFT-VR SW and SAFT-VR Mie EOSs are displayed in Fig. 18. While one is able to correlate the experimental data at low pressures with both approaches, it is striking to see the marked improvement obtained with the SAFT-VR Mie EOS at high pressures near the critical region. To better illustrate this point, we display the corresponding vapour-liquid critical lines of both mixtures in Fig. 19 together with the corresponding experimental data. 137,138 The significant improvement seen near the critical region is not only a consequence of the choice of a more-realistic pair potential, but also to the fact that the new theory makes use of a perturbation expansion to higherorder than the SAFT-VR SW EOS providing a much better description particularly in the critical region. This is depicted in Figs. 18 and 19 where calculations obtained with the full SAFT-VR Mie perturbation expansion are compared with those obtained when only a first-order expansion of the RDF at contact is used to evaluate the chain term.

VI. CONCLUSIONS
We have developed a new, accurate SAFT-VR Mie EOS based on a reference monomer system characterized by a Mie potential of variable repulsive and attractive range. A number of important modifications have been made in the theory for both the equation of state and radial distribution function of the reference fluid compared to previous work. 72 In our current paper the Barker and Henderson 84, 85 high-temperature perturbation expansion is applied up to third order for the Helmholtz free energy of the reference Mie fluid, and up to second order for the corresponding RDF at contact. Analytical expressions for the second-and third-order contributions to the Helmholtz free energy are obtained by making use of Monte Carlo simulation data (to obtain a highly accurate fluctuation term) and VLE data up to the critical point (to determine the form of the third-order contribution). A key feature of the new SAFT-VR Mie EOS is the use of a thirdorder contribution which allows one to reproduce accurately the critical point of the Mie family of potentials. A morerigorous derivation of the contact value of the RDF for Mie fluids is also presented. The high sensitivity of the Wertheim TPT1 description of the chain contribution to approximations that are made for the structure of the reference system is emphasized. The expression proposed here for the RDF is shown to capture accurately the non-conformal behaviour for chains of tangent Mie spheres for a broad range of values of the repulsive exponent. Good agreement is obtained with molecular simulation for chains of up to 8 segments; the description is found to become slightly less accurate for very long chains of up to 100 segments. This is not surprising however, and is mainly a consequence of the use of the linear approximation in the TPT1 treatment to evaluate the m-body distribution function at contact. A new parameterization of the meanattractive energy of the hard-core Sutherland potential of variable range λ is developed which allows for the treatment of very short-ranged potentials with attractive/repulsive exponents spanning (5<λ<100). In view of the accuracy of the theory in the representation of the thermophysical properties of model systems, the adequacy of the SAFT-VR Mie EOS in describing the vapour-liquid equilibria and second-derivative properties of selected real substances (n-alkanes, fluorinated molecules, n-alcohols, and other compounds including carbon dioxide, benzene, and toluene) is investigated. In particular we demonstrate the versatility of varying the steepness of the Mie potential to represent correctly the softness/harness of the force field, e.g., the steep repulsive interactions between perfluorinated hydrocarbons. As Mie 19 and others had recognized over a century earlier, the second-derivative thermodynamic properties (in Mie's case the compressibility) are found to be very sensitive to the precise form of the repulsive interactions. Our SAFT-VR Mie approach based on potentials of a Mie (rather than a fixed LJ) form therefore offers an enhanced capability in the simultaneous description of both VLE and derivative properties. It should also be noted that, in general, the new theory allows for a much better description of the critical region of non-associating molecules compared to previous work. We highlight this significant improvement by examining the fluid-phase behaviour of two prototypical binary mixtures which exhibit large regions of fluid-fluid criticality, namely, ethane + n-decane and carbon dioxide + n-decane. In both cases, the SAFT-VR Mie EOS provides a very satisfactory description of the entire phase envelope up to the critical point, despite the fact that no special treatment of the critical behaviour is considered here. An in-depth study of the performance of the theory for more-complex mixtures, including challenging associating substances, will be presented in future publications.
The main aims that we set ourselves in Sec. I have now been achieved: to develop a highly accurate analytical description of the thermodynamic and structural properties of associating chain molecules formed from Mie segments; to use the resulting SAFT-VR Mie EOS to describe accurately the properties of real substances and thus uncover the empirical form of the intermolecular potential which arises from subtle quantum-mechanical effects for a number of representative examples. The homonuclear formulation inherent in the SAFT-VR Mie approach is being extended to heteronuclear molecules (formed from segments of different types) to represent chemically distinct groups, which is the basis of the SAFT-γ Mie group contribution approach 139 originally developed for heteronuclear molecules formed from SW segments. 140,141 Our accurate SAFT-VR Mie approach is also proving invaluable in the development of coarse-grained intermolecular potential models for a wide variety of apolar and polar systems. [126][127][128] The use of a versatile potential based on the Mie interaction of variable form offers a very significant advantage in the description of coarse-grained interactions. The repulsive and attractive exponents obtained for coarse-grained (multi-atomic) models turn out to be very different to those characterizing the smaller united-atom segments (traditionally represented with a LJ form), something that Maxwell 13 had alluded to as early as 1867 in his analysis of the form of the repulsive interactions from studies of the transport properties of gases.

ACKNOWLEDGMENTS
We are very grateful to Alejandro Gil-Villegas and Fernando del Río for useful input into aspects of the historical perspective, to Felix Llovell for his assistance with the soft-SAFT calculations, and Patrice Paricaud for helpful discussions regarding the fluctuation contribution. T.L. and C.A. thank the Engineering and Physical Sciences Research  The same treatment is applied to the second-order fluctuation term a 2 of the mixture: x s,i x s,j a 2,ij , where a 2,ij = 1 2 2λ a,ij 0,ij a S 1,ij (ρ s ; 2λ a,ij ) + B ij (ρ s ; 2λ a,ij ) − 2x λ a,ij +λ r,ij 0,ij a S 1,ij (ρ s ; λ a,ij + λ r,ij ) + B ij (ρ s ; λ a,ij + λ r,ij ) + x 2λ r,ij 0,ij a S 1,ij (ρ s ; 2λ r,ij ) + B ij (ρ s ; 2λ r,ij ) . (A20) In this expression K HS is the isothermal compressibility of the mixture of hard spheres (cf. Eq. (16)), and χ ij is given by x s,i x s,j σ 3 ij , and Finally, the third-order perturbation term is obtained from the following simple expression: x . (A25) The functions f k (k = 1, . . . , 6) in Eqs. (A22) and (A25) are obtained with the generalized form of Eq. (20) using the coefficients defined in Table II:

Chain contribution for the formation of mixtures of chains of Mie segments
The Helmholtz free energy due to chain formation can be obtained as a sum of contributions from the RDF of the Mie fluid evaluated at contact, 101 corresponding to tangentially bonded segments: An analytical expression for g Mie ij (σ ij ) can be obtained by applying the MX1b 70 mixing rule to the perturbation expansion, Eq. (72), described in Sec. III B, so that g Mie ij (σ ij ) = g H S d,ij (σ ij )exp β g 1,ij (σ ij )/g H S d,ij (σ ij ) The expression for g H S d,ij (σ ij ) is obtained from g H S d,ij (σ ij ) = exp k 0 + k 1 x 0,ij + k 2 x 2 0,ij + k 3 x 3 0,ij , (A29) where the density-dependent coefficients k i (for i = 0, . . . , 3) are calculated from the pure-component relations given in Eqs. (47)-(50) using a VDW-1 treatment: The first-order contribution to the RDF is obtained in terms of the mean-attractive energy of the corresponding Sutherland term by expressing Eq. (64) in its generalized VDW-1 form: ∂a 1,ij ∂ρ s − C ij λ a,ij x λ a,ij 0,ij a s 1,ij (ρ s ; λ a,ij ) + B ij (ρ s ; λ a,ij ) ρ s + C ij λ r,ij x λ r,ij 0,ij a s 1,ij (ρ s ; λ r,ij ) + B ij (ρ s ; λ r,ij ) ρ s .
The second-order RDF is similarly obtained (cf. Eqs. (65) and (66)) in terms of the fluctuation of Sutherland potentials as