Spherically averaged versus angle-dependent interactions in quadrupolar fluids
Bortolo Matteo Mognetti, Martin Oettel, Leonid Yelash, Peter Virnau, Wolfgang Paul, Kurt Binder
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] M a r Spherically averaged versus angle-dependent interactions inquadrupolar fluids
B. M. Mognetti, M. Oettel, L. Yelash, P. Virnau, W. Paul and K. BinderInstitut f¨ur Physik, Johannes Gutenberg–Universit¨at Mainz,Staudinger Weg 7, D-55099 Mainz
Employing simplified models in computer simulation is on the one hand often en-forced by computer time limitations but on the other hand it offers insights intothe molecular properties determining a given physical phenomenon. We employ thisstrategy to the determination of the phase behaviour of quadrupolar fluids, wherewe study the influence of omitting angular degrees of freedom of molecules via aneffective spherically symmetric potential obtained from a perturbative expansion.Comparing the liquid-vapor coexistence curve, vapor pressure at coexistence, inter-facial tension between the coexisting phases, etc., as obtained from both the modelswith the full quadrupolar interactions and the (approximate) isotropic interactions,we find discrepancies in the critical region to be typically (such as in the case ofcarbon dioxide) of the order of 4%. However, when the Lennard-Jones parametersare rescaled such that critical temperatures and critical densities of both modelscoincide with the experimental results, almost perfect agreement between the above-mentioned properties of both models is obtained. This result justifies the use ofisotropic quadrupolar potentials. We present also a detailed comparison of our sim-ulations with a combined integral equation/density functional approach and showthat the latter provides an accurate description except for the vicinity of the criticalpoint.
I. INTRODUCTION
The study of fluid phase equilibria by computer simulation methods [1, 2, 3, 4, 5] hasbecome an extremely active field, since accurate information on thermodynamic propertiesof simple and complex fluids and their mixtures is of enormous importance for a varietyof applications [6, 7, 8]. The problem of understanding the phase behavior of such fluidsis a fundamental problem of statistical mechanics, too [9, 10]. While such properties inprinciple can be found from experiments, particularly for mixtures such data still are ratherincomplete, since a cumbersome study of a large space of control parameters (temperature T , pressure p , and mole fraction(s) x ( x α ) in the case of binary (multicomponent) mixtures)needs to be made.In simulations an economical use of computational resources often dictates the use ofmodels that are as simple as possible. The standard approach for the simulation of fluidsis to apply classical Monte Carlo and Molecular Dynamics methods [1, 2, 11] that requireeffective potentials (usually of pairwise type). Thus, both quantum effects associated withthe finite mass of the nuclei are ignored as well as the degrees of freedom of the electrons(sometimes the latter are considered, when the effective potentials are derived by “ab initio”quantum chemistry methods, see e.g. [12, 13], but sometimes these potentials are postulatedon purely empirical grounds [14, 15]).Now, even when the above approximations are accepted, there the question remainswhether an all-atom model is needed for the description of intermolecular forces, or whetherfurther degrees of freedom may be eliminated. For example, consider carbon dioxide (CO ),which is an extremely important fluid due to its use as supercritical solvent [6, 8]. Modelsused for the simulation for CO are truly abundant [14, 16, 17, 18, 19, 20, 21, 22, 23, 24,25, 26]; they include all-atom models with either flexible or rigid intermolecular distances,as well as models where a CO molecule is reduced to a point particle, with [27] or evenwithout [28] a quadrupolar moment. While the latter case, where the molecules are describedas simple point particles interacting with Lennard-Jones forces, is computationally mostefficient, it is also the least accurate. Recently it was suggested [27, 29] that a significantgain in accuracy with almost no loss in computational efficiency can be obtained by usingperturbation theory to construct an effective isotropic quadrupolar potential [30]. While verypromising results for a variety of molecular fluids, including carbon dioxide and benzene, [27]have been obtained, it still needs to be established to what extent the isotropic quadrupolarpotential actually reproduces the physical effects of the actual angle-dependent quadrupolarinteractions.In the present paper we fill this gap, using the case of CO as an archetypical example.Applying the same grand-canonical Monte Carlo techniques in conjunction with successiveumbrella sampling [31] and finite-size scaling analysis [32, 33] that were used for the workapplying the isotropic quadrupolar potential [27], we obtain for the present model (which isdescribed in Sec. 2) the phase diagram in the temperature-density and temperature-pressureplanes, as well as the interfacial tension (Sec. 3). In Sec. 4 we discuss the optimized choice ofthe Lennard-Jones parameters, and show that the differences between the models with thefull and averaged quadrupolar interactions are rather minor, when the critical temperaturesand densities are matched. Sec. 5 presents a comparison with integral equation/densityfunctional calculations and shows that the latter approach can describe isotherms as well asthe coexistence curve of the model very accurately, but is not applicable very near to thecritical point. Sec. VI summarizes our conclusions. II. MODEL AND SIMULATION METHODSA. Models with Full and Averaged Quadrupolar Interaction
We start from a system of uncharged point particles which have a quadrupole moment Q and interact also with Lennard-Jones (LJ) forces, U LJij = 4 ǫ "(cid:16) σr ij (cid:17) − (cid:16) σr ij (cid:17) , (1)with r ij = | ~r i − ~r j | being the distance between particles i and j at sites ~r i , ~r j , and therange and strength of the LJ potential are denoted as σ and ǫ , respectively.The quadrupole-quadrupole interaction is U QQij = 3 Q r ij f QQij , (2)with [34] f QQij = 1 − Θ i − Θ j + 17 cos Θ i cos Θ j + 2 sin Θ i sin Θ j cos (Φ i − Φ j ) −
16 sin Θ i cos Θ i sin Θ j cos Θ j cos(Φ i − Φ j ) , (3)where (Θ i , Φ i ) are the polar angles characterizing the orientation of the uniaxial moleculerelative to the axis connecting the sites ~r i , ~r j of the two particles.In order to speed up the Monte Carlo simulation, we wish to introduce a cutoff r c suchthat the total potential is zero for r > r c . This needs to be done such that the total potentialis continuous for r = r c , and the same condition should apply for the equivalent isotropicquadrupolar potential, resulting from treating the quadrupole-quadrupole interaction in sec-ond order thermodynamic perturbation theory in the partition function [29, 30]. As a result,we use the following potential (“F” stands for “full potential”) U Fij = ǫ "(cid:16) σr ij (cid:17) − (cid:16) σr ij (cid:17) − q F r(cid:16) σr ij (cid:17) − (cid:16) σr c (cid:17) f QQij + S , r ≤ r c , r ≥ r c , (4)where S = ( σ/r c ) − ( σ/r c ) . For r ij ≪ r c (and large enough r c ) Eq. (4) reduces toEqs. (1)-(3), noting the abbreviation q F = Q ( ǫσ ) . (5)Following previous investigations [5, 27, 28] we use r c = 2 · √ σ . Indeed in this work we willuse some results of Ref. [27], like the critical lines, which depend on the choice of the cutoff.Differences in the thermodynamic properties arising from a different choice of the cut-off areminor at least for simple Lennard-Jones potentials (and a suitable renormalization of theLJ parameters) [27].Note that the potential in Eqs. (1)-(3) is cut off in Eq. (4) in such a way that the cutoff r c does not depend on the angles φ i , Θ i . The corresponding isotropic spherically averaged,potential is (“A” stands for “averaged potential”) [29] U Aij = ǫ A "(cid:16) σ A r ij (cid:17) − (cid:16) σ A r ij (cid:17) − q A (cid:16) σ A r ij (cid:17) + S A , r ≤ r c , r ≥ r c , (6)where again the constant S A is chosen such that the potential is continuous for r = r c [27]. Note that for the potential, Eq. (6), the forces at r c are discontinuous but do notdiverge there. This is a requirement if one wishes to estimate the pressure from the virialtheorem, when one does a simulation in the µ VT or NVT ensemble, respectively [1, 2].In Ref. [27] we investigate extensively the use of potential (6) for modeling quadrupolarsubstances like carbon dioxide and benzene. These results were compared with prior inves-tigations of a simple LJ model [28] without quadrupolar moment in which we only matchcritical temperature and density with corresponding experimental values to obtain ǫ A and σ A . (In principle, ǫ A and σ A could also be matched at any point in the phase diagram,which by definition would improve agreement around this point - however, at the cost ofimposing an inaccurate description of the critical region.) The introduction of a sphericallyaveraged quadrupolar moment improves agreement with the experimental phase behavioursignificantly, especially for carbon dioxide. (Compare with Figs. 6, 7 and 8).The optimal choice of the parameters ǫ A , σ A and q A is somewhat subtle. A straightforwardchoice simply requires that the two potentials are strictly equivalent for temperatures T →∞ , where the perturbation expansion becomes exact. This would imply ǫ A = ǫ, σ A = σ, q A = Q k B T ǫσ = ǫk B T q F . (7)Note that for the potential U Aij the parameter q A is not a constant, but inversely propor-tional to temperature T .However, the physically most interesting region of the system is clearly not the regime T → ∞ , but rather the vicinity of the critical temperature T c . Thus, in our previous workon carbon dioxide (CO ) [27] where Eq. (6) was used, we have chosen the parameters ǫ A , σ A ,such that both T c and the critical density ρ c of the model precisely coincide with theirexperimental counterparts [35] T c, exp and ρ c, exp . Using also the experimental value of thequadrupole moment for CO , Q = 4 . ǫ A = 3 . × − J , σ A = 3 .
785 ˚A , q
A,c ≡ q A ( T c ) = 0 . , q F = 0 . . (8)In Sec. 3, we shall present numerical results for thermodynamic properties of the fullmodel, Eq. (4) and (5) with this choice of parameters, Eq. (8), and compare them to thecorresponding results based upon Eq. (6) (some of the latter results have been compared in[27] to both experiment and simulations of CO using other models).As we shall see in Sec. 3, in the critical region both models, Eqs. (4) (5) and Eq. (6) areno longer strictly equivalent to each other, as expected since the accuracy of perturbationtheory deteriorates the lower the temperature. Being interested in the critical region, it ismore natural to choose the parameters ǫ, σ and ǫ A , σ A of both models such that T c , ρ c ofboth models match their experimental counterparts. This requires necessarily a choice of ǫ and σ different from the choice implied by Eqs. (7), (8), since the latter choice would yield(choosing Lennard-Jones units ε , σ of Eq. 8) T ∗ ( F ) c = 1 . , ρ ∗ ( F ) c = 0 . , (9) T ∗ ( A ) c = 1 . , ρ ∗ ( A ) c = 0 . , (10)as shall be discussed in more detail in Sec. 3.Since the relation between q F and Q (Eq. 5) or q A and Q (Eq. 7) depends on ǫ and σ as well, finding the choice of the latter parameters which yields T c = T c, exp and ρ c = ρ c, exp is a self consistency problem [27]. In principle, one needs to record the functions T ∗ c ( q F ) /T ∗ c ( q F = 0) and ρ ∗ c ( q F ) /ρ ∗ c ( q F = 0), similarly as described in [27]. However, sincethe differences between the results of Eqs. (9) and (10) are rather small, it is a very goodapproximation to simply keep the value of q F as found in Eq. (8) and just recompute theappropriate values of ǫ and σ , which we shall denote as ǫ F and σ F , in order to distinguishthem from the choice of Eqs. (7,8). This procedure immediately yields ǫ F = 3 . × − J , σ F = 3 .
760 ˚A . (11)A good test of possible errors introduced by this approximation is provided by usingEq. (5) together with Eq. (11) to check the precise value of the physical quadrupole momentstrength Q this corresponds to. This yields Q new = 4 .
292 D˚A instead of the value used in[27], Q = 4 .
300 D˚A (note that the actual quadrupole moment of CO is negative, but sinceonly the square of Q actually matters, cf. Eq. (2), we suppress the sign of Q throughout). Ifone uses this slightly modified value Q new instead of Q in Eq. (7) together with the master-curves T ∗ c ( q A ( T c )) /T ∗ c (0) and ρ ∗ c ( q A ( T c )) /ρ ∗ c (0) calculated in [27], instead of Eq. (8) slightlyrevised estimates of ǫ A and σ A would result ǫ A = 3 . × − J , σ A = 3 .
784 ˚A , q A ( T c ) = 0 . . (12)But in view of the large error with which the actual quadrupole moment strength of CO is known [35], Q = 4 . ± . T c , ρ c , Q ) thatcoincide with their experimental counterparts, have precisely the same values. B. Simulation Methods and Tools for the Analysis of the Simulation Data
In this section, we summarize our procedures for carrying out the simulations and theiranalysis only very briefly, since more detailed descriptions for similar models can be foundin the literature [5, 27, 28].The estimation of vapor-liquid coexistence curves and critical parameters is done in thegrand canonical ( µ VT) ensemble, varying the chemical potential µ and recording the densitydistribution P L ( ρ ), and analyzing carefully the dependence on the linear dimension L of thesimulation box (as usual, we take V = L at the critical point, while deep into the coexistenceregion we use an elongated box V = 2 · L . For both geometries periodic boundary conditionsare applied). For T < T c , the value of µ coex ( T ) where phase coexistence between vaporand liquid occurs is found from the “equal weight rule” [31, 32, 33, 36]. For an accuratesampling of P L ( ρ ) including the densities inside the two phase coexistence region that alsoneed to be studied for an accurate estimation of the weights of the vapor and liquid phases,successive umbrella sampling methods [28, 31] are used, as well as re-weighting procedures[31, 32, 33, 37]. Note that the presence of the orientational degrees of freedom in Eqs. (3,4)does not constitute any principal difficulty here. The acceptance rate for the insertionof particles with a randomly chosen orientation is of the same order as for the isotropicpotential, Eq. (6), where this degree of freedom has been eliminated. This fact is understoodeasily, since the strength of the quadrupolar interaction, Eq. (2), is distinctly smaller than thestrength of the Lennard-Jones interaction, Eq. (1), for the present choice of q F . However, thetime required to compute the energy change caused by such a particle insertion or deletionis about an order of magnitude larger when Eq. (4) rather than Eq. (6) is used, due to thecomplicated angular dependence of the quadrupole-quadrupole interaction (Eq. 3).Nevertheless it is still feasible for this model, Eq. (4), to obtain sufficiently accurateinformation on P L ( ρ ) for a variety of temperatures T and lattice linear dimensions L ,following a path along the coexistence curve µ = µ coex ( T ) in the ( µ, T ) plane, and itscontinuation for T > T c (there the path is defined by the condition that the derivative( ∂ρ/∂µ ) T = L ( h ρ i T,µ − h ρ i T,µ ) is maximal). Fig. 1 shows, as an example, second andfourth order cumulants U and U along such a path as a function of temperature. Thesecumulants are defined by U = h M i / h| M |i , U = h M i / h M i , M = ρ − h ρ i , (13)where now we have omitted the subscripts ( T, µ ) from the averages h . . . i . As is well known[4, 5, 11, 32, 33], accurate estimates for T c can be obtained from the common intersectionpoint of either U ( L, T ) or U ( L, T ) for different L . The justification of this simple recipefollows from the theory of finite size scaling [38, 39, 40]. Note also that the ordinate values U ∗ , U ∗ of these common crossing points should be universal for all systems belonging to theIsing model universality class, to which both models Eq. (4) and (6) should belong, and areknown with very good accuracy [4, 41].From Fig. 1 it is evident that this intersection property does not work out perfectlywell, in particular, the curve for L = 6 . σ is somewhat off. However, finite size scalingshould become exact only in the limit where both L → ∞ and T → T c , while otherwisecorrections come into play. Systematic improvements (taking the so-called “field mixing”[4] and “pressure mixing” effects [42] into account) are possible, but are not consideredto be necessary here, since the relative accuracy of our estimate for T c extracted fromFig. 1 is clearly not worse than 3 · − , and this suffices amply for our purposes. A recentcomparative study of different finite size scaling based approaches for the study of criticalpoint estimation of Lennard Jones models [43] is in full agreement with this conclusion. Notethat the accuracy of the data in Fig. 1 is comparable to data taken for a pure LJ model[28] or for Eq. (6) [27], respectively (we estimate the relative accuracy of the curves in Fig.1 to be of the order of 0.5% or better). We have used 7 · , 3 · , 6 · and 9 · MonteCarlo steps (respectively for the
L/σ = 6 .
74, 9, 11.3 and 13.5 system) for each simulationpoint ( T ∗ ,∆ µ coex ( T ∗ )), for which the data for P L ( ρ ) were sampled, and applied histogramextrapolation methods (see [27, 28, 33] for details) to obtain the smooth curves drawn inFig. 1. In every step 100 attempts to insert or delete particle plus local moves are done.For T < T c the densities ρ (1)coex , ρ (2)coex of the two coexisting vapor and liquid phases can besimply read off from the peak positions of P L ( ρ ), and from the density minimum in betweenthe peaks the interfacial tension γ ( T ) can be estimated, using the relation γ ( T ) /k B T = 0 . L − ln[ P L ( ρ (1 , ) /P L ( ρ d )] , L → ∞ (14)where ρ d = ( ρ (1)coex + ρ (2)coex ) / δr around its old position. Simultaneously the orientation of its molecular axis is chosen insidea cone of angle δ Θ around its old orientation [1, 2]. The choices of δr and δ Θ were adjustedto have an acceptance rate of 40% for such moves.
III. DIRECT COMPARISON BETWEEN RESULTS FOR THE FULL ANDCORRESPONDING AVERAGED QUADRUPOLAR POTENTIAL
As noted in Sec. 2.1, Eq. (6) results from Eqs. (3,4) when one carries out a second orderthermodynamic perturbation theory and interprets the result as being due to an averagepotential [29]. Since thermodynamic perturbation theory is basically a high temperatureexpansion in powers of 1 /T , it is a matter of concern how accurate such a procedure reallyis in the temperature region around criticality and below.As a first test, we have carried out a NVT simulation using the averaged model, Eq. (6),at a density that is much larger than the critical density, namely ρ ∗ = 0 . p ∗ ( T ) from the virial formula (Fig. 2a). Thispressure then was used as an input for a NpT simulation of the full model, Eqs. (3,4).Fig. 2(b) shows the corresponding comparison: one sees that for large T ∗ (i.e., T ∗ ≥ . ρ ∗ that waschosen, while for T ∗ ≤ . ρ ∗ = 0 .
544 is reproduced overthe entire temperature interval shown in Fig. 2(b) with negligibly small errors, since withthe chosen volume (
L/σ A = 10 .
3) systematic discrepancies between the different ensemblesof statistical mechanics are completely negligible (although such discrepancies will occur inthe two-phase coexistence region or near the critical point). The deviations seen in Fig. 2(b)simply represent the higher order terms of the 1 /T expansion, by which the averaged and0full potentials differ. Similar discrepancies between the full and averaged potential were alsoseen on the vapor side of the coexistence curve (not shown here).As a result, coexistence densities in the ( T ∗ , ρ ∗ ) plane are rather different as well for thetwo models, Fig. 3, as expected from the differences noted in Eqs. (9,10), and a similardiscrepancy occurs between the predictions for the pressure along the coexistence curve(Fig. 4) and the interfacial tension (Fig. 5). IV. COMPARISON BETWEEN RESULTS FOR THE FULL AND AVERAGEDQUADRUPOLAR POTENTIALS WITH OPTIMIZED PARAMETERS
Now we present a comparison between the full quadrupolar plus Lennard-Jones potential(Eqs. 4, 5) and the spherically averaged one (Eq. 6) choosing the parameters as given in Eq.(11) for the full potential and in Eq. (12) for the averaged one, for which critical temperatures T c and critical densities ρ c coincide with their experimental counterparts in both cases.Fig. 6 shows that along the vapor branch of the vapor-liquid coexistence curve the av-eraged potential slightly underestimates the experimental vapor densities, while the fullpotential slightly overestimate them. However, these deviations are of the same order inboth cases, and hardly visible (on the scale of Fig. 6) anyway. Recalling also the fact thatthe coexistence curves (and other data extracted from the simulation) still may suffer fromsystematic effects (residual finite size effects) and statistical errors, see Sec. II, of the order ofup to 0.5%, one should not pay too-much attention to the residual differences. We concludethat both models describe the vapor branch of the coexistence curve equally well, over thestudied range of temperatures (which extends from T c down to about T = 250 K).However, for the liquid branch of the coexistence curve the model with the averagedinteraction performs distinctly better. Of course, there is no physical reason known tous why this should be the case. We believe that this more accurate description of theisotropically averaged model is only due to a fortunate compensation of errors.With respect to the vapor pressure at phase coexistence (Fig. 7), we see, however, thatat low temperature (250 K ≤ T ≤
280 K) the model with the full quadrupolar interactionperforms slightly better than the isotropically averaged model. Near the critical point,however, the isotropic quadrupolar interaction performs slightly better, since it predictsthe critical pressure a bit more accurately. Thus, we conclude that the vapor pressure at1coexistence is predicted by both models about equally well.Fig. 8 finally compares both model predictions with the data [35] for the surface ten-sion between both phases. In this case there is a clear preference for the model with theisotropically averaged quadrupolar interaction. Taking the results from Figs. 6-7 together,we conclude that for a description of phase coexistence the model Eq. (6) is clearly thebetter ”effective” model. Also for a supercritical isobar (Fig. 10 presents an example) themodel Eqs. (4,5,11) does not have an advantage. The comparisons presented in this sectionthus fully justify the use of Eq. (6) for practical purposes.An additional interesting test now concerns the temperature dependence of the densitythat results when we compare NVT simulations for the averaged potential with NpT simu-lations for the full potential (similarly to what was done in Fig. 2), choosing the parametersof Eq. (11) and (12). We have also included a comparison with two analytical approaches,namely an integral equation/density functional theory (IE/DF), see the following section,and perturbation theory combined with mean spherical approximation (MSA), see Ref. [27]for a description of this method in the current context. Both approaches agree with ourresults very well. Fig. 9(a) is the counterpart of Fig. 2(a); again it is seen that the pressureat ρ = 0 .
733 g/cm for the full model is in very good agreement with the correspondingexperimental data and averaged potential. However in this case the full model is superiorwith respect to the averaged model. Fig. 9(b) shows that indeed the NpT results for thefull potential now converge rapidly to a somewhat higher density (near ρ ≈ .
75 g/cm )as the temperature is raised from the critical region to higher temperatures. Of course,as expected, it is not possible to fit the critical region (as done in Fig. 6-9) and the hightemperature region (as done in Figs.2-5) simultaneously. V. INTEGRAL EQUATIONS WITH REFERENCE FUNCTIONALS
In this section we summarize a combined integral equation/density functional methodto calculate equations of state. A novel approach to avoid unphysical no–solution domainsnear the critical point is outlined and data for the equation of state of the averaged modeldefined by Eq. (6) are compared with the simulation results.The pair correlation function h ( r ) and the direct correlation function (of second order) c ( r ), which contain all thermodynamic information of a given homogeneous model system of2density ρ and interacting with an isotropic pair potential U ij ( r ), fulfill the following generalrelations [44, 45, 46]: h ( r ) − c ( r ) = ρ Z d r ′ h ( | r − r ′ | ) c ( r ′ ) , (15)ln[1 + h ( r )] − βU ij ( r ) = h ( r ) − c ( r ) − b ( r ) . (16)The first relation is known as the Ornstein–Zernike equation, while the second is the generalclosure relation in terms of a yet unknown bridge function b ( r ). For a solution, it is necessaryto specify the bridge function b in terms of h and c . Thermodynamics is obtained eitherthrough the virial route [47], giving the pressure p , βpρ = 1 − π ρ β Z ∞ dr r (1 + h ( r )) dU ij ( r ) dr , (17)or through the compressibility route [47] β ∂p∂ρ = 1 − πρ Z ∞ dr r c ( r ) . (18)Both routes are not necessarily identical for a given approximation for the bridge function.For one–component systems, advanced methods exist which yield good agreement withsimulations for the pressure in the whole ρ − T plane and for the whole coexistence curve, e.g.the SCOZA (self–consistent Ornstein–Zernike approximation) [48] or the HRT (hierarchicalreference theory) [49]. A drawback of these methods is their “non–locality” in the ρ − T plane, i.e. for a SCOZA solution a partial differential equation has to be solved on thisplane and for a HRT solution a renormalization flow equation has to be solved on a specifiedisotherm.Therefore we treat our problem within a formalism which is close to the referencehypernetted-chain (RHNC) method. In its original version [50] (developed for repulsivecore fluids), the bridge function b was taken from a reference hard sphere system withsuitably optimized reference packing fraction (or hard sphere diameter). Since the RHNCclosure equation (16) can be derived from an approximate bulk free energy functional, aclosed expression for the chemical potential µ also exists. Near the critical point, however,there exists a no–solution domain (extending into the supercritical region T > T c ) whereno real solution to the RHNC closure can be found. As shown below, this problem can beeliminated by a method (FHNC) where a generating functional for the bridge function isadopted from a suitable reference system.3 A. Bridge functions from a reference functional
The FHNC generalization of RHNC within the framework of density functional theorywas proposed in Ref. [51]. By making a Taylor expansion of the free energy functionalaround bulk densities the general closure equation is derived and the bridge function in thebulk is determined via a density functional for a suitably chosen reference system of hardspheres. To be more explicit, the excess free energy functional F ex [ ρ ( r )] = F HNC [ ρ ( r )] + F B [ ρ ( r )] (19)is split into a part ( F HNC ) which generates the hypernetted-chain (HNC) closure ( b = 0)and a remainder ( F B ) which generates a non-zero bridge function. The HNC part F HNC isa Taylor expansion of the excess free energy up to second order in the deviations from thebulk density ρ , ∆ ρ ( r ) = ρ ( r ) − ρ : F HNC = F ( ρ ) + µ ex Z d r ∆ ρ ( r ) − β Z d r Z d r ′ c ( | r − r ′ | ) ∆ ρ ( r ) ∆ ρ ( r ′ ) . (20)There the defining relations for the excess chemical potential µ ex ( ρ ) and the direct correla-tion function c have been used: δ F ex δρ ( r ) (cid:12)(cid:12)(cid:12)(cid:12) ρ ( r )= ρ = µ ex ( ρ ) , β δ F ex δρ ( r ) δρ ( r ′ ) (cid:12)(cid:12)(cid:12)(cid:12) ρ ( r )= ρ ( r ′ )= ρ = − c ( | r − r ′ | ; ρ ) . (21)The general closure equation follows by employing the test particle method: the grandpotential Ω[ ρ ( r )] = F id [ ρ ( r )] + F ex [ ρ ( r )] − Z d r ( µ − V ( r )) ρ ( r ) (22) (cid:18) β F id [ ρ ( r )] = Z d r ρ ( r )(ln[Λ ρ ( r )] − (cid:19) is minimized in the presence of a fixed test particle of the same species which exerts theexternal potential V ( r ) ≡ U ij ( r ) onto the fluid particles [51, 52] (Λ is the thermal de-Brogliewavelength). The precise form of the closure equation (Eq. 16), is recovered upon thefollowing identifications: h ( r ) = ∆ ρ ( r ) ρ , b ( r ) = β δ F B δρ ( r ) (cid:12)(cid:12)(cid:12)(cid:12) ρ ( r )= ρ ( h ( r )+1) . (23)In general, the excess free energy functional beyond second order, the bridge functional F B ,is not known. Therefore the key step of the present method is to approximate F B by a4density functional for a reference system in the following manner: F B [ ρ ] ≈ F B , ref [ ρ ] = F ref [ ρ ] − F HNC , ref [ ρ ] , (24)where the second order HNC contribution is defined as in Eq. (20) with the replacements F ( ρ ) → F ref ( ρ ), µ ex → µ ex , ref and c → c ref . For fluids with repulsive cores, the referencefunctionals of choice are hard–sphere functionals which are known with high accuracy (suchas e.g. the ones in Refs. [53, 54]). In such a manner the system of equations (15) and (16) isclosed and amenable to numerical treatment. According to Ref. [52] the optimal referencehard sphere packing fraction η ref = ( π/ ρ d is determined through the (local) condition ∂∂d ref (cid:0) F B , ref [ ρ ( h ( r ) + 1); d ref ] − F B , ref [ ρ ( h ref ( r ) + 1); d ref ] (cid:1) ! = 0 , (25)which corresponds to extremizing the free energy difference between the fluid and the refer-ence system with respect to the reference hard sphere diameter d ref . The chemical potentialof the fluid can also be expressed as a functional of h locally in the ρ − T plane [52]. Thusa coexistence curve for a given fluid can be determined straightforwardly by the equality of p and µ on the fluid and the gas side, respectively, and no thermodynamic integrations arenecessary. B. The critical region
Similarly to RHNC and HNC, the FHNC method outlined above, together with theoptimization criterion for η ref , Eq. (25), exhibits a no–solution domain in the ρ − T planewhich stretches into the supercritical region ( T > T c ). This can be attributed to a failure ofthe optimization criterion in the critical region which assigns a wrong long–range behavior tothe direct correlation function c . Consider the asymptotic expansion of the closure, Eq. (16),where h , c and b are small: − h ( r )2 + h ( r )3 + · · · = − c ( r ) − b ( r ) ( r → ∞ ) . (26)(Here we assume that the potential U ij ( r ) is cut off.) In HNC (RHNC) the bridge function b is zero (short–ranged), therefore we find to leading order c ( r ) ≈ h ( r ) /
2. However this isinconsistent with the critical behavior of the correlation functions. In three dimensions, thiscritical behaviour is approximately described by the Ornstein–Zernike form [55] h OZ ( r ) → exp( − r/ξ ) r (27)5Through the Ornstein–Zernike equation (15) it follows that in this limit c OZ ( r ) stays shortranged in the sense that its Fourier transform ˜ c OZ ( q ) = c + c q + . . . permits an expansionaround q = 0. In Eq. (27), ξ is the correlation length which goes to infinity upon approachingthe critical point and is related to ˜ c OZ through ξ = − ρ c / (1 − ρ c ). Assuming h → h OZ ,the asymptotic HNC (RHNC) closure c ≈ h / c staying short–ranged upon approaching the critical point, i.e. it cannot be a solution of theOrnstein–Zernike equation. On the other hand, within FHNC the bridge function b itselfdepends on h as b ( r ) ≈ ( η ref ) β ∂ µ ex , ref ∂ ( η ref ) h ( r ) + O ( h ) ( r → ∞ ) . (28)and the asymptotic closure, Eq. (26) reads c ( r ) = 12 (cid:18) − ( η ref ) β ∂ µ ex , ref ∂ ( η ref ) (cid:19) h ( r ) + O ( h ) . (29)Thus we see that upon requiring β ( η ref ) ∂ µ ex , ref ∂ ( η ref ) ! = 1 the closure is consistent with h → h OZ and c staying short–ranged. This condition is fulfilled for η ref ≈ .
13 and in the criticalregion, this condition on the reference system packing fraction replaces Eq. (25). Incidentally,this condition is consistent with the intuition that the reference hard sphere diameter d ref isroughly equal to the Lennard–Jones diameter σ for densities close to the critical density. Wechecked for various supercritical isotherms that the modified optimization criterion indeedremoves the no–solution domains. (For a different approach to this problem within FHNC,see Ref. [56].) C. Numerical results
Numerical data for the coexistence curve (virial route) of the averaged model (Eq. 6)are given in Fig. 3 ( q A = 0 . q A = 0 . T c where a noticeable shift of the coexisting gas densities towards higher values can beobserved. This is also the reason why the pressure on the coexistence curve is somewhatlarger than the pressure determined by the simulations (Fig. 7). Analyzing the behavior ofthe solutions in the critical region more closely, we find that there (as many other integralequation approaches) the FHNC method suffers from the inconsistency between the virial6and the compressibility route to the equation of state. Indeed, the compressibility route (for q A = 0 . T com c ≈ .
14 and critical density ρ ∗ c ≈ .
36 with abehavior of the correlation function h consistent with the Ornstein–Zernike form, Eq. (27),confirming the applicability of the reasoning in the previous subsection. Therefore, the virialroute coexistence data for temperatures above T com c are not reliable.We remark that the FHNC method is not particularly designed to capture the correctcritical behavior. With respect to the prediction of the coexistence curve and coexistencepressure only, it is not particularly superior in accuracy to the (simpler and faster) per-turbative Mean Spherical Approximation (MSA) whose results have been discussed in [27](for a detail description of the Equation of State used in [27], we refer to [57, 58]). Clearly,in this respect renormalization–group based methods such as HRT perform much better.However, supercritical properties of CO are reproduced quite accurately as the comparisonof the p − ρ isotherms ( T = 316 .
36 K) between the experimental data for CO , FHNC andMSA shows (Fig. 11). While the experimental data and the FHNC results are almost ontop of each other, the perturbative MSA results exhibit a van–der–Waals loop due to theunderlying mean–field approximation which results in a too large T c [27]. Additionally weobserve very good agreement with simulations for the supercritical isochore ρ = 0 .
733 g/cm (Fig. 9) and the supercritical isobars p=100 bar and 200 bar (Fig. 10).A problem of the FHNC approach seems to be the accurate prediction of surface tensions.Although the technique can be extended to compute this quantity, the results are much lesssatisfactory, since the simulation results are about 40% lower than the FHNC results, inthe temperature region around T=270 K where the coexistence curve is predicted rathersatisfactorily by FHNC (Fig. 6). A similar problem was noted in a recent comparison ofMonte Carlo and density functional theory results for phase separation in colloid-polymermixtures [59].In concluding this section, we find that the presented FHNC method allows a fast andprecise determination of the equation of state except for the vicinity of the critical point.Within FHNC, the pressure p and the chemical potential µ are obtained through local re-lations in the T − ρ plane. It appears as an advantage that FHNC is straightforwardlygeneralizable to mixtures since the functionals for the reference hard–sphere mixtures areknown. First studies of Lennard–Jones mixtures [60] confirm the accuracy of the approach.Besides the computation of the pair structure in fluids, the connection to density func-7tional theory makes FHNC also a versatile tool to study wetting/drying phenomena [52]and effective depletion potentials in dilute colloidal solutions [61, 62, 63, 64]. VI. CONCLUSION
In the present work, two coarse-grained models of quadrupolar fluids were comparedwith each other (and with both experiment and a theory (IE/DF) combining an integralequation approach with density functional theory). The aim of our works is to develop someunderstanding for which region of parameters such coarse grained models are accurate, andalso to clarify the applicability of the analytic theory. While we use experimental data forcarbon dioxide as a prototype example of a quadrupolar fluid for comparison, our aim is notto provide a chemically realistic modeling description of this substance or any other material.Recalling that there is a lot of interest to use supercritical carbon dioxide as a solvent andfor chemical processing [6, 7, 8], and that other quadrupolar fluids such as benzene also findwidespread applications, there is a need for efficient coarse-grained models of such simplemolecular fluids. (A chemically realistic modeling of systems like polystyrene-carbon dioxidemixtures would be far beyond reach). Due to the fact that critical fluctuations invalidatesimple analytic theories extending the Van-der-Waals approach (see [65] for a discussion),a simulation approach as presented here is well-suited to include a sufficiently accuratedescription of the critical region.We model the quadrupolar fluid by spherical point particles carrying a quadrupole mo-ment (the strength of this quadrupole moment being taken from experiment), so that thetotal interaction between two molecules is the sum of a Lennard-Jones interaction (Eq. 1)and the quadrupole-quadrupole interaction (Eqs. 2, 3). Possible three-body forces are notat all included explicitly, but to some extent implicitly, since the Lennard-Jones parametersof our effective potential are chosen such that the actual critical temperature and density ofthe material (CO in the chosen example) are reproduced. For the sake of computationalefficency, the potential is truncated at a cutoff distance r c and shifted to zero there (Eq. 4).As a second model we chose a closely related one, where the angular dependence of thequadrupolar interaction is averaged perturbatively, Eq. (6). This isotropic potential can aseasily be treated numerically (sec. V) as other simple isotropic pairwise potentials.By construction, the two models have to agree at very high temperatures, but this is not8the region of interest for applications. In the critical region, discrepancies of the order of4% are typically found, in the case of carbon dioxide.However, when we determine the effective Lennard-Jones parameters for both modelssuch that they reproduce the same critical temperature and density (namely the criticalparameters of carbon dioxide in our case), we find that both models give a similarly accuratedescription of the equation of state over a rather wide region of parameters (Fig. 6, 7,10). Considering also the surface tension (Fig. 8), the simpler model with the averagedinteractions is in better agreement with experiments, despite the fact that at high enoughtemperatures, relative deviations of the averaged model from experiment of the order of1-2% can be identified (see Fig. 9 for example). But even for quantities such as isobars atp=100 bar and 200 bar (Fig. 10), about three times the critical pressure, the experimentalresults are very well reproduced over the full density region of interest (from 0.1 g/cm to 1.0 g/cm ). For such data outside of the critical region, our IE/DF theory yields aquantitatively accurate description without any adjustable parameter whatsoever, providedwe use the Lennard-Jones parameters obtained from the Monte Carlo study as an input(fitting the Lennard-Jones parameter to the critical parameters of the analytic theory is notappropriate, of course, since the latter is inaccurate).Of course the fact that the simple isotropic model is even slightly ”better” than themore complicated one, as far as the comparison with experiment goes, must be attributedto some lucky cancellation of errors. In particular, the temperature dependence of theinterface tension (Fig. 8) suggests that the full model might itself be too simple for anoptimal description of the fluid, and it may be necessary to include additional effects likethe non-spherical shape of the molecule. (The model with the angular-dependent interactionis still far from a full description of chemical reality, of course: but the comparison presentedin [27] shows that many very sophisticated atomistic models do not perform better than thecurrent simple isotropic model either). A possible explanation of this could stay in the factthat the physical quadrupolar moment used in coarse grained interactions (Eqs. 4 and 6)could require some effective corrections related to the truncation of our potentials. However,also in view of our need of efficient coarse grained models, the fact that the averaged modelis definitively faster than the model in which angular degrees of freedom are taken intoaccount, leads to the conclusions of the present work, namely we strongly support the useof the averaged models [27, 29].9We hence suggest that the present approach using effective isotropic models forquadrupolar fluids in spite of the angular dependence of the interactions in such systems isuseful and we plan to extend it to binary fluid mixtures as well. ACKNOWLEDGEMENTS
B.M.M. thanks the BASF AG (Ludwigshafen) for financial support, while M.O. wassupported by the Deutsche Forschungsgemeinschaft via the Collaborative Research Cen-tre (Sonderforschungsbereich) SFB-TR6 ”Colloids in External Fields” (project sectionD6-NWG). CPU times was provided by the NIC J¨ulich and the ZDV Mainz. Usefuland stimulating discussion with F.Heilmann, L.G.MacDowell, M.M¨uller and H.Weiss aregratefully acknowledged. [1] M. P. Allen and D. J. Tildesley,
Computer Simulation of Liquids (Clarendon Press, Oxford,1987).[2] D. Frenkel and B. Smit,
Understanding Molecular Simulation: From Algorithms to Applica-tions , 2nd ed. (Academic Press, San Diego, 2002).[3] A. Z. Panagiotopoulos, Mol. Phys. , 813 (1987); Mol. Sim. , 1 (1992).[4] N. B. Wilding, J. Phys.: Condens. Matter , 585 (1997).[5] K. Binder, M. M¨uller, P. Virnau, and L. G. MacDowell, Adv. Polym. Sci. , 1 (2005).[6] E. Kiran and J. M. H. Levelt-Sengers (eds) Supercritical Fluids (Kluwer, Dordrecht, 1994).[7] J. J. McKetta (ed)
Encyclopaedia of Chemical Processes and Design (Marcel Dekker, NewYork, 1996).[8] M. F. Kemmere and Th. Meyer (eds.)
Supercritical Carbon Dioxide in Polymer ReactionEngineering (Wiley-VCH, Weinheim, 2005).[9] I. M. Prausnitz, R. N. Lichtenthaler, and E. G. de Azevedo,
Molecular Thermodynamics ofFluid Phase Equilibria , 2nd ed. (Prentice Hall, Englewood Cliffs, 1986).[10] J. S. Rowlinson and F. L. Swinton,
Liquids and Liquid Mixtures (Butterworths, London, 1982).[11] K. Binder and G. Ciccotti,
Monte Carlo and Molecular Dynamics of Condensed Matter (So-cieta Italiana di Fisica, Bologna, 1996).[12] S. Bock, E. Bich and E. Vogel, Chem. Phys. , 147 (2000). [13] R. Bukowsky, J. Sadlej, B. Jeziorski, P. Jankovski, K. Szalewicz, A. Kucharski, H. L. Williamsand B. M. Rice, J. Chem. Phys. , 3785 (1999).[14] C. S. Murthy, S. Singer and I. R. McDonald, Mol. Phys. 44, 135 (1981).[15] G. C. Maitland, M. Rigby, E. B. Smith, and W. A. Wakeham, Intermolecular forces, theirorigin and determination (Clarendon Press, Oxford, 1981).[16] H. J. B¨ohm, C. Meissner, and R. Ahlrichs, Mol. Phys. , 651 (1984).[17] H. J. B¨ohm and R. Ahlrichs, Mol. Phys. , 445 (1985).[18] S. B. Zhu and G. W. Robinson, Comp. Phys, Commun. , 317 (1989).[19] R. D. Etters and B. Kuchta, J. Phys. Chem. , 4537 (1989).[20] L. C. Geiger, B. M. Ladanyi and M. E. Clapin, J. Chem. Phys. , 4533 (1990).[21] B. J. Palmer and B. C. Garrett, J. Chem. Phys. , 4047 (1993).[22] J. G. Harris and K. H. Yung, J. Phys. Chem. , 12021 (1995).[23] J. Vorholz, V. I. Harismiadis, B. Rumpf, A. Z. Panagiotopoulos and G. Maurer, Fluid PhaseEquilib. , 203 (2000).[24] J. Stoll, J. Vrabec, H. Hasse and J. Fischer, Fluid Phase Equilib. , 339 (2001).[25] J. Vrabec, J. Stoll and H. Hasse, J. Phys. Chem. B , 12126 (2001).[26] Z. Zhang and Z. Duan, J. Chem. Phys. , 214507 (2005).[27] B. M. Mognetti, L. Yelash, P. Virnau, W. Paul, K. Binder, M. M¨uller, and L. G. MacDowell,J. Chem. Phys. , 104501 (2008).[28] P. Virnau, M. M¨uller, L. G. MacDowell and K. Binder, J. Chem. Phys. , 2169 (2004).[29] E. A. M¨uller and L. D. Gelb, Ind. Eng. Chem. Res. , 4123 (2003); S. Albo and E. A. M¨uller,J. Phys. Chem. B , 1672 (2003).[30] G. Stell, J. C. Rasaiah and H. Narang, Mol. Phys. , 1393 (1974)[31] P. Virnau and M. M¨uller, J. Chem. Phys. , 10925 (2004).[32] K. Binder, Rep. Progr. Phys. , 487 (1997).[33] D. P. Landau and K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics , 2nded. (Cambridge Univ. Press, Cambridge, 2005).[34] C. G. Gray and K. E. Gubbins,
Theory of Molecular Fluids, Vol. I: Fundamentals . ClarendonPress Oxford (1984).[35] NIST website: http://webbook.nist.gov/chemistry/[36] K. Binder and D. P. Landau, Phys. Rev. B , 1477 (1984); C. Borgs and R. Kotecky, J. Stat. Phys. , 79 (1990).[37] A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett. , 2635 (1988).[38] M. E. Fisher, in Critical Phenomena , ed. M. S. Green (Academic Press, London, 1977) p. 1.[39] K. Binder, Z. Physik B , 119 (1981).[40] K. Binder, in Computational Methods in Field Theory , ed. C. B. Lange and H. Gausterer(Springer, Berlin, 1992).[41] K. Binder and E. Luijten, Phys. Rep. , 179 (2001).[42] Y. C. Kim, M. E. Fisher, and E. Luijten, Phys. Rev. Lett. , 065701 (2003); Y. C. Kim andM. E. Fisher, Computer Phys. Commun. , 295 (2005).[43] J. P´erez-Pellitero, P. Ungerer, G. Orkoulas, and A. D. Mackie, J. Chem. Phys. , 054515(2006).[44] T. Morita and K. Hiroike, Prog. Theor. Phys , 1003 (1960).[45] G. Stell, in The Equilibrium Theory of Classical Fluids edited by H. L. Frisch J. L. Lebowitz(Benjamin, New York 1964) p. II-171.[46] J. K. Percus, in
The Equilibrium Theory of Classical Fluids edited by H. L. Frisch andJ. L. Lebowitz (Benjamin, New York 1964) p. II-33.[47] J.–P. Hansen and I. R. McDonald,
Theory of simple liquids (Academic Press, London 2006).[48] D. Pini and G. Stell, Physica A , 270 (2002).[49] A. Parola and L. Reatto Adv. Phys. , 211 (1995).[50] F. Lado, S. M. Foiles, and N. W. Ashcroft, Phys. Rev. A , 2374 (1983).[51] Y. Rosenfeld, J. Chem. Phys. , 8126 (1993).[52] M. Oettel, J. Phys.: Condens. Matter , 429 (2005).[53] Y. Rosenfeld, Phys. Rev. Lett. , 980 (1989).[54] R. Roth, R. Evans, A. Lang, and G. Kahl, J. Phys.: Condens. Matter , 12063 (2002); Y.-X.Yu and J. Wu, J. Chem. Phys. , 10156 (2002).[55] H. E. Stanley, Introduction to Phase Transitions and Critical Phenomena , Oxford UniversityPress (Oxford 1987), chapter 7.[56] S. Amokrane, A. Ayadim, and J. G. Malherbe, Mol. Phys. , 3419 (2006).[57] L. G. MacDowell, M. M¨uller, C. Vega, and K. Binder, J. Chem. Phys. , 419 (2000).[58] L. G. MacDowell, P. Virnau, M. M¨uller, and K. Binder, J. Chem. Phys. , 6360 (2002).[59] R. L. C. Vink and J. Horbach, J. Phys.: Condens. Matter , S3807 (2004). [60] G. Kahl, B. Bildstein, and Y. Rosenfeld, Phys. Rev. E , 5391 (1996).[61] S. Amokrane and J. G. Malherbe, J. Phys.: Condens. Matter , 7199 (2001).[62] M. Oettel, Phys. Rev. E , 041404 (2004).[63] A. Ayadim A, J. G. Malherbe, and S. Amokrane S, J. Chem. Phys. , 234908 (2005).[64] S. Amokrane, A. Ayadim, and J. G. Malherbe, J. Phys. Chem. C , 15982 (2007).[65] A. Kostrowicka Wyczalkowska, J. V. Sengers and M. A. Anisimov, Physica A , 482 (2004). T ∗ = T /ǫ for q F = 0 .
682 usingthe model defined by Eqs. (3,4), for four choices of L , namely L/σ = 6 . , , . L ). The dotted horizontal lines indicatethe theoretical values for the 3d-Ising universality class [41].FIG. 2: (a) Reduced pressure p ∗ ( T ) (in units of the parameters of Eq. 8) plotted vs.reduced temperature T ∗ for the reduced density ρ ∗ = 0 . × ). (b) Reduced density ρ ∗ plotted versus reduced temperature,when one takes the pressure from part (a), as input for an NpT simulation, using the fullpotential, Eqs. (3), (4), with the parameters chosen as given in Eq. (8) (see ⋄ ). Note thatthe statistical errors are estimated not to exceed the size of the symbols.FIG. 3: Vapor-liquid coexistence curve in the ( T ∗ , ρ ∗ ) plane as predicted by Eq. (6)(full line), using the parameters as quoted in Eq. (8), and as predicted by Eqs. (3), (4)(broken line), using the corresponding parameters (Eqs. 7, 8) implying exact agreementbetween both models in the limit T → ∞ . For the averaged model, we also report resultsof the integral equation/density functional theory described in Sec. V (see ∗ ). The relativeaccuracy of the curves representing simulation results in this figure and in the following isestimated to be better than 0.5%.FIG. 4: Vapor-liquid coexistence curve in the ( p ∗ , T ∗ ) plane, for the same choices asin Fig. 3. For the averaged model, we report also results of the integral equation/densityfunctional theory described in Sec. V (see ∗ ).FIG. 5: Interfacial tension plotted vs. T ∗ , for the two models as specified in Fig. 3.FIG. 6: Same as Fig. 3, but choosing the parameters of Eq. (11) for the model withfull quadrupolar interaction (broken line) and of Eq. (12) for the model with the averagedinteraction (full line). Experimental data from Ref. [35] are included (broken-dotted line).With respect to Fig. 3, critical temperature and density for both models now coincidewith the experimental values. We also include the LJ predictions of Ref. [28] (dotted line).We notice that the spherical and averaged model in this reparametrized plot produces4coexistence densities that are in good agreement. Indeed differences are comparable tothe two models used to describe CO in our previous paper [27] (i.e. q A , c = 0 .
387 and q A , c = 0 . ∗ ).FIG. 7: Coexistence pressure. We report the results for the full model with simula-tion parameters given in Eq. (11) (broken line), the results for the averaged model withsimulation parameters reported in Eq. (12) (full line) and the experimental results [35](broken-dotted line). We also include the LJ predictions of Ref. [28] (dotted line). Weobserve that at low temperatures the full model gives better results with respect to theaveraged model. This is due to the fact that for high densities the orientational part ofthe quadrupolar interaction becomes more important. On the other hand, near the criticalpoint the averaged model performs better than the full model. For the averaged model, wereport also results of the integral equation/density functional theory described in Sec. V(see ∗ ).FIG. 8: Prediction of interface tension. We report the results for the full modelwith simulation parameters given in Eq. (11) (broken line), the results for the averagedmodel with simulation parameters reported in Eq. (12) (full line) and the experimentalresults [35] (broken-dotted line). We also include the LJ predictions of Ref. [28] (dottedline). The averaged model is in perfect agreement with experimental results.FIG. 9: (a) Supercritical isochore ( ρ =0.733 g/cm ) for the averaged model with pa-rameters given in Eq. (11) (full line) and for the full model with parameters given in Eq.(12) (broken line). Experimental results are also included [35] (broken-dotted line). Weobserve a very nice agreement between the experimental results and the full potential. Thesmall discrepancy between the averaged model and the full model at high temperaturescan be understood in the light of the results of Fig. 2. Indeed, due to the fact that theuse of the same ǫ and σ for the averaged and full models (see Eq. 8) produces the sameequilibrium states at high temperature (see Fig. 2), the new choice of parameters Eqs. (11)5(12) produces a systematic discrepancy at high temperature. For the averaged model, wealso report results of the integral equation/density functional theory described in Sec. V ( ∗ )and the results of the perturbative MSA theory described in [27] (gray line). The insertedpicture shows as both the theory are in very nice agreement with the MC results.(b). Similarly as in Fig. 2(b) we report the prediction for the densities for the full model(diamond) obtained in an NpT simulation taking as input the pressures plotted in Fig. 9(a)for the averaged model (full line). The resulting densities agree well with experiments [35].The small deviation at high temperature (between the two models) can be explained byconsiderations discussed for Figs. 2 and 9(a).FIG. 10: Supercritical isobars (p=100 and p=200 bar). We report the results forthe full model with simulation parameters given in Eq. (11), the results for the averagedmodel with simulation parameters reported in Eq. (12) and the experimental results [35].The agreement is very good. For a comparison with two other averaged models we referto Fig. 8 of our previous paper [27]. For the averaged model, we report also results of theintegral equation/density functional theory described in Sec. V.FIG. 11: Comparison between the Integral Equation/Density Functional (IE/DFT)theory ( ∗ ) and an equation of state in the perturbative Mean Spherical Approximation(MSA) described in [27, 57, 58] (full line) for the supercritical isotherm T=316.36 K.IE/DFT, if compared to experiments [35], performs better for intermediate densities 0.3g/cm < ρ < . However for ρ > . the two theory predict almost the sameequilibria states. F I G . : T* U U L=11.3 σ L=6.74 σ L=9 σ L=13.5 σ FIG. 2: T* P * ρ *=0.544, q A,c =0.387 T* ρ * q A,c =0.387q F =0.682 F I G . : ρ * T * q F =0.682 q A,c =0.387 (IE/DF)q
A,c =0.387 F I G . : T* P * q A,c =0.387q F =0.682q A,c =0.387 (IE/DF) F I G . : T* γ * q A,c =0.387q F =0.682 F I G . : Density [g/cm ] T e m pe r a t u r e [ K ] q A,c =0.385q F =0.682Exp. (NIST)q A,c =0.385 (IE/DF)LJ (q
A,c =0) F I G . :
240 260 280 300
Temperature [K] P r e ss u r e [ ba r ] q A,c =0.385Exp. (NIST)q F =0.682q A,c =0.385 (IE/DF)LJ (q
A,c =0) F I G . :
240 260 280 300
Temperature [K] I n t e r f a c e T en s i on [ m N / m ] Exp. (NIST)q F =0.682q A,c =0.385LJ (q
A,c =0) FIG. 9:
400 500 600
Temperature [K] P r e ss u r e [ ba r ] q A,c =0.385 (IE/DF)q
A,c =0.385 (MSA)q
A,c =0.385 q F =0.682Exp. (NIST)
560 600110012001300
300 400 500 600
Temperature [K] D en s i t y [ g / c m ] q A,c =0.385q F =0.682Exp. (NIST) F I G . : Density [g/cm ] T e m pe r a t u r e [ K ] Exp. (NIST)q F =0.682q A,c =0.385q
A,c =0.385 (IE/DF) F I G . : Density [g/cm ] P r e ss u r e [ ba r ]]