Discrete perturbation theory for continuous soft-core potential fluids , , , , , L. A. Cervantes , G. Jaime-Muñoz , A. L. Benavides , J. Torres-Arenas , and F. Sastre

Citation: The Journal of Chemical Physics 142, 114501 (2015); doi: 10.1063/1.4909550 View online: http://dx.doi.org/10.1063/1.4909550 View Table of Contents: http://aip.scitation.org/toc/jcp/142/11 Published by the American Institute of Physics

THE JOURNAL OF CHEMICAL PHYSICS 142, 114501 (2015)

Discrete perturbation theory for continuous soft-core potential fluids L. A. Cervantes,1,a) G. Jaime-Muñoz,2,b) A. L. Benavides,2,c) J. Torres-Arenas,2,d) and F. Sastre2,e) 1

Departamento de Infraestructura, Universidad de Guanajuato, Noria Alta, Guanajuato, CP 36000, México División de Ciencias e Ingenierías, Universidad de Guanajuato, Loma del Bosque 103, Colonia Lomas del Campestre, León, Guanajuato, CP 37150, México 2

(Received 30 September 2014; accepted 9 February 2015; published online 16 March 2015) In this work, we present an equation of state for an interesting soft-core continuous potential [G. Franzese, J. Mol. Liq. 136, 267 (2007)] which has been successfully used to model the behavior of single component fluids that show some water-type anomalies. This equation has been obtained using discrete perturbation theory. It is an analytical expression given in terms of density, temperature, and the set of parameters that characterize the intermolecular interaction. Theoretical results for the vapor-liquid phase diagram and for supercritical pressures are compared with previous and new simulation data and a good agreement is found. This work also clarifies discrepancies between previous Monte Carlo and molecular dynamics simulation results for this potential. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4909550]

I. INTRODUCTION

Since the beginning of this century, the interest in understanding the phenomenon of liquid and glass polymorphism in single component fluids has attracted the attention of the scientific community. For instance, it is well known that helium exhibits two critical points: the common gas-liquid and a liquid-liquid one. Recently, experimental and/or theoretical evidence has been published, showing that other important substances could also show this anomalous behavior. Besides, the appearance of more than one type of glass has also been observed in real single component systems.1 For the case of some single components substances, there are some good atomistic potential models that are able to reproduce some of the water-type anomalies and polymorphism, but they can only be studied by simulation techniques that require high performance computers.2,3 Another way of studying these atypical systems is by using some effective isotropic pair potentials models that can be handled by different statistical-mechanics techniques. It seems that a good effective potential, able to predict some of these anomalies must have two characteristic lengths, as, for example, the discrete and soft-core continuous versions of the squareshoulder + square-well (SS + SW) and the Jagla ramp potentials.4–16 In the context of fluids, a good and efficient statisticalmechanics technique able to generate analytical equations of state is the perturbation theory.17–19 An example of this approach is Discrete Perturbation Theory (DPT) which has been used to study non-polar and polar systems.20,21 Examples of DPT applications for systems that exhibit water-type a)Electronic mail: [email protected] b)Electronic mail: [email protected] c)Electronic mail: [email protected] d)Electronic mail: [email protected] e)Electronic mail: [email protected]

0021-9606/2015/142(11)/114501/8/$30.00

anomalies are the discontinuous SS + SW potential10 and the Jagla ramp potential.11 For the discontinuous SS + SW potential, DPT has predicted multiple fluid-fluid transitions but not some other water-like anomalies, e.g., the maximum density anomaly. Besides, it has been shown that DPT can be applied not only to hard-core discontinuous potentials but it also has a good performance for hard-core continuous potentials.22,23 The main advantages of using DPT are (1) it provides an analytical expression for the Helmholtz free-energy as a function of density, temperature, and the set of parameters that characterize the intermolecular potential, which allows to obtain all the thermodynamic properties in a straightforward way, and (2) it is a very efficient tool to generate thermodynamic properties when compared, for instance, with simulation techniques. Using a molecular dynamics (MD) study, Franzese8,9,24,25 has found that a soft-core continuous version of this SS + SW potential can predict both multiple critical points and the maximum density anomaly for some particular set of potential parameters. Besides, recently, integral equations theory has been used for this potential together with Monte Carlo (MC) simulations. Some discrepancies between MD and MC data have been found for the critical points.16 To our knowledge, no analytical equation of state has been developed for this potential. The DPT application to soft-core continuous potentials has only been done for the Lennard-Jones potential.20,26 The Franzese potential is a suitable example to test the performance of the theory for other soft-core potentials. In this first study, DPT will be used to analyze the vapor-liquid phase diagram near the critical point and single phase pressures at the super-critical region. We have selected this study region because DPT is an inverse temperature perturbation expansion and as a consequence it converges faster at high temperatures. Even though we are limiting our

142, 114501-1

© 2015 AIP Publishing LLC

114501-2

Cervantes et al.

J. Chem. Phys. 142, 114501 (2015)

study to the critical and super-critical zone, super-critical fluids have interesting properties that can be used in the design of new processes: they behave near the critical point as solvents, exhibiting a high compressibility and with transport properties intermediate between gases and liquids. Besides, in an astronomical context, they could model the physical properties inside the bigger gaseous planets.27 The supercritical region has also been recently studied with more accurate experiments and simulations and its topology is more intriguing than expected.28 The work is organized as follows. In Sec. II, two different versions of the Franzese potential are presented. Section III provides a brief description of DPT. For both potentials, in Sec. IV, DPT predictions for critical data and super-critical pressures are given together with a comparison with new and available simulation data. The article is closed in Sec. V with the main conclusions of this work.

III. DISCRETE PERTURBATION THEORY

II. THE FRANZESE POTENTIAL

The Franzese potential8,9,24 is a continuous version of a discontinuous hard-core SS + SW potential and can be expressed reduced with respect to an appropriate energy scale as UR∗ U(x) = u F (x) = |U A | 1 + exp(∆(x − R∗R )) − U A∗

exp *− ,

(x − R∗A)2 2δ2A

FIG. 1. The Franzese potential u F (solid line) and the modified Franzese potential u FMD (dashed line).

( ) 24 + + U∗ 1 , A x -

(1)

l n A AHS (η)   m = + β N kT N kT m=1 i=1   × aSm (η, λi , ϵ i ) − aSm (η, λi−1, ϵ i ) ,

(3)

HS

where x = r/σ, with σ being the corresponding SS + SW hard-core diameter. The reduced energy parameters are defined as UR∗ = UR/|U A |, U A∗ = U A/|U A | = 1. For this potential, the reduced characteristic lengths are R∗R = RR/σ and R∗A = R A/σ, ∆ is a parameter that controls the steepness of the repulsive part of the potential, and δ A is the approximate attractive width for the Gaussian function defining the attractive well. In the limit of high ∆ values, the SS + SW potential that originates the Franzese potential is recovered. As can be seen, this potential represents a family of potentials characterized by a set of parameters. In previous works, the selected parameters were those of the original paper,8 and for different ∆ values. In this work, we will consider the same cases: UR∗ = 2.0, R∗R = 1.6, R∗A = 2.0, δ∗2 A = 0.1 and ∆ = 15, 30, 100, 300, 500. It is important to remark that Franzese, in order to treat this potential with MD, added a constant and a linear term to the original potential to have both the potential and its derivative equal to zero at the cutoff (x c = 3.0). This new version of the potential for the case of ∆ = 15 is uFMD = u F (x) + C + Bx,

For a system of N spherical particles of diameter σ contained in a volume V , which are interacting with an arbitrary discrete potential, DPT expresses the excess Helmholtz freeenergy as a high-temperature expansion to l th-order20 as

(2)

with C = 0.208 876 and B = −0.067 379 4. The same constants were used for different ∆ values since all satisfy with a good precision the same conditions at the cutoff.29 In this work, we will consider both potentials, u F (x) and uFMD(x). In Figure 1, both potentials are shown for the case ∆ = 15. As can be noticed, the potentials are not identical.

where A is the free energy for the hard-sphere reference fluid, β = 1/kT, k is the Boltzmann’s constant, T is the temperature, n is the total number of steps conforming the discrete potential, η is the hard-sphere packing fraction η = (π/6)ρ∗, where ρ∗ = (N/V )σ 3 is the reduced density, and aSm are the m th-order perturbation terms for a SW (ϵ i < 0) or square shoulder (ϵ i > 0) fluid. Equation (3) only requires the knowledge of the SW freeenergy perturbation terms since we assume that the SW and SS potentials are related through aiSS(η, λ, ϵ) = (−1)i aiSW(η, λ, ϵ).

(4)

The excess Helmholtz free-energy can be rewritten in a more compact expression as an expansion of the inverse of the reduced temperature, T ∗ = kT/|U A |, which is given by Ae x a1(η, λi , ϵ i ) = a H S (η) + N kT T∗ a2(η, λi , ϵ i ) a3(η, λi , ϵ i ) + + + · · ·. (5) T ∗2 T ∗3 From this expression, one can obtain any thermodynamic property. ( e x )For instance, the compressibility factor, Z = 1 + η ∂a∂η , can also be expressed as a high-temperature expansion, ae x ≡

1 z1(η, λi , ϵ i ) T∗ 1 1 + ∗2 z2(η, λi , ϵ i ) + ∗3 z3(η, λi , ϵ i ) + · · ·, T T ( ) i where z i = η ∂a ∂η . Z = 1 + z H S (η) +

(6)

114501-3

Cervantes et al.

J. Chem. Phys. 142, 114501 (2015) 3

Pσ A reduced pressure, P∗ = |U , can be obtained from A| ∗ ∗ ∗ Z, since P = Z ρ T . The reduced chemical potential can be obtained by adding the Helmholtz free-energy and the compressibility factor, µ∗ = µ/kT = a + Z. These expressions are useful to obtain the coexistence curve. We solved numerically the system of non-linear equations obtained from gasliquid equilibrium conditions at a fixed reduced temperature T ∗,

P∗(ρ∗gas) = P∗(ρ∗liquid), µ∗(ρ∗gas) = µ∗(ρ∗liquid).

(7)

Up to this point, DPT expressions are suitable for any hard-core plus an arbitrary discrete potential. Since in this work we are interested in two soft-core potentials, DPT requires a new approximation to account of the potential softness. One possibility is to use the Barker and Henderson relation18 that provides an expression for a temperature dependent diameter  σ d= (1 − e−βu(r ))dr. (8) 0

New reduced quantities can be defined with this new Pd 3 . All above DPT effective diameter, ρ △ = ρd 3 and P △ = |U A| thermodynamic properties remain the same, except that they will be evaluated using this new reduced density. The implementation of DPT to the Franzese potentials requires a discrete version of the potentials in terms of step functions. (See Figure 2 for a picture of this kind of discretization.) The number of partitions n was estimated as n = (x c − 1)/b, where b defines the width of the steps and we selected b = 0.14. In this work, the middle point on each step was used to evaluate the potential. Besides, DPT requires a cutoff selection. Since the range of both potentials is fixed to 3 so the same value was used in DPT. In Eq. (5), the hard-sphere term has been evaluated using the Carnahan-Starling equation.30 The terms aiSW are the main ingredient in DPT and can be obtained from different sources

(simulation data, integral equations, perturbation theory, etc.). DPT requires the knowledge of these expressions to any order and for any λ value up to the selected cutoff, however, this information is limited in the literature. There are only a few available analytical SW expressions according to the range of the potential and at most to fourth-order.31–33 Third- and fourth-order terms in some of these expressions are obtained from the correlation of simulation data with high uncertainties and in general are not so accurate for all the ranges and densities required by DPT. For the potentials under study, because the selected cutoff is 3, no long-range terms ai are required. Besides, since the selection of the step width was 0.14 and because the middle point of the first step is at 1.07, DPT only requires terms ai for 1.07 ≤ λ ≤ 3.0. In this work, the expressions for a1, a2, and a3 from Espíndola-Heredia et al.32 have been used. We have assumed the validity of these expressions for 1.07 ≤ λ ≤ 3.0. The explicit form of these ai terms is given in the Appendix.

IV. RESULTS

In order to test the performance of DPT applied to potentials u F and uFMD, supercritical pressures and vaporliquid phase diagrams will be compared against available and new complementary simulation data. The new simulations were MC type only for the potential uFMD for the case ∆ = 15 in the NVT ensemble in order to obtain supercritical pressures. The number of particles was N = 1372 in a cubic box with periodic boundary conditions. In our simulations, the system reached equilibrium after 25 000 cycles. Each cycle consisted of N attempts of particle movements. The averaged results were obtained over 50 000 cycles. In all runs, the trial move acceptance ratios were always around 40%. To estimate the Barker and Henderson diameter, we solved numerically Eq. (8). The effective diameter was estimated for both potentials for different ∆ values and for 0.8 ≤ T ∗ ≤ 10. For this calculation, we have used the continuous version of the potentials. Results were very similar for both potentials and all ∆ values considered so we used the same polynomial d∗ =

FIG. 2. Discrete version of the Franzese potential. Notice that discretization starts at x = r /σ = 1.

d = 1 − 0.002 853T ∗ − 0.001 046T ∗2 + 0.000 077T ∗3. σ (9)

Although thermodynamic properties were calculated in terms of ρ △ and P △, in order to compare them with simulation data, we returned to the original reduced densities ρ∗ and pressures P∗ with the transformations: ρ∗ = ρ △/d ∗3 and P∗ = P △/d ∗3. Before exploring the super-critical region, DPT predictions for the critical points were analyzed. Since for obtaining critical data from perturbation theories, it is important to include high-order perturbation terms35,36 we have used DPT to second-order (DPT2) and to third-order (DPT3). Vaporliquid equilibrium equations for each potential in regions close to their critical temperatures have been solved. Barker and Henderson diameter was included in DPT2 and DPT3,

114501-4

Cervantes et al.

J. Chem. Phys. 142, 114501 (2015)

TABLE I. Critical vapor-liquid data for the Franzese potential from DPT and MC simulations of Huš and Urbic.16 MC ∆ 15 30 100 300 500

DPT2

DPT3

T∗

ρ∗

P∗

T∗

ρ∗

P∗

T∗

ρ∗

P∗

1.13 ± 0.02 1.19 ± 0.02 1.25 ± 0.02 1.24 ± 0.02 1.24 ± 0.02

0.085 ± 0.005 0.084 ± 0.005 0.085 ± 0.005 0.086 ± 0.005 0.085 ± 0.005

0.0255 ± 0.002 0.0290 ± 0.002 0.0310 ± 0.002 0.0291 ± 0.002 0.0304 ± 0.002

1.13 1.17 1.24 1.25 1.25

0.0711 0.0704 0.0732 0.0737 0.0737

0.0249 0.0256 0.0280 0.0283 0.0283

1.072 1.113 1.185 1.191 1.191

0.0718 0.0723 0.0758 0.0761 0.0761

0.0239 0.0251 0.0278 0.0281 0.0281

however, in this region it had no effect (similar results were found with and without it). DPT2 and DPT3 critical data are given in Tables I and II together with reported simulation data9,16 for each potential. In Table I, MC and DPT results for the u F potential are shown. The MC simulation errors on critical temperatures, densities, and pressures are at most 0.02, 0.005, and 0.002, respectively.34 Table II includes MD and DPT results for the uFMD potential. Different ∆ values were considered. We have separated MC and MD simulation data intentionally, in order to try to give a possible explanation about the discrepancies found by Huš and Urbic16 between MC and MD critical data, as will be discussed at the end of this section. In Figure 3, the vapor-liquid phase diagram (T ∗ vs ρ∗) is presented for different ∆ values. The case for ∆ = 500 was not included in this figure since it is almost equal to the case ∆ = 300, as can be seen in both tables. In this figure, DPT2 and DPT3 predictions for both potentials together with MC and MD simulation data of Hus and Urbic16 and Vilaseca and Franzese,9 respectively, are shown. As it can be observed, for each ∆ case, it is clear that MC and MD data for the critical temperature are different, the same happens with DPT predictions. DPT2 and DPT3 predictions for u F potential are closer to the central values of MC data and DPT2 and DPT3 predictions for uFMD are closer to the central values of MD data. Comparing DPT2 and DPT3 results for each potential, on average for all the ∆ cases considered, DPT2 predicts better the critical temperature when compared with its corresponding simulation data. This result was not expected since normally the inclusion of higher order terms in perturbation theories improves the predictions of critical data. For both potentials, the critical temperature predicted by DPT increases as ∆ increases, but the change ratio is clearly

slower for larger ∆ values, this behavior is consistent with the simulation results. Considering the critical density, it can be seen that MC data are within the error bars of MD data. For each potential, differences in the critical density are less pronounced between DPT2 and DPT3 than for the case of critical temperature. In all cases, DPT3 predictions are slightly better than DPT2 ones when compared with simulation values. In Figure 4, the vapor-liquid phase diagram (P∗ vs ρ∗) is presented for different ∆ values. Since MD data have big uncertainties, it is not possible to obtain conclusive information from these results. DPT2 and DPT3 predictions for u F are different from those of uFMD. Again DPT2 and DPT3 results for u F are closer to MC data. The effect of including a third-order term in DPT is not relevant for the critical pressures. For u F , and ∆ ≥ 100, DPT2 and MC data predict critical data very close to the critical data of the SS + SW potential that originates this family of potentials with critical values: Tc = 1.24 ± 0.01, ρ∗ = 0.09 ± 0.02, and P∗ = 0.03 ± 0.01.37 This behavior was expected since as ∆ increases potentials should reach as a limiting case this SS + SW potential. However, DPT and MD predictions for uFMD do not recover the critical data of this limiting potential. As can be observed from the vapor-liquid phase diagrams (Figures 4 and 5), the introduction of the third-order term in DPT has not improved the critical values when compared with simulation data. This effect of third-order terms in a perturbation expansion was not expected, but as mentioned in the introduction, it is difficult to find accurate analytical expressions for 3rd- or higher-order SW perturbation terms available in the literature. The a3 SW term used in this work comes from simulation data correlation with large uncertainties difficult to fit in a single expression as a function

TABLE II. Critical vapor-liquid data for the FMD potential from DPT and MD simulations of Vilaseca and Franzese.9 MC ∆ 15 30 100 300 500

DPT2

DPT3

T∗

ρ∗

P∗

T∗

ρ∗

P∗

T∗

ρ∗

P∗

0.95 ± 0.06 1.01 ± 0.07 1.06 ± 0.04 1.06 ± 0.05 1.08 ± 0.06

0.08 ± 0.03 0.08 ± 0.03 0.08 ± 0.03 0.09 ± 0.02 0.09 ± 0.03

0.019 ± 0.008 0.022 ± 0.008 0.025 ± 0.005 0.027 ± 0.009 0.027 ± 0.008

0.989 1.029 1.083 1.087 1.087

0.0604 0.0604 0.0626 0.0628 0.0628

0.0188 0.0196 0.0213 0.0214 0.0214

0.871 0.910 0.975 0.980 0.980

0.0653 0.0662 0.0689 0.0691 0.0691

0.0182 0.0193 0.0214 0.0216 0.0216

114501-5

Cervantes et al.

J. Chem. Phys. 142, 114501 (2015)

FIG. 3. Phase diagram (T ∗ vs ρ ∗) for the potentials considered in this work for different ∆ values. The two curves to top correspond to DPT2 (solid line) and DPT3 (dashed line) predictions for u F . The couple of curves to bottom are predictions of DPT2 (solid line) and DPT3 (dashed line) for u FMD. MD simulation data from Vilaseca and Franzese9 (solid square) and MC simulation data of Huš and Urbic16 (solid diamond).

of λ and density (see Fig. 3 of Espíndola-Heredia et al.32). Besides, DPT approximation could also be responsible of this not expected effect. In order to analyze the performance of DPT in the supercritical region, reduced pressures as a function of density for four super-critical temperatures are presented in Figure 5 for both potentials. DPT2 and DPT3 predictions overlap in all cases, so we just show DPT2 results. We also have included MC data of Huš and Urbic15 and our new simulation data. As can be seen, the agreement between DPT2 and simulation is good as long as they are compared with their

corresponding potential simulation data. In general, the theory predicts the same tendency of the simulation data, pressure rises as temperature, and density raises. As expected from an inverse temperature perturbation expansion, the agreement between DPT2 and simulation data improves as the temperature increases. In this figure, for each potential, the corresponding predictions without Barker and Henderson diameter are shown as dotted lines. As can be seen, an improvement is obtained with the inclusion of the Barker and Henderson diameter for the higher densities and temperatures, but systematic overestimation of the theory with respect to

FIG. 4. Phase diagram (P ∗ vs ρ ∗) for the potentials considered in this work for different ∆ values. The two curves to top correspond to DPT2 (solid line) and DPT3 (dashed line) predictions for u F . The couple of curves to bottom are predictions of DPT2 (solid line) and DPT3 (dashed line) for u FMD. MD simulation data from Vilaseca and Franzese9 (solid square) and MC simulation data of Huš and Urbic16 (solid diamond).

114501-6

Cervantes et al.

J. Chem. Phys. 142, 114501 (2015)

FIG. 5. Reduced pressures as a function of reduced density for the potentials considered in this work for several supercritical reduced temperatures T ∗ for the case ∆ = 15. For each case, from top to bottom, DPT2 results for u F (solid line) and DPT2 results for u FMD (dashed line). Dotted lines represent DPT2 predictions without Barker and Henderson diameter for each potential. MC simulation data of Huš and Urbic15 are shown with solid circles and MC data of this work with solid triangles up.

the simulation data is observed. A better approximation than Barker and Henderson diameter to account for the softness of the potentials is required. Discrepancies between some simulation techniques for a given potential are sometimes due to the way the original potential is modified to avoid technical simulation problems.38 For instance, Trokhymchuk and Alejandre39 discussed this point for the case of a Lennard-Jones fluid and found that a possible reason, among others, for the discrepancies between interfacial tension MC and MD predictions is that MC uses a truncated potential, while MD uses a truncated force, i.e., a derivative of the potential. The difference between potentials u F and uFMD is the term C + Bx, and this is possible the origin for the discrepancies between the two versions of the Franzese potential, which is reinforced by DPT predictions.

V. CONCLUSIONS

Discrete perturbation theory has been applied to a family of soft-core continuous potentials. The vapor-liquid phase diagram near the critical point and super-critical pressures for two versions of the Franzese potential have been obtained and compared with simulation results and a good agreement was found. Discrepancies between simulation data previously obtained by Monte Carlo and molecular dynamics for these potentials have been clarified. We are working on improving this theory to study the low temperature region where a second critical point and some interesting anomalies appear. ACKNOWLEDGMENTS

We thank Dr. T. Urbic for providing us his simulation data. We also thank the financial support from CONACYT (México): Project No. 152684 and Universidad de Guanajuato (México) Grant No. 56-060. A.L.B. also thanks CONACYT

(México) Convocatoria 2014 de Estancias Sabáticas Nacionales, Estancias Sabáticas al Extranjero y Estancias Cortas para la Consolidación de Grupos de Investigación.

APPENDIX: SQUARE-WELL HELMHOLTZ FREE-ENERGY PERTURBATION TERMS

The first three-order perturbation terms as obtained by Espíndola-Heredia et al.32 used in this work are presented. The first coefficient from Eq. (3) can be expressed as a1 =

3 

( α i,1(λ)

i=2

6η π

) (i−1) +

4 

( γi (λ)

i=1

6η π

) (i+2) ,

(A1)

here, α2,1 is given by 2π 3 (λ − 1). 3i!

(A2)

( π )2  P (λ), λ ≤ 2   1 ,  6  P1(2), λ > 2

(A3)

α2,i = − The α3,1 is given by α3,1 = −

where P1(λ) is the 6th order polynomial P1(λ) = λ6 − 18λ4 + 32λ3 − 15.

(A4)

As for the γn , their expressions are γn = γn,1λ + γn,2λ2 + Rn (λ)/Q n (λ),

(A5)

which are valid for n ≤ 4, Rn (λ) is Rn = γn,3 +

8  j=4

γn, j (λ3 − 1) j−2

(A6)

114501-7

Cervantes et al.

J. Chem. Phys. 142, 114501 (2015)

TABLE III. γ n, j coefficients.32 j

γ 1, j

γ 2, j

γ 3, j

γ 4, j

1 2 3 4 5 6 7 8 9 10 11 12 13 14

−59.046 4 26.098 26.445 4 7.401 36 11.074 3 −5.491 52 0.781823 −0.031 975 1 0.827 621 0.605 635 −0.254 959 0.037 711 1 −0.002 108 96 0.000 045 232 8

214.316 −88.139 4 273.3 95.975 9 71.122 8 −40.265 6 5.94069 −0.238 42 −2.175 58 −1.292 55 0.554 993 −0.085 754 3 0.004 925 11 −0.000 107 067

−225.479 88.820 2 250.472 90.260 6 57.027 4 −33.237 6 4.99527 −0.195 714 1.846 77 0.998 13 −0.440 314 0.070 879 3 −0.004 162 74 0.000 091 729 1

65.050 4 −25.096 74.309 5 26.215 3 18.439 7 −10.089 1 1.50243 −0.057 694 −1.871 54 −1.016 82 0.445 247 −0.072 510 7 0.004 278 62 −0.000 094 972 3

and Q n (λ) is Q n = γn,9 +

14 

( π )2  P2(λ) − P1(λ)/2, λ ≤ 2    17 . (A12)  − + P (λ), 6  λ>2 4  2 P1 is the polynomial given by Eq. (A4), P2 and P4 are also 6th order polynomials defined as α3,2 =

γn, j (λ3 − 1) j−7.

(A7)

j=10

The coefficients γn, j are listed in Table III. The second coefficient from Eq. (3) is given by ( )3  6η 6η  a2 = χ(η, λ) exp ξ2(λ) + ϕ1(λ) π π  ( ) 4 6η  + ϕ2(λ) , π  

P2(λ) = −2λ6 + 36λ4 − 32λ3 − 18λ2 + 16,

(A13)

P4(λ) = 32λ − 18λ − 48.

(A14)

3

2

Finally, the third coefficient is given by   6η 6η a3 = α2,3 exp ξ3 + K3(η, λ) , π π

(A8)

(A15)

where K3 is given by

with ( )2 6η * 6η χ = α2,2(λ) /1.5129+ , 1− (A9) π , π where α2,2 is given by Eq. (A2). The functions ϕi are given by ϕi =

7 

ϕi, n λn ,

K3(η, λ) =

whose coefficients ϕi, n are listed in Table IV. The function ξ2 is given by (A11)

where α3,2 is

1

π + 6η π

4 n n=1 θ 3, n λ

n=5 θ 3, n λ

7

n−4

.

(A16)

The coefficients θ 3, n are contained in Table IV. The function ξ3 is given by ξ3 = α3,3(λ)/α2,3(λ),

(A10)

n=0

ξ2 = α3,2(λ)/α2,2(λ),

( 6η ) 2 

(A17)

where α2,3 is given by Eq. (A2) and α3,3 is  P2(λ) − P1(λ)/6 − P3(λ), λ ≤ 2 ( π )2     α3,3 = . (A18) 17 6    − + P4(λ) − P5(λ), λ>2  6 P1, P2, and P4 are given by Eqs. (A4), (A13), and (A14), respectively. While P3 and P5 are defined by

TABLE IV. ϕ n, j and θ 3, j coefficients.32

P3(λ) = 6λ6 − 18λ4 + 18λ2 − 6,

(A19)

j

ϕ 1, j

ϕ 2, j

θ 3, j

P5(λ) = 5λ − 32λ + 18λ + 26.

(A20)

0 1 2 3 4 5 6 7

−1 320.19 5 124.1 −8 145.37 6 895.8 −3 381.42 968.739 −151.255 9.985 92

1 049.76 −4 023.29 6 305.95 −5 265.42 2 553.84 −727.3 113.631 −7.562 66

... −945.597 1 326.61 −471.688 ... 23.227 1 −2.634 77 ...

6

1H.

3

2

E. Stanley and P. Debenedetti, Advances in Chemical Physics, Liquid Polymorphism (Wiley, 2012), Vol. Vol. 152. 2C. Vega and J. L. F. Abascal, Phys. Chem. Chem. Phys. 13, 19663 (2011). 3C. Vega, E. Sanz, J. L. F. Abascal, and E. G. Noya, J. Phys.: Condens. Matter 20, 153101 (2008). 4P. C. Hemmer and G. Stell, Phys. Rev. Lett. 24, 1284 (1970). 5E. A. Jagla, J. Chem. Phys. 111, 8980 (1999). 6G. Malescio, G. Franzese, A. Skibinsky, S. V. Buldyrev, and H. E. Stanley, Phys. Rev. E 71, 061504 (2005). 7J. B. Caballero and A. M. Puertas, Phys. Rev. E 74, 051506 (2006).

114501-8 8G.

Cervantes et al.

Franzese, J. Mol. Liq. 136, 267 (2007). Vilaseca and G. Franzese, J. Chem. Phys. 133, 084507 (2010). 10L. A. Cervantes, A. L. Benavides, and F. del Río, J. Chem. Phys. 126, 084507 (2007). 11A. L. Benavides, L. A. Cervantes, and J. Torres, J. Phys. Chem. C 111, 16006 (2007). 12S. A. Egorov, J. Chem. Phys. 129, 024514 (2008). 13J. Y. Abraham, S. V. Buldyrev, and N. Giovambattista, J. Phys. Chem. B 115, 14229 (2011). 14Yu. D. Fomin, E. N. Tsiok, and V. N. Ryzhov, J. Chem. Phys. 135, 124512 (2011). 15M. Huš, M. Zalar, and T. Urbic, J. Chem. Phys. 138, 224508 (2013). 16M. Huš and T. Urbic, J. Chem. Phys. 139, 114504 (2013). 17R. W. Zwanzig, J. Chem. Phys. 22, 1420 (1954). 18J. A. Barker and D. Henderson, J. Chem. Phys. 47, 4714 (1967). 19J. R. Solana, Perturbation Theories for the Thermodynamic Properties of Fluids and Solids (CRC Press, 2013). 20A. L. Benavides and A. Gil-Villegas, Mol. Phys. 97, 1225 (1999). 21A. L. Benavides and F. Gámez, J. Chem. Phys. 135(13), 134511 (2011). 22J. Torres-Arenas, L. A. Cervantes, A. L. Benavides, G. A. Chapela, and F. del Río, J. Chem. Phys. 132, 034501 (2010). 23L. Escamilla, J. Torres-Arenas, and A. L. Benavides, J. Mol. Liq. 185, 20 (2013). 9P.

J. Chem. Phys. 142, 114501 (2015) 24A.

B. De Oliveira, G. Franzese, P. A. Netz, and M. C. Barbosa, J. Chem. Phys. 128, 064901 (2008). 25P. Vilaseca and G. Franzese, J. Non-Cryst. Solids 357, 419 (2011). 26G. A. Chapela, F. Del Río, A. L. Benavides, and J. Alejandre, J. Chem. Phys. 133(23), 234107 (2010). 27P. F. McMillan and H. E. Stanley, Nat. Phys. 6, 479 (2010). 28D. M. Heyes and L. V. Woodcock, Fluid Phase Equilib. 356, 301 (2013). 29G. Franzese, private communication (2014). 30N. F. Carnahan and K. E. Starling, J. Chem. Phys. 51, 635 (1969). 31S. Rivera-Torres, F. del Río, R. Espíndola-Heredia, J. Kolafa, and A. Malijevský, J. Mol. Liq. 130, 44 (2013). 32R. Espíndola-Heredia, F. del Río, and A. Malijevský, J. Chem. Phys. 130, 024509 (2009). 33A. L. Benavides and F. del Río, Mol. Phys. 68, 983 (1998). 34T. Urbic, private communication (2014). 35T. Lafittte, A. Apostolakou, C. Avendaño, A. Galindo, C. S. Adjiman, E. A. Muller, and G. Jackson, J. Chem. Phys. 139, 154504 (2013). 36A. F. Ghobadi and J. R. Elliot, J. Chem. Phys. 139, 234104 (2013). 37A. Skibinsky, S. V. Buldyrev, G. Franzese, G. Malescio, and H. E. Stanley, Phys. Rev. E 69, 061504 (2004). 38D. Frenkel and B. Smit, Understanding Molecular Simulation From Algorithms to Applications (Academic Press, 2002). 39A. Trokhymchuk and J. Alejandre, J. Chem. Phys. 111, 8510 (1999).