A comprehensive comparison between thermodynamic perturbation theory and first-order mean spherical approximation: Based on discrete potentials with hard core

Chemical Physics 493 (2017) 1–11

Contents lists available at ScienceDirect

Chemical Physics journal homepage: www.elsevier.com/locate/chemphys

A comprehensive comparison between thermodynamic perturbation theory and first-order mean spherical approximation: Based on discrete potentials with hard core Shiqi Zhou ⇑, Run Zhou School of Physics and Electronics, Central South University, Changsha, Hunan 410083, China

a r t i c l e

i n f o

Article history: Received 22 February 2017 In final form 24 May 2017 Available online 29 May 2017 Keywords: Radial distribution function Thermodynamic perturbation theory

a b s t r a c t Using the TL (Tang and Lu, 1993) method, Ornstein-Zernike integral equation is solved perturbatively under the mean spherical approximation (MSA) for fluid with potential consisting of a hard sphere plus square-well plus square-shoulder (HS + SW + SS) to obtain first-order analytic expressions of radial distribution function (RDF), second-order direct correlation function, and semi-analytic expressions for common thermodynamic properties. A comprehensive comparison between the first-order MSA and high temperature series expansion (HTSE) to third-, fifth- and seventh-order is performed over a wide parameter range for both a HS + SW and the HS + SW + SS model fluids by using corresponding ‘‘exact” Monte Carlo results as a reference; although the HTSE is carried out up to seventh-order, and not to the first order as the first-order MSA the comparison is considered fair from a calculation complexity perspective. It is found that the performance of the first-order MSA is dramatically model-dependent: as target potentials go from the HS + SW to the HS + SW + SS, (i) there is a dramatic dropping of performance of the firstorder MSA expressions in calculating the thermodynamic properties, especially both the excess internal energy and constant volume excess heat capacity of the HS + SW + SS model cannot be predicted even qualitatively correctly. (ii) One tendency is noticed that the first-order MSA gets more reliable with increasing temperatures in dealing with the pressure, excess Helmholtz free energy, excess enthalpy and excess chemical potential. (iii) Concerning the RDF, the first-order MSA is not as disappointing as it displays in the cases of thermodynamics. (iv) In the case of the HS + SW model, the first-order MSA solution is shown to be quantitatively correct in calculating the pressure and excess chemical potential even if the reduced temperatures are as low as 0.8. On the other hand, the seventh-order HTSE is less modeldependent; in most cases of the HS + SW and the HS + SW + SS models, the seventh-order HTSE improves the fifth- and third-order HTSE in both thermodynamic properties and RDF, and the improvements are very demonstrable in both the excess internal energy and constant volume excess heat capacity; for very limited cases, the seventh-order HTSE improves the fifth-order HTSE only within lower density domain and even shows a bit of inadaptation over higher density domain. Ó 2017 Elsevier B.V. All rights reserved.

1. Introduction Radial distribution function (RDF) g(r) plays two important roles in representing the structure of fluid and relating bulk thermodynamic properties to microscopic molecular interactions, and investigation of g(r) and thermodynamic properties of liquid and solid substances has always been a crucial field [1]. The bulk structural and thermodynamic properties are not only fundamental to phase equilibrium calculations of bulk fluid, but also a basic precondition

⇑ Corresponding author. E-mail address: [email protected] (S. Zhou). http://dx.doi.org/10.1016/j.chemphys.2017.05.018 0301-0104/Ó 2017 Elsevier B.V. All rights reserved.

for constructing the grand potential functional approximation, the most basic element in a classical density functional theory (DFT) approach. Although the bulk Helmholtz free energy and bulk second-order direct correlation function c(r) are well known for a long time to play important roles in constructing the free energy functional approximation [2], in a recent publication [3], one of the present authors illustrates that direct incorporation of the RDF information into the functional approximation improves obviously the validity of the latter. Efficient ways to find the fluid g(r) with reliable accuracy and with less calculation include (i) a socalled Ornstein-Zernike (OZ) integral equation [4], which relates a total correlation function denoted by h(r)=g(r)-1 to the c(r); in

2

S. Zhou, R. Zhou / Chemical Physics 493 (2017) 1–11

this approach, accuracy of the RDF depends on the bridge function approximation used. (ii) Via Percus’ test particle method [5]. Here one solves for the equilibrium density profiles in the presence of an external potential that represents a test bulk particle located at the origin. It can be shown that these profiles, when normalized by the bulk density, are the RDF. In this approach, the accuracy of the RDF depends on the grand potential functional approximation, whose minimization gives the Euler–Lagrange equations used to calculate the inhomogeneous equilibrium density profiles. (iii) A coupling parameter series expansion (CPSE) [6], whose higher-order truncation could be easily achieved using numerical solution of the OZ integral equation. Analytic expression of the RDF is of significance as it can, by help of the OZ integral equation, give analytic expression of the c(r), and can, by help of integrating the molecular potential function weighted by the RDF, give analytic expressions of excess internal energy Uex, excess Helmholtz free energy Fex, and other thermodynamic quantities. Approximate and analytic methods are given enough attention in statistical mechanics community [7]. Analytically solving the OZ integral equation under certain bridge function approximation is reliable and common approach to acquire analytic RDF and c(r) [4]. Tang and Lu (TL) propose a method to solve the OZ integral equation for both pure fluids and mixtures [8] with potential consisting of a hard core repulsion with an arbitrary tail under a Percus-Yevick (PY) approximation or a mean spherical approximation (MSA) by combining the perturbation theory with the application of the Hilbert transform. They find that each perturbation term can be solved analytically through the Hilbert transform with operations in the k-space or s-space, representing the Fourier and Laplace transforms, respectively. It is found that the first-order MSA solution dominates in the thermodynamics and structures of the hard sphere (HS)+Yukawa and HS + square well (HS + SW) model fluids, and is found rather successful to approximate the full MSA solution. Note that these studies would have been much more cumbersome had they been performed by the full MSA solution. In the TL method, there are three of the most encountered functions in the solution—g(r), c(r), and a factorization function, which are internally convoluted. At the initial stage, the first-order MSA solution is not fully completed as the studies were all targeting g (r) and ignoring the two other functions. By applying the Hilbert transform, the first-order factorization and c(r) are analytically obtained [9,10], with emphasis on the MSA for the HS + Yukawa and HS + SW fluids. The TL method is elaborated upon further, and a new approach for analytic Laplace inversion to obtain explicit g(r) expression in a compact and consistent manner is proposed [11]; the proposed approach can yield g(r) values directly in any number of shells corresponding to any r values and provide an extraordinary contrast with analytic efforts reported in literature [12–17], all of which are confined to a limited number of shells and to the system of hard spheres, and whose extension to a large number of shells is prohibitive in practice. The TL method had been applied to several typical model potentials [8,9,11,18,19], including the hard sphere, sticky hard sphere, HS + Yukawa, HS + SW, Lennard-Jones and Kihara; in these works, the authors focus on the g(r) and several of the thermodynamic properties. Due to lack of comprehensive comparison between the first-order analytic expressions and computer simulations, it is still unclear how far and how precise the first-order MSA solution accords with the ‘‘exact” simulation results. Potential models consisting of a HS repulsion at close inter-particle separation accounting for the essential impenetrability of the particles, an attractive well at short distances arising from the polymer-induced depletion attraction, plus a soft repulsion at larger distances (inhibiting the vapor-liquid transition which would take place if only the attractive well were present) can serve to

model the effective interactions in solution of colloidal particles in the presence of nonadsorbing polymers or electrolyte solutions. In some alkali metals [20] and alloys [21], the effective interionic pair potentials have a similar form, to some extent. It has been found that fluids with potential models having this general shape may be in an inhomogeneous state with particles arranged in clusters or stripes [22] at low enough temperatures. These kinds of patterns have been observed [23] experimentally and also have been found [24] in computer simulations. The simplest model potential, which has the above mentioned general shape, is the so-called HS + SW + square shoulder (SS) potential. The HS + SW + SS model reduces to the HS + SW model by adjusting the parameters of the former model. Because of its simplicity, the HS + SW model has served as models of a wide variety of physical systems including, e.g., He, Ne, Ar, H2, CO2, CH4, C2H6, n-pentane, and n-butane, and to capture the essential features of the interactions found in colloidal systems; moreover, simulation study had confirmed the presence in the HS + SW system of the Yang-Yang anomaly expected and experimentally found for asymmetric fluids [25]. Aim of the present work is to apply the TL method to solve the OZ integral equation for a single component fluid interacting through the HS + SW + SS potential. The first-order MSA analytic expressions for g(r), c(r), and semi-analytic expressions for common thermodynamic properties such as compressibility factor Z, Uex, Fex, constant volume excess heat capacity C ex V , excess chemical potential lex , and excess enthalpy Hex are obtained and compared with corresponding simulation data available in literature and results based on the CPSE to third-, fifth-order and a high temperature series expansion (HTSE) to seventh-order. The HS + SW + SS fluid is a particular case of the so-called discrete potential fluids (DPF), which have been extensively used within the context of both simple and complex fluids. The structural and thermodynamic properties of the DPF as a function of the interaction range have been recently studied within the context of the OZ integral equation theory and computer simulations [26]. Structure of the present paper is organized as follows: in Section 2 application of the TL method to the HS + SW + SS model is briefly summarized and the obtained analytic expressions are presented; Comparisons with the simulation results and the CPSE to third-, fifth-order and the HTSE to seventh-order are presented and analyzed in Section 3. Our conclusions are summarized in Section 4. In Appendix A and B, some expressions, which are essential to the calculation of basic formulae in the text, and the solving procedure and final expression for c(r) are recorded. 2. Theoretical method 2.1. Analytic and semi-analytic expressions of the first order MSA solution for the HS+SW+SS model The HS + SW + SS model potential is of following form:

8 1; > > > < e; uðrÞ ¼ > a e; > > : 0;

r=r < 1

1 < r=r < k1

k1 < r=r < k2

;

ð1Þ

r=r > k2

where, r is diameter of the HS, and used as length unit throughout the paper; e is energy parameter, k1 and k2 are potential range parameters, and a measures the relative strength of the SW and SS. By defining u1 ðrÞ and u2 ðrÞ:

u1 ðrÞ ¼

8 > < > :

1;

r=r < 1

0;

r=r > k1

e  ae; 1 < r=r < k1 ;

ð2Þ

3

S. Zhou, R. Zhou / Chemical Physics 493 (2017) 1–11

8 r=r < 1 > < 1; u2 ðrÞ ¼ ae; 1 < r=r < k2 ; > : 0; r=r > k2

ð3Þ

u(r) can be rewritten as

uðrÞ ¼ u1 ðrÞ þ u2 ðrÞ:

ð4Þ

If the tail potential is considered as a perturbation on that of the HS, a series expansion for g(r) and c(r) can be written in a power series of k:

gðrÞ ¼ g 0 ðrÞ þ kg 1 ðrÞ þ k2 g 2 ðrÞ þ . . . ; 2

cðrÞ ¼ c0 ðrÞ þ kc1 ðrÞ þ k c2 ðrÞ þ . . . ;

ð5Þ ð6Þ

where the subscript 0 represents the quantities of the HS. The subscripts 1, 2,. . ., indicate the first-order, the second-order perturbation, respectively. k is a perturbation parameter which may vary from 0 to 1. The MSA is given as:

cðrÞ ¼ buðrÞ;

r P r;

ð7Þ

here and throughout the paper, b ¼ 1=kB T with kB Boltzmann constant and T absolute temperature, reduced temperature is defined as T⁄ = kBT/e. Substituting Eqs. (5) and (6) into the OZ integral equation, the MSA relation yields:

c1 ðrÞ ¼ buðrÞ;

r P r;

2 X 1 X rg 11 ðrÞ ¼ ða þ 1Þð1  gÞ8 be ð1 þ nÞð12gÞn j¼0 n¼0

8 9  a1 ðjÞDð6; n; n þ 2; tj ; r  n  1Þ ðk 1Þt r > > > e 1 j > > > < = þa2 ðjÞEð6; n; n þ 2; tj ; r  n  1Þ  ;   > > a3 ðjÞDð6; n; n þ 2; t j ; r  n  k1 Þ > > > > :  ; þa2 ðjÞEð6; n; n þ 2; t j ; r  n  k1 Þ 2 X 1 X ð1 þ nÞð12gÞn rg 12 ðrÞ ¼ að1  gÞ8 be j¼0 n¼0

8 0 9  a1 ðjÞDð6; n; n þ 2; tj ; r  n  1Þ ðk 1Þt r > > > e 2 j > > > < = 0 þa2 ðjÞEð6; n; n þ 2; tj ; r  n  1Þ  :  0  > > a3 ðjÞDð6; n; n þ 2; t j ; r  n  k2 Þ > > > > :  ; þa02 ðjÞEð6; n; n þ 2; t j ; r  n  k2 Þ

r P r and i P 2:

Z 1 U ex gðxÞuðxÞx2 dx ¼ 2pqb NkT 0 Z 1 ½g 0 ðxÞ þ g 11 ðxÞ þ g 12 ðxÞðu1 ðxÞ þ u2 ðxÞÞx2 dx ¼ 12gb 0

The hard core condition is expanded into

g i ðrÞ ¼ 0;

r < r and i P 1:

ð10Þ

The TL analysis method [8] includes three steps: (i) acquiring explicit expression for g L1 ðsÞ, the Laplace transform of rg1(r); (ii) obtaining g1(r) by inverting g L1 ðsÞ analytically; (iii) the Hilbert transform is applied to solving the two remaining functions, c(r), and the factorization function, which are internally convoluted. As the HS + SW + SS potential can be considered as a summary of two HS + SW potentials of potential range k1 and k2 , respectively, the first-order term g 1 ðrÞ consists of g 11 ðrÞ and g 12 ðrÞ, contributed by u1 ðrÞ and u2 ðrÞ, respectively. As a result, the first-order MSA RDF of the HS + SW + SS model is:

gðrÞ ¼ g 0 ðrÞ þ g 1 ðrÞ ¼ g 0 ðrÞ þ g 11 ðrÞ þ g 12 ðrÞ;

ð11Þ

where the following GMSA expression [11] is used for the HS RDF g 0 ðrÞ: 1 X rg 0 ðrÞ ¼ ð12gÞn Cð1; n þ 1; n þ 1; r  n  1Þ n¼0

1 g ð1  gÞ X ð1 þ nÞð12gÞn Dð6; n; n þ 2; z0 ; r  n  1Þ; þ 2 n¼0 2

ð12Þ where, expression of z0 is given in Appendix A. Both g 11 ðrÞ and g 12 ðrÞ have the same mathematical form, as that of the HS + SW model, there is only a different coefficient between them, which is caused by the different coefficient occurring in integrating u1 ðrÞ and u2 ðrÞ into the HS + SW + SS potential:

Z

¼ 12gbðae þ eÞ

k1

1

Z

k1



g 0 ðxÞx2 dx  12gbðae þ eÞ

g 11 ðxÞx2 dx  12gbðae þ eÞ

1

þ 12gbae

Z

k1

g 12 ðxÞx2 dx

1

Z

k2

g 0 ðxÞx2 dx þ 12gbae

1

ð8Þ ð9Þ

ð14Þ

Once the first-order MSA RDF g(r) is acquired, the excess internal energy U ex can be obtained by integrating the potential energy weighted by g(r):

þ 12gbae ci ðrÞ ¼ 0;

ð13Þ

Z

Z

k2

g 11 ðxÞx2 dx

1 k2

g 12 ðxÞx2 dx;

ð15Þ

1

hereafter, x = r/r, expressions of the six integrals involved can be obtained by further utilizing the Hilbert transform [19], and are given in Appendix A. According to definition of constant volume heat capacity, the excess constant volume heat capacity C ex V is calculated by

      C ex @ U ex @ U ex V ¼ ; ¼  T   @T @T Nk Ne NkT q q

ð16Þ

Substituting the expression of U ex into the above equation, one has,

Z k1 C ex V g 11 ðxÞx2 dx þ 12gbðae þ eÞ ¼ 12gbðae þ eÞ Nk 1 Z k1 Z k2  g 12 ðxÞx2 dx  12gbae g 11 ðxÞx2 dx 1

 12gbae

1

Z

k2

g 12 ðxÞx2 dx:

ð17Þ

1

As to the excess Helmholtz free energy Fex, one can readily obtain it by means of the relation

F ex  F 0ex 1 ¼ b NkT

Z

b

0

U ex db NkT

¼ 12gbðae þ eÞ

Z

k1

g 0 ðxÞx2 dx  6gbðae þ eÞ

1

Z

k1



g 11 ðxÞx2 dx  6gbðae þ eÞ

1

Z

þ 12gbae þ 6gbae

Z 1

Z

k1

g 12 ðxÞx2 dx

1 k2

1 k2

g 0 ðxÞx2 dx þ 6gbae

Z

k2

g 11 ðxÞx2 dx

1

g 12 ðxÞx2 dx;

ð18Þ

4

S. Zhou, R. Zhou / Chemical Physics 493 (2017) 1–11

where F0ex is the HS excess Helmholtz free energy, we calculate it by means of a Carnahan-Starling (CS) equation of state [27]. The compressibility factor Z is calculated by

Z  Z hs

¼g

@

Fex F  0ex NKT

@g

R k1 @½ g 0 ðxÞx2 dx F 0ex ¼ F exNKT  12g2 bðae þ eÞ 1 @g Rk Rk @½ 1 g 11 ðxÞx2 dx @½ 1 g 12 ðxÞx2 dx : 6g2 bðae þ eÞ 1 @ g  6g2 bðae þ eÞ 1 @g R k2 R k @½ g 0 ðxÞx2 dx @½ 2 g 11 ðxÞx2 dx þ12g2 bae 1 @ g þ 6g2 bae 1 @g Rk @½ 2 g 12 ðxÞx2 dx þ6g2 bae 1 @ g ð19Þ

Expressions of the six integral derivatives involved are recorded in Appendix A, and the HS compressibility factor Zhs is calculated by the CS equation of state [27]. Finally, the excess chemical potential and excess enthalpy are acquired by basic thermodynamic equalities:

blex ¼ Z  1 þ F ex =NkT;

ð20Þ

Hex =Ne ¼ U ex =Ne þ ðZ  1ÞT  :

ð21Þ

2.2. Seventh-order high temperature series expansion In the HTSE using hard sphere fluid as reference fluid (HS-HTSE) [28], the reduced excess Helmholtz free energy per particle F ex =NkT can be expressed as an expansion in power series of an inverse of reduced temperature T⁄,

F ex =NkT ¼ F exref =NkT þ

X n¼1

Fig. 1. Excess Helmholtz free energy, pressure and excess chemical potential of HS + SW fluids with two sets of combination of potential range k and reduced temperature T*, as shown in figures. Points: simulation data from literature [33]. Black dashed curves: results based on the CPSE to third-order calculated in literature [33]; Black solid curves: the results based on the present first-order MSA semi-analytic expressions; Red short dashed curves: results based on the present 7th-order HS-HTSE.

/

an ðqÞ ; T n

ð22Þ

Fig. 2. Same as in Fig. 1 except that the results are for the excess internal energy, excess enthalpy, and constant-volume excess heat capacity.

5

S. Zhou, R. Zhou / Chemical Physics 493 (2017) 1–11

where, F exref =NkT is a reduced excess Helmholtz free energy per particle of the reference HS fluid. For the present HS + SW + SS potential, the HS potential uHS ðr; rÞ with a HS diameter r is a natural choice as the reference potential, and correspondingly, F exref =NkT is obtained by integrating the CS equation of state [27]. an is the nthorder perturbation term coefficient of the HS-HTSE, and is related to fluctuations of the energy due to the perturbation potential uper ðrÞ calculated in the canonical ensemble of the reference system. As the sum of the reference potential and the perturbation potential is equal to the target potential uðrÞ considered, uper ðrÞ is of following form,

for calculations of all ai of the HS-HTSE in principle. The CPSE is characterized by the following two equations only involved with numerically available two-particle RDF of the reference fluid and imaginary fluid and derivatives of the latter with respect to a coupling parameter n,

8 0; > > > < e; uper ðrÞ ¼ > a e; > > : 0;

F pern =NkT ¼ ð1=n!Þ2pq

r=r < 1

1 < r=r < k1

k1 < r=r < k2

:

F exZhou =NkT ¼ F exref =NkT þ

X

/ F pern =NkT;

ð24Þ

n¼1

and

Z

 @ ðn1Þ g imag ðr; n q; TÞ drr buper ðrÞ   @nðn1Þ 2

; n¼0

ð23Þ

r=r > k2

The expressions of an are involved with 2n -order correlation functions [28] of the reference fluid, consequently their accurate estimations has not been settled for a long time except that a1 is amenable to accurate calculation as a1 is only concerned with two-order correlation function of the reference fluid, which is accurately known for a long time. The CPSE [6] offers a novel route

ð25Þ where, the imaginary fluid, whose RDF is denoted by g imag ðr; n q; TÞ, interacts with a pair potential uimag ðr; n; rÞ:

uimag ðr; nÞ ¼ uHS ðr; rÞ þ nuper ðrÞ:

ð26Þ

@ 0 g imag ðr; n q; TÞ=@n0 jn¼0 ¼ g imag ðr; 0 q; TÞ is simply RDF of the HS bulk fluid as uimag ðr; nÞ with n ¼ 0 reduces to the HS potential. Evaluating F pern =NkT at T⁄ = 1 gives an [29], and this has been practiced over the past few years [6,29,30]. It may be necessary to recall that to numerically acquire g imag ðr; n q; TÞ at several discere values of n   ðnÞ around n ¼ 0 to calculate the derivatives g imag ðr; n qb ; TÞ , we n¼0

close the OZ IE by approximating the imaginary non-hard sphere fluid bridge function by a hard-sphere Malijevsky-Labik (ML) bridge function [31]. Using numerical derivative method combined with the use of double precision real number, the calculation of

Fig. 3. Excess Helmholtz free energy of the HS + SW + SS fluids with two sets of combination of potential parameters, as summarized in Table 1. Symbols: the simulation data from literature [32]; black dashed lines: 5th-order CPSE calculated in literature [32]; black solid + symbol curves: the results based on the present firstorder MSA semi-analytic expression; red short dashed curves: results based on the present 7th-order HS-HTSE. In the case of the second model, squares, circles, up triangles, down triangles, and diamonds correspond to T* = 0.95, 1.05, 1.12, 1.25, and 1.5 for model II; in the case of the fourth model, squares and circles correspond to T* = 0.95 and 1.2 for model IV, respectively.

Fig. 4. As in Fig. 3 but for the pressure.

6

S. Zhou, R. Zhou / Chemical Physics 493 (2017) 1–11

 @ n g imag ðr; n q; TÞ=@nn n¼0 can be carried out to sixth-order, and the HS-HTSE can be thus carried out to seventh-order. By expanding g imag ðr; n qÞ around n ¼ 0 into Taylor series, one has

g imag ðr; n qÞ ¼ g imag ðr; 0 qÞ þ

X

/

 @ n g imag ðr; n qÞ=@nn n¼0 n!

n¼1

nn ;

ð27Þ

taking n ¼ 1, and noticing that uimag ðr; 1Þ ¼ uHS ðr; rÞ þ uper ðrÞ ¼ uðrÞ, the HS + SW + SS under the present consideration, one has g imag ðr; 1 qÞ ¼ gðrÞ, the RDF of the HS + SW + SS. Consequently, Eq. (27) reduces to the expression of the HS-HTSE RDF:

gðrÞ ¼ g HS ðr; q; rÞ þ

X n¼1

/

 @ n g imag ðr; n qÞ=@nn n¼0 n!

:

ð28Þ

3. Comparison with simulation results and CPSE results To evaluate the accuracy of the analytic and semi-analytic firstorder MSA expressions for the RDF and thermodynamic properties, and how the accuracy differs from model to model, we perform a comprehensive comparison between the first-order MSA and simulations [32,33] in dealing with the structures and thermodynamics of both the HS + SW + SS model fluids and its a ¼ 0 limit, i.e. the HS + SW fluids (It is noted that the HS + SW + SS model reduces to the HS + SW model when a is equal to zero and k1 is equal to the HS + SW model potential range parameter k); meanwhile we carry out the seventh-order HS-HTSE calculations for the two model

Fig. 5. As in Fig. 3 but for the excess chemical potential.

fluids and the present seventh-order HS-HTSE results combined with those of the CPSE to third- and fifth-order calculated in literature [32,33] are also exhibited to evaluate the relative performance of the first-order MSA and the thermodynamic perturbation theories, and analyze the convergence of the perturbation theory. It is noted that the inequivalence between the CPSE and HTSE under the general condition is proved [34]; but if the HS potential is used as reference, as in the case of the present study, both the CPSE and HTSE are equivalent, and the CPSE can be rewritten in the form of HS-HTSE [35,29]. The comparisons are presented in Figs. 1 and 2 in the text and Figs. 1 and 2 in the supporting information for the HS + SW model thermodynamics, in Figs. 3–7, in the text and in the supporting information, respectively, for the HS + SW + SS model thermodynamics, in Fig. 8 in the text, for the HS + SW + SS model constant-volume excess heat capacity, in Figs. 9 and 10 in the text and in Figs. 8-10 in the supporting information, respectively, for the RDF of the HS + SW + SS model. For the HS + SW model, two values of k and two temperatures are considered; whereas for the HS + SW + SS model, four sets combination of parameters are considered, as listed in Table 1; the reduced temperatures considered are as follows. Model I: T⁄ = 0.76, 0.82, 0.90, 1.00, and 1.20, Model II: T⁄ = 0.95, 1.05, 1.12, 1.25, and 1.50, Model III: T⁄ = 1.05, and 1.30, Model IV: T⁄ = 0.95, and 1.20. Analysis of the comparison appears below. (i) For the HS + SW fluid with a shorter potential range of k ¼ 1:05, the first-order MSA solution is rather reliable even if the reduced temperature is lowered down to 0.8; upon increasing the potential range to k ¼ 1:15, the temperature lower limit rises to T  ¼ 0:9. However, for the two values

Fig. 6. As in Fig. 3 but for the excess internal energy.

S. Zhou, R. Zhou / Chemical Physics 493 (2017) 1–11

Fig. 7. As in Fig. 3 but for the excess enthalpy.

of the potential range parameter, even if the reduced temperature is lowered to 0.5 and 0.68, respectively, the firstorder MSA is still qualitatively correct. Among the six thermodynamic quantities considered, the accuracy of the first-order MSA for the pressure calculation is the highest, whereas the constant volume excess heat capacity is the most difficult quantity to be predicted accurately. Encouragingly, for the HS + SW fluid, the first-order MSA is still qualitatively correct even for the constant volume excess heat capacity. (ii) For the HS + SW + SS fluid, the performance of the first-order MSA solution suffers significant deterioration. More specifically, only the pressure of model III can be predicted quantitatively over entire density and temperature ranges considered; whereas in the cases of model I and IV, the first-order MSA solution only can reliably predict the pressures over entire density range at higher temperatures, and at lower temperatures, the pressure can only be predicted reliably at lower densities. On the other hand, in the case of model II, the first-order MSA solution completely fails to predict the pressure at high densities whether the temperatures are low or high. As for the excess Helmholtz free energy, excess chemical potential, and excess enthalpy, from Figs. 3, 5, and 7 in text and supporting information, respectively, one can conclude that the first-order MSA solution only gets more reliable with increasing temperatures, and moreover, the temperature lower limit, beyond which, the reliability can be assured, is higher than that in the case of pressure. Finally, both the excess internal energy and constant volume excess heat capacity are the hardest cases, which the first-order MSA almost completely fails to predict

7

even qualitatively over the temperature range considered. In other words, the temperature lower limits of the two quantities are well above those of other thermodynamic quantities, especially the pressure. (iii) By comparison, for the case of HS + SW fluid, the HTSE to third and to seventh order are very accurate for the calculation of the six thermodynamic quantities considered. For the case of the relatively longer potential range, i.e. k ¼ 1:15, the seventh-order HTSE always improves the third-order HTSE over lower density range and both have the same accuracy over higher density range; whereas for the case of the relatively shorter potential range, i.e. k ¼ 1:05, the seventhorder HTSE is superior to the third-order HTSE over lower density range, but inferior to the third-order HTSE over the higher density range. In the case of the harder HS + SW + SS fluid, the HTSE still demonstrates a strong adaptability in comparison with the first-order MSA; generally speaking, the HTSE to fifth and to seventh order can rather accurately calculate all the thermodynamic quantities of the HS + SW + SS fluid except that the performance deteriorates a little at lower temperatures and higher densities in the case of the hardest model II. More specifically, in most cases of the HS + SW + SS models considered, the seventh-order HTSE improves the fifth-order HTSE in calculating all of the six thermodynamic properties, and the improvements are very clear in both the excess internal energy and constant volume excess heat capacity; in particular, the fifth-order HTSE has done very badly for the case of Model II heat capacity, and the seventh-order version is truly very necessary for this quantity. Besides, for the calculations of the excess Helmholtz free energy, excess chemical potential, excess internal energy, and excess enthalpy of higher density Model II fluid, even the seventh-order version is not so satisfactory. For very limited cases, the seventh-order HTSE improves the fifth-order HTSE only at lower density domain and even shows a bit of inadaptation within higher density domain. (iv) Concerning the RDF of the HS + SW + SS fluid, by analyzing the Figs. 9-10 in the text and Figs. 8-10 in the supporting information, respectively, one concludes that the first-order MSA is not as disappointing as it displays in the case of thermodynamic calculations. On the whole, the performance difference between the first-order MSA and the HTSE is somewhat reduced, and is mainly concentrated in the region near the discontinuities in the potential function. After the second discontinuity, with increasing of the inter-particle separation, the performance difference decreases continuously. It is worth noting that for model I, the first-order MSA shows a little abnormality that the calculated RDF becomes negative over the inter-particle separation corresponding to the SS domain of the potential function at lower densities and three lower temperatures. It should be noted that this kind of abnormality also occurs in RDF of a HS + SS model fluid, and can be cured by an exponential enhancement of the first-order MSA RDF [36], and it is found that in the region of low densities and low temperatures, where the first-order MSA fails, the exponential-based theory besides being qualitatively correct also provides with a notable quantitative improvement of the theoretical description. Unfortunately, in the case of the HS + SW + SS fluid, we find that the exponential enhancement of the first-order MSA RDF, is even worse than the first-order MSA RDF itself, and consequently, we have not presented the corresponding results. As for the influence of the truncating order of the HTSE on the RDF of the HS + SW + SS fluid, one can observe that the seventh-order truncation improves the fifth-order truncation over the lower and middle density

8

S. Zhou, R. Zhou / Chemical Physics 493 (2017) 1–11

Fig. 8. As in Fig. 3 but for the constant-volume excess heat capacity of the four sets of combination of potential parameters, as summarized in Table 1.

domain, and the difference tends to reduce with the density. This observation is in agreement with that made in the case of the thermodynamic properties, and this is understandable as the accuracy of the HTSE depends on the validity of the expansion coefficients ai whose estimation by the CPSE is   ðnÞ directly connected with the derivatives g imag ðr; n qb ; TÞ ; n¼0   ðnÞ it is noted that one sum of all orders g imag ðr; n qb ; TÞ n¼0

divided by one appropriate factorial is the RDF in the HTSE. It is observed that the improvement of the seventh-order HTSE on the first-order MSA covers the inter-particle separations corresponding to the SW and SS domains of the potential function, and after the two domains, the improvement continues over an ignorable domain. It may be interesting to have a discussion about the self consistency of the perturbation theory. In liquid distribution function theory, ‘‘self consistency” is normally taken to mean an equivalence between thermodynamic properties from different routes based on structural information of the liquid. As an example, if the pressures computed from the virial and the compressibility routes are equal, the pressure self consistency is considered to be satisfied. The self consistency problem is usually used to improve the distribution function theory. For the computation of thermodynamic properties, in the present perturbation theory, ‘‘self consistency” issue does not arise; this is because one first gains the

Helmholtz free energy, one characteristic function using temperature and volume as natural variables, and all of other thermodynamic properties are obtained by differentiating the Helmholtz free energy. As a result, the thermodynamic ‘‘self consistency” of any thermodynamic perturbation theories has already been fully met. However, the structural functions are also available in the present perturbation theory, the ‘‘self consistency” issue of the structural function in principle does arise, such as the pressure self consistency problem. In the present perturbation theory, the structural function ‘‘self consistency” issue and its role in improving the present perturbation theory are pretty interesting and complex problem, and deserve in-depth studies; we will report the relevant results in a separate publication [37]. Here, we only can point out that the structural function ‘‘self consistency” problem is not met in the present rdf calculations because the pressure self consistency is not forced to be met in solving the imaginary fluid OZ equation. 4. Concluding remarks To summarize, in this paper the OZ integral equation closed with the MSA are analytically and perturbatively solved for the HS + SW + SS fluid. The first-order MSA analytic expressions are obtained for the RDF g(r), second-order DCF c(r), and semi-analytic expressions for thermodynamic properties are derived. The relative performance of the first-order MSA for the HS + SW and HS

9

S. Zhou, R. Zhou / Chemical Physics 493 (2017) 1–11

Fig. 9. Radial distribution function for model II at several reduced densities and two reduced temperatures. Symbols: the simulation data from literature [32]; black dashed lines: 5th-order CPSE calculated in literature [32]; black solid + symbol curves: the results based on the present first-order MSA analytic expression; red short dashed curves: results based on the present 7th-order HS-HTSE. Curves (from the bottom up successively corresponding to q* = 0.1,0.3,0.5,0.7,0.9, respectively) corresponding to higher densities are shifted by an appropriate distance along with vertical axis for clarity.

+ SW + SS fluids in comparison with the HTSE to 7th, 5th and 3rd order is analyzed by one comprehensive comparison between the these theories and ‘‘exact” Monte Carlo simulation. The comparison covers a wide range of bulk properties from six basic thermodynamic properties to RDF, and the bulk states considered are distributed in various domains of the fluid phase space. The conclusions reached and relevant discussions are summarized as follows: (i) Although it is shown that the results of the first-order MSA RDF for the HS + Yukawa and HS + SW fluids are found rather successful to approximate the full MSA solution, and moreover, in the case of the HS + SW model, the present study indicates that the first-order MSA solution is quantitatively correct in calculating the pressure and excess chemical potential even if the reduced temperatures are as low as 0.8, the first-order MSA is model-dependent and sensitive. At non-high temperatures of the HS + SW + SS fluid, the observed significant deviation of the first-order MSA thermodynamic properties from simulation data indicates that the higher order terms in the perturbation series may not be neglected in the case of the HS + SW + SS model. In principle, the higher order terms also can be obtained analytically; however, the algebra in the first-order MSA expressions is rather tedious and frightening, and this may stop such attempt.

Fig. 10. As in Fig. 9 but for models I, III, and IV and that the uppermost curve of the bottom sub-graph is for q* = 0.8 instead of 0.9.

Table 1 Potential model parameters values considered. Model

k1

k2

a

1 2 3 4

1.2 1.2 1.1 1.15

1.35 1.6 1.3 1.4

1 1 0.6 0.85

(ii) Concerning the RDF, the first-order MSA is not as disappointing as it displays in the case of thermodynamic calculations. However, it is noted that for few cases, the first-order MSA shows a little abnormality that the first-order MSA RDF becomes negative over the inter-particle separations corresponding to the SS domain of the potential function at low

10

S. Zhou, R. Zhou / Chemical Physics 493 (2017) 1–11

Fig. 11. Excess Helmholtz free energy of the HS + SW fluid with potential range of k = 1.04 and at one reasonably low temperature of T* = 0.37. Symbols are for computer simulation results [40] and solid lines with different colors are for the HTSE truncated at different orders as marked in the figure.

densities and three lower temperatures; moreover, the exponential enhancement [36] of the first-order MSA RDF, which is a very effective trick in curing the abnormality associated with the HS + SS model fluid, is not right at all for the HS + SW + SS potential. (iii) Due to lack of simulation data for the second-order DCF c(r) of the two models, we have not tested the validity of the first-order MSA c(r). In one recent publication [38], one of the present authors illustrates that scaling an approximate c(r) by the bulk pressure can enhance distinctly the performance of the approximate c(r) in formulating accurate density functional approximation. Consequently, by using the scaling method, even an approximate c(r), just as the present Eq. (B20), may be of great value in constructing the density functional approximation. (iv) The higher-order HTSE theories are shown to be surprisingly accurate for the HS + SW fluid, and the accuracy can be improved continuously with the truncation order changing from third-order, fifth-order to seventh-order; thus, the HTSE shows a good convergence. One remarkable feature of the HTSE is that its seventh-order truncation predicts the constant-volume excess heat capacity (the hardest quantity to deal with by most liquid state theories) rather well in the majority of cases. The model-sensitiveness of the HTSE is much weaker than that shown by the first-order MSA, but is never completely gone, and this is embodied in the Model II, for which the HTSE predictions at higher density domain are generally unsatisfactory. On the other hand, The present HTSE shows an undesirable feature. For a certain potential parameter combination, i.e. the model I, the convergence of the HTSE depends on the bulk density; in lower density domain, the convergence seems to be guaranteed, whereas in higher density domain, there is sign of divergence. To summarize, the HTSE may not be a convergence series, but an asymptotic series. The asymptotic behavior depends on many factors, such as density, model potential under consideration, and specific combination of the model potential parameters, as analyzed above. However, the most important factor determining the asymptotic behavior of the HTSE may be the temperature. To illustrate this point, we calculate the F ex =NkT for the HS + SW fluid with a shallow well width and at one reasonably low temperature using the HTSE truncated at different orders (see Fig. 11). It is seen that although at

higher temperatures the HTSE shows overall convergence behavior even up to seventh order (as analyzed above), the most optimal truncating order for the low temperature situation is third-order in the particular case considered; this fully demonstrates that the HTSE is an asymptotic series. The continuing model-sensitivity of the HTSE in the case of the HS + SW + SS potential indicates that there’s room for improvement in the estimation of the coefficients an by the CPSE; by improving the bridge function approximation used for the calculation of the imaginary fluid RDF, the accuracy of the acquired an is expected to rise. However, in the framework of the first-order MSA, there seemed little room for improvement. (v) It should be pointed out that one probably thinks that the comparison between the first-order MSA and the seventhorder HS-HTSE could be unfair as the HTSE theory is carried out to seventh-order, and not to the first-order as the perturbed-OZ-MSA scheme. However, from a practical point of view, the comparison can be considered to be impartial. First, the expressions for the first-order MSA are modeldependent; as a result, one needs to invest time and invent tricks in acquiring the expressions for different model potentials. On the other hand, the key of the present HS-HTSE is acquiring the coefficients an by solving the CPSE theory at T⁄ = 1; the numerical solution of the CPSE simply consists of numerically solving the OZ integral equation for the imaginary fluid at several discrete values of n around n ¼ 0 and then using numerical derivative method to calculate the   ðnÞ derivatives g imag ðr; n qb ; TÞ . As a result of the numerical n¼0

method, solution of the present HS-HTSE does not differ between different model fluids. To summarize, from the point of availability of the theory, the present HS-HTSE is at least not inferior to the first-order MSA. Second, the expressions of the first-order MSA are frightening, as shown in the Appendix, the expressions of the second-order MSA will be more terrible, let alone those of the seventh-order MSA. Nevertheless, the complexity of the present HS-HTSE will almost not rise with the order of the truncation. Consequently, from the perspective of the complexity of implementing the theory, the present seventh-order HS-HTSE is not at a disadvantage when compared with the first-order MSA. Third, for a given potential parameter combination, the coefficients an in the HS-HTSE is only function of the particle packing fraction g; one only needs to solve the CPSE at T⁄ = 1 for some discrete values of g ranging from approaching zero to a little beyond the solid-liquid transition point, the thermodynamic and structural properties over entire fluid phase space are at hand. Thanks to use of a rapid and stable numerical algorithm for the OZ integral equation such as the LMV acceleration algorithm [39], on a desktop PC approaching ten years old one can obtain ai (with i up to 7) over a reduced density range from zero to 1.0, with step 0.01, with less than 8 min. To summarize, from the point of computational cost, the present seventh-order HS-HTSE is not in an unfavorable position when compared with the analytic first-order MSA as the time-consuming of the two theories are the same order of magnitude. Acknowledgments The authors would like to thank genuinely the anonymous reviewers for their constructive comments which help in deepening the discussions. This project is supported by the National Natural Science Foundation of China (No. 21373274 and 21673299).

S. Zhou, R. Zhou / Chemical Physics 493 (2017) 1–11

Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.chemphys.2017. 05.018. References [1] R.K. Mishra, R. Lalneihpuii, R. Pathak, Chem. Phys. 457 (2015) 13; C. Yigit, J. Heyda, J. Dzubiella, J. Chem. Phys. 143 (2015) 064904; Y. Ebato, T. Miyata, Aip Adv. 6 (2016) 055111; H. Thomsen, M. Bonitz, Phys. Rev. E 91 (2015) 043104; D. Stojiljkovic, J.R. Scepanovic, S.B. Vrhovac, N.M. Svrakic, J. Stat. Mech.-Theory E (2015), P06032; A.P. Lyubartsev, A. Naome, D.P. Vercauteren, A. Laaksonen, J. Chem. Phys. 143 (2015) 243120; A. Sarkar, M.K. Dixit, B.L. Tembe, Chem. Phys. 447 (2015) 76; M. Kurban, O.B. Malcioglu, S. Erkoc, Chem. Phys. 464 (2016) 40; S. Zhou, J. Phys. Chem. Solids 89 (2016) 53; S. Zhou, Phys. Rev. E 92 (2015) 052317; S. Zhou, J. Stat. Mech.-Theory E Paper ID/P11030, 2015.; A. Sengupta, J. Adhikari, Chem. Phys. 469 (2016) 16; E. Keshavarzi, F. Namdari, S.R. Jildani, Chem. Phys. 468 (2016) 15; S. Zhou, G. Liu, AIP Adv. 6 (2016) 045307. [2] D. Henderson, Fundamentals of Inhomogeneous Fluids Dekker, New York, 1992, and reference literature therein; S. Zhou, J. Stat. Mech.-Theory E Paper ID/P11039, 2010; S. Zhou, J. Chem. Phys. 132, 194112, 2010. [3] S. Zhou, J. Chem. Phys. 115 (2001) 2212. [4] J.-P. Hansen, I.R. McDonald, Theory of Simple Liquids, 3rd ed., Academic Press, London, 2006. [5] J.K. Percus, in: H.L. Frisch, A.L. Lebowitz (Eds.), The Equilibrium Theory of Classical Fluids, Benjamin, New York, 1964, p. 113; D. Henderson, Fundamentals of Inhomogeneous Fluids, Marcel Dekker, New York, 1992. [6] S. Zhou, Phys. Rev. E 74 (2006) 031119; S. Zhou, J. Chem. Phys. 125 (2006) 144518. [7] M.Y. Belyakov, E.E. Gorodetskii, V.D. Kulikov, V.P. Voronov, B.A. Grigoriev, Chem. Phys. 445 (2014) 53; T. Nemec, Chem. Phys. 467 (2016) 26; J.M.G. Palanco, L.G. MacDowell, Molec. Phys. 113 (2015) 1076; R. Fantoni, G. Pastore, Molec. Phys. 113 (2015) 2593; Y.V. Kalyuzhnyi, O.A. Vasilyev, P.T. Cummings, J. Chem. Phys. 143 (2015) 044904; M. Klawunn, Phys. Lett. A 380 (2016) 2650; K. Jimenez, J. Reslen, Phys. Lett. A 380 (2016) 2603. [8] Y. Tang, B.C.-Y. Lu, J. Chem. Phys. 99 (1993) 9828. [9] Y. Tang, J. Chem. Phys. 118 (2003) 4140. [10] Y. Tang, J. Chem. Phys. 127 (2007) 164504. [11] Y. Tang, B.C.-Y. Lu, Molec. Phys. 90 (1997) 215. [12] M.S. Wertheim, Phys. Rev. Lett. 10 (1963) 321. [13] G.J. Throop, R.J. Bearman, J. Chem. Phys. 42 (1965) 2408.

[14] [15] [16] [17] [18] [19] [20] [21] [22]

[23] [24] [25]

[26]

[27] [28] [29] [30]

[31] [32] [33] [34] [35] [36] [37]

[38] [39] [40]

11

F. Mandel, R.J. Bearman, M.Y. Bearman, J. Chem. Phys. 52 (1970) 3315. J.W. Perram, Molec. Phys. 30 (1975) 1505. J. Chang, S.I. Sandler, Molec. Phys. 81 (1994) 735. S.P. Goodwin, J.D. Boughey, J.R. Heritage, Molec. Phys. 75 (1992) 917. Y. Tang, B.C.-Y. Lu, J. Chem. Phys. 100 (1994) 3079. Y. Tang, B.C.-Y. Lu, J. Chem. Phys. 100 (1994) 6665. J.-F. Wax, R. Albaki, J.-L. Bretonnet, Phys. Rev. B 62 (2000) 14818. J.A. Moriarty, M. Widom, Phys. Rev. B 56 (1997) 7905. T. Kolokolnikov, H. Sun, D. Uminsky, A.L. Bertozzi, Phys. Rev. E 84 (2011) 015203; A. Hess, M. Pretzl, L. Heymann, A. Fery, N. Aksel, Phys. Rev. E 84 (2011) 031407. R.P. Sear, S.-W. Chung, G. Markovich, W.M. Gelbart, J.R. Heath, Phys. Rev. E 59 (1999) R6255. S.J. Mejía-Rosales, A. Gil-Villegas, B.I. Ivlev, J. Ruiz-García, J. Phys. Condens. Matter 14 (2002) 4795. J.O. Hirschfelder, C.F. Curtiss, R.B. Bird, Molecular Theory of Gases and Liquids, Wiley, New York, 1954; R.V. Gopala, D. Debnath, Colloid Polym. Sci. 268 (1990) 604; G. Foffi, E. Zaccarelli, F. Sciortino, P. Tartaglia, K.A. Dawson, J. Stat. Phys. 100 (2000) 363; G. Orkoulas, M.E. Fisher, A.Z. Panagiotopoulos, Phys. Rev. E 63 (2001) 051507. E. Schöll-Paschinger, A.L. Benavides, R. Castaneda-Priego, J. Chem. Phys. 123 (2005) 234513; I. Guillen-Escamilla, M. Chavez-Paez, R. Castaneda-Priego, J. Phys.: Condens. Matt. 19 (2007) 086224; I. Guillen-Escamilla, E. Schöll-Paschinger, R. Castaneda-Priego, Mol. Phys. 108 (2010) 141; I. Guillen-Escamilla, E. Schöll-Paschinger, R. Castaneda-Priego, Physica A 390 (2011) 3637; M. López de Haro, S.B. Yuste, A. Santos, Mol. Phys. 114 (2016) 2382. N.F. Carnhan, K.E. Starling, J. Chem. Phys. 51 (1969) 635. R.W. Zwanzig, J. Chem. Phys. 22 (1954) 1420. S. Zhou, AIP Adv. 1 (2011) 040703. S. Zhou, J. Chem. Phys. 130 (2009) 014502; S. Zhou, J. Phys. Chem. B 113 (2009) 8635; S. Zhou, J.R. Solana, J. Chem. Phys. 131 (2009) 134106; S. Zhou, J.R. Solana, J. Chem. Phys. 138 (2013) 244115; S. Zhou, J.R. Solana, J. Phys. Chem. B 117 (2013) 9305; S. Zhou, J. Chem. Phys. 139 (2013) 124111; S. Zhou, J.R. Solana, J. Chem. Phys. 141 (2014) 244506. A. Malijevsky, S. Labik, Molec. Phys. 60 (1987) 663. S. Zhou, J.R. Solana, AIP Adv. 3 (2013) 102103. S. Zhou, J.R. Solana, Phys. Rev. E 78 (2008) 021503. S. Zhou, Phys. Rev. E 79 (2009) 011126. A.S.V. Ramana, J. Chem. Phys. 140 (2014) 154106. S.P. Hlushak, P.A. Hlushak, A. Trokhymchuk, J. Chem. Phys. 138 (2013) 164107. S. Zhou, and R. Zhou, An investigation about the structural function ‘‘self consistency” and its role in improving the coupling parameter series expansion perturbation theory, in preparation for Chem. Phys., 2018. S. Zhou, J. Colloid Interface Sci. 298 (2006) 31. S. Labik, A. Malijevsky, P. Vonka, Molec. Phys. 56 (1985) 709. S. Zhou, J.R. Solana, Phys. Chem. Chem. Phys. 11 (2009) 11528.