液態鎵熱力學及結構之研究

(1)

行政院國家科學委員會補助專題研究計畫

期中進

度報告

液態鎵熱力學及結構之研究(1/2)

計畫類別：個別型計畫

計畫編號：NSC 97－2112－M－009－0005－MY2

執行期間： 97 年 8 月 1 日至 99 年 7 月 31 日

計畫主持人：吳天鳴

共同主持人：

計畫參與人員：蔡昆憲

成果報告類型(依經費核定清單規定繳交)：精簡報告

執行單位：交通大學物理研究所

中華民國 97 年 5 月 25 日

(2)

中文摘要

在流體熱力學微擾理論中，多利用硬球流體做參考系統，再加上粒

子間吸力做微擾，以計算流體的自由能。這是過去 Mansoori-Canfield 及

Rasaiah-Stell(MCRS)理論所用的方法。在硬球粒子間所禁止的相空間部分應

對流體的自由能有貢獻；這貢獻對非常軟的粒子流體應是特別明顯。最近，依

據 MCRS 理論對軟粒子流體的熱力學微擾理論已被提出，而計算的結果和直接

從電腦模擬的結果非常符合，數據非常接近。

在一液態鎵的位能模型，兩鎵離子間的作用力包括短程的斥力與長程的

Friedel 振盪位能，而此短程的斥力是非常“軟＂的。此計劃的第一部份，是

利用此修正的 MCRS 理論計算液態鎵的熱力學函數，其一目的是用此液態鎵的

位能模型檢驗此新的軟粒子流體熱力學微擾理論。此一成果已發表於 Journal

of Chemical Physics, Vol. 129, 024503 (2008)。詳細內容請參閱所付的論

文。

(3)

Hard-sphere perturbation theory for a model of liquid Ga

K. H. Tsai1and Ten-Ming Wu2,a兲

1

Department of Electrophysics, National Chiao-Tung University, HsinChu, Taiwan 300, Republic of China

2

Institute of Physics, National Chiao-Tung University, HsinChu, Taiwan 300, Republic of China

共Received 22 April 2008; accepted 30 May 2008; published online 8 July 2008兲

Investigating thermodynamic properties of a model for liquid Ga, we have extended the application of the hard-sphere共HS兲 perturbation theory to an interatomic pair potential that possesses a soft repulsive core and a long-range oscillatory part. The model is interesting for displaying a discontinuous jump on the main-peak position of the radial distribution function at some critical density. At densities less than this critical value, the effective HS diameter of the model, estimated by the variational HS perturbation theory, has a substantial reduction with increasing density. Thus, the density dependence of the packing fraction of the HS reference fluid has an anomalous behavior, with a negative slope, within a density region below the critical density. By adding a correction term originally proposed by Mon to remedy the inherent deficiency of the HS perturbation theory, the extended Mansoori–Canfield/Rasaiah–Stell theory 关J. Chem. Phys. 120, 4844 共2004兲兴 very accurately predicts the Helmholtz free energy and entropy of the model, including an excess entropy anomaly. Almost occurring in the same density region, the excess entropy anomaly is found to be associated with the anomalous packing faction of the HS fluid. © 2008 American Institute of

Physics.关DOI:10.1063/1.2948950兴

I. INTRODUCTION

A general realization about simple liquids is that the liq-uid structure is primarily determined by the repulsive core of particles, while the attractions between particles determine the cohesion of the liquid.1 Since the repulsions and attrac-tions in a liquid play different roles in thermodynamics of the liquid, the two portions of interatomic interactions are indi-vidually approximated with different strategies. The structure of a liquid can be described by a reference fluid with the repulsive core only, whereas the attractions are treated as a perturbation. This forms the central theme of thermodynamic perturbation theory, which has played an important role in the development of statistical mechanism for simple liquids.2,3 Some refinements of the theory progress recently, with an aim of applications in soft condensed systems.4–9

Owing to the benefit of available analytical expressions of the hard-sphere共HS兲 fluid in calculating thermodynamic and structural properties, the reference fluid is further ap-proximated to be a HS one with the HS diameter determined by requiring that the two fluids are equivalent in some ther-modynamic quantity. For example, in the Weeks–Chandler– Andersen共WCA兲 theory,10–13the HS diameter is determined by the equation of equal compressibility between the two fluids so that their free energies differ in the fourth order of the softness parameter, which is a measure for the difference between the repulsive core and the HS potential. Generally, the WCA theory is accurate for the Lennard-Jones共LJ兲 fluids at high densities. However, as the repulsive core is very soft, the softness parameter is so large that the inaccuracy of the WCA theory becomes serious. With a different choice for the

HS diameter, Lado14gave some improvement in the predic-tion of the WCA theory. Alternatively, according to the Gibbs–Bogolubov inequality, Mansoori–Canfield15 and Rasaiah–Stell16 共MCRS兲 developed independently a varia-tional approach, based on a HS reference system and the first-order perturbation. By treating the HS diameter as a variational parameter, this approach provides an upper bound of the Helmholtz free energy. As the repulsive core is softer than that of the LJ potential, the prediction of this variational approach is found to be more accurate than that of the WCA theory. This variational HS perturbation theory has been ap-plied to calculate the structures and thermodynamics of liq-uid metals, which can be described by the effective pair po-tentials with very soft repulsive cores.17–20

Pointed out by Mon21–23in a few years ago, the inaccu-racy of the HS perturbation theory is resulted from an inher-ent deficiency of the theory: A region in the phase space of a realistic fluid is forbidden by a HS fluid due to the singularity of the HS potential. To make up the deficiency of the HS perturbation theory, Mon proposed a correction term that should be added into the variational function to improve the accuracy of the variational approach. By testing the corrected variational function with the inversed-power fluids and a model of liquid Na,24the deficiency of the HS perturbation theory is accounted. Recently, including Mon’s correction term in the two-body approximation, Ben-Amotz and Stell have extended the MCRS theory.25Different from other HS perturbation theories, this extended theory needs two effec-tive HS diameters. The cavity distribution function of a HS fluid, which is in the correction term, is evaluated by the HS-diameter value determined by the Lado-WCA method; other HS-diameter parameters, appearing either explicitly or implicitly in the variational function, are still treated as the variational parameter. Also tested with the inversed-power a兲_{Author to whom correspondence should be addressed. Electronic mail:}

[email protected].

THE JOURNAL OF CHEMICAL PHYSICS 129, 024503共2008兲

(4)

fluids, the extended-MCRS共E-MCRS兲 theory gives more ac-curate predictions than other first-order HS perturbation theories.

In terms of an interatomic pair potential generated by a first-principles pseudopotential theory,26 the structural and dynamic properties of liquid Ga close to the triple point are well described, including the well-known shoulder in static structure factor and the anomaly in the linewidth function of dynamic structure factors at high wavevectors.27 This inter-atomic pair potential is characterized by two features: a ledge-shape repulsive core and the long-range oscillations induced by conduction electrons. Although generated nu-merically, the interatomic pair potential is continuous and smooth, and the behavior of the repulsive core is generally like the two-scale Jagla potential,28which have been used to investigate the waterlike anomalies.29–31 It is interesting to examine the accuracy of the HS perturbation theories men-tioned above with this liquid-Ga model, by comparing their predictions with the simulation results. We also estimate the effective HS diameter of this model fluid by these perturba-tion theories and compare their density variaperturba-tions. We find that the E-MCRS theory gives extraordinarily accurate pre-dictions on the thermodynamic properties of this model from low to high densities, including an excess entropy anomaly at intermediate densities. In the next section, we briefly out-line the MCRS and E-MCRS theories. In Sec. III, our model is described and the conditions of our simulations are given. In Sec. IV, we present the calculated results. Our conclusions are given in Sec. V.

II. THEORETICAL BACKGROUND

The excess Helmholtz free energy Aex _{of a fluid is}

de-fined as the difference between the total free energy of the fluid and that of an ideal gas at the same number density␳ and temperature T. Aex_{is attributed to interactions between}

particles. Here, we give two methods for how Aex_is

calcu-lated. The first method is based on the HS perturbation theory and the second is evaluated by a thermodynamic integration of pressure.

Relative to a reference fluid of HSs with diameter␴HSat

the same temperature and number density, A˜ =Aex_/Nk BT of a

fluid with N particles can be expressed as

A

˜ = A˜_HS+⌬A˜, 共1兲 where A˜HS= AHS

ex_/Nk

BT is the corresponding quantity of the

HS fluid and⌬A˜ is the difference. A˜HS is well described by

the Carnahan–Starling expression,32

A ˜

HS=

␩共4 − 3␩兲

共1 −␩兲2 , 共2兲

where␩=␲␳␴_HS3 /6 is the packing fraction of the HS fluid. In the MCRS theory,15,16⌬A˜ of a fluid with a pair potential␾共r兲 is given as

⌬A˜ ⬇ 2␲␳␤

冕

␴HS

⬁

gHS共r兲␾共r兲r2dr, 共3兲

where gHS共r兲 is the radial distribution function of the HS

fluid and␤=共kBT兲−1. By treating␴HSas a variational

param-eter, a measured value of A˜ is obtained by minimizing a sum of A˜HSand⌬A˜ given in Eqs.共2兲and共3兲, respectively.

In principle, the entire configuration space ⍀ is acces-sible to a realistic fluid. Because of the singularity of a HS potential, part of ⍀ is not available to the HS fluid, which configuration space is denoted as ⍀_HS. This is due to the excluded volume effect of the HS particles. Mon pointed out that due to the difference between ⍀_HS and ⍀ a correction term that should be added into⌬A˜ is formulated as21–23

⌬_␩= 1

Nln

冋

兰⍀HSe

−␤⌽

兰⍀e−␤⌽

册

, 共4兲

where⌽ is the total interaction energy of the fluid of interest. To the lowest order of fluid density, Ben-Amotz and Stell have shown that the correction term can be approxi-mated to be25

⌬_␩= − 2␲␳

冕

0 ␴HS

g共r兲r2dr, 共5兲

where g共r兲 is the radial distribution function of the fluid with

␾共r兲. As commonly used in the HS perturbation theory, g共r兲 in the repulsive-core region is approximated as

g共r兲 = yHS共r兲exp共−␤␾共r兲兲, 共6兲

where y_HS共r兲 is the cavity distribution function of the HS reference fluid. By using this approximation, a new varia-tional formula of ⌬A˜, refereed as the E-MCRS theory, is given as ⌬A˜ = 2␲␳␤

冕

␴HS ⬁ gHS共r兲␾共r兲r2dr − 2␲␳

冕

0 ␴HS yHS共r兲exp共−␤␾共r兲兲r2dr. 共7兲

It has been shown that the E-MCRS theory gives better pre-dictions in thermodynamic properties of the inverse-power fluids than other first-order HS perturbation theories.

How-FIG. 1. 共Color online兲 Interatomic pair potential␾共r兲共solid line兲. The dot-dashed line is the LJ potential with the same␧ and␴. The dashed line is a linear function a1共r/␴兲+a2with a1= −115␧ and a2= 109␧; the dotted line is the function c1exp共−c2r/␴兲共r/␴兲c3, where c1= 4.42⫻106␧, c2= 16.4, and

c3= 1.5. The inset shows the oscillatory part of␾共r兲 at large distances.

(5)

ever, a notice should be given: The HS diameter for describ-ing y_HS共r兲 in the second integral is determined by the Lado-WCA method and may have a different value from that of the variational parameter ␴_HS, which implicitly determines

g_HS共r兲 in the integrand of the first integral and appears

ex-plicitly in the lower and upper limits of the first and second integrals, respectively.

In thermodynamics, A˜ of a fluid at density␳is related to the pressure P共␳

⬘

兲 of the fluid at lower density ␳

⬘

via the following integration:2 A ˜ =

冕

0 ␳

冉

␤P共␳

⬘

兲 ␳

⬘

− 1

冊

d␳

⬘

␳

⬘

. 共8兲

P共␳

⬘

兲 can be obtained via the pressure equation, in which

two inputs are the radial distribution function of the fluid at density␳

⬘

and the first derivative of the pair potential. With the radial distribution functions generated by computer simu-lation, A˜ obtained via the integration in Eq.共8兲is recognized as the simulation results, which have no approximations in thermodynamics.

III. MODEL OF LIQUID Ga AND SIMULATIONS

The interatomic pair potential␾共r兲 used in this paper is generated by a first-principles pseudopotential theory for liq-uid Ga at T = 323 K.26Although the numerical data of␾共r兲 are generated at discrete points, the overall shape of ␾共r兲, shown in Fig. 1, is continuous and smooth. Generally,␾共r兲 has a ledge-shape repulsive core and a long-range oscillatory part with a behavior following the Friedel oscillations.27Two parameters of ␾共r兲 are ␧, the depth of the main attractive well, and ␴, which is the shortest distance at ␾共r兲=0. The corresponding temperature of ␧ is about 47 K and ␴ = 4.04 Å. The repulsive core of␾共r兲, denoted as␾0共r兲, is the

potential for r smaller than r0, the position of the main

at-tractive well, where r0is about 1.07␴. According to the

be-havior of the repulsive core, ␾0共r兲 can be roughly divided

into three sections. As r is larger than 0.9␴,␾0共r兲 with

en-ergy less than 6␧ behaves like the repulsive core of the LJ potential with the same␴and␧.33In the intermediate region between 0.9␴and 0.7␴, the value of␾0共r兲 varies from 6␧ to

30␧ and the shape of␾0共r兲 becomes soft and has a ramplike

behavior, with a reflection point around 0.8␴.34At r less than 0.7␴, ␾0共r兲 increases almost exponentially with decreasing

distance so that the repulsive core changes to be extremely stiff.

With ␾共r兲, we have carried out a series of molecular dynamics simulations at constant NVE ensemble by using a predictor-corrector algorithm35 in a time step of 5.45 fs. In each simulation, the temperature is set to be 6.85␧, corre-sponding to T = 323 K, the box size is fixed at 10.2␴, and the range of␾共r兲 is terminated at half of the box size. All quan-tities given in this paper are in units of␴,␧ and the mass of Ga atom. The number of simulated particles starts at N = 3500. Then, in each simulation for a new thermodynamic state, N is reduced by 100, with the lowest N being 100. At

N = 3500, the simulated system has a number density ␳

= 0.05 Å−3, with a reduced density␳ⴱ⬅␳␴3 equal to 3.305, which is close to that of liquid Ga at T = 323 K and pressure of about 1 bar. At the simulated conditions of N = 3500, the

FIG. 2.共Color online兲 Radial distribution function of the model fluid at the indicated particle number. In each panel, the solid and dashed lines are, respectively, the simulation result and the approximation given by Eq.共6兲 with yHS共r兲 evaluated at the Lado-WCA value of␴HS. The dotted line indi-cates the␴HSvalue obtained by the E-MCRS theory.

FIG. 3. 共Color online兲 Potential of mean force w共r兲 of the model fluid at

N = 3500 共dashed line兲, 2500 共short-dashed line兲, 2000 共dot-dashed line兲,

1000共dot-dot-dashed line兲, and 100 共dotted line兲. The solid line is the in-teratomic pair potential␾共r兲. All potentials are scaled with ␧ and distance is scaled with␴.

024503-3 Perturbation theory for liquid Ga J. Chem. Phys. 129, 024503共2008兲

(6)

static and dynamic structure factors of the system reproduce well the experimentally observed structural and dynamic anomalies of liquid Ga at temperatures close to the melting point共Tm= 303 K兲.

27

The radial distribution functions of our simulated sys-tems at some particle numbers are shown in Fig. 2. At N = 3500, the main peak of the radial distribution function is located at 0.686␴, which is in the extremely stiff region of the repulsive core. As N is decreased, the main peak is low-ered and shifts outwardly, while a shoulder near ␴ appears and then grows. Around N = 2100, corresponding to ␳ⴱ⬇2, the roles of the main peak and the shoulder exchange: the main peak changes to occur near␴ and the shoulder resides in the inner part of the repulsive core. As N keeps on de-creasing, the inner shoulder gradually disappears, associated with a continuously growing main peak near␴as a compen-sation. As N is less than 1000, the number density of the system is rather low and the main peak moves toward the first attractive well, with a final position located at the mini-mum of the attractive well. A similar density variation of the radial distribution function is also observed for the one-scale and two-scale ramp potentials31 introduced by Jagla28 and the core-softened continuous potential with an attractive well.36

Another way to manifest the density effect on the radial distribution function of this model is to compare the potential of mean force, which is defined as w共r兲=−kBT ln g共r兲, with

␾共r兲. w共r兲 is the potential that gives the mean force acting on a particle in the fluid. In concept, the mean force between two neighboring particles in a dense fluid includes the direct force due to the interatomic pair potential and an effective force indirectly intermediated through other particles. In gen-eral, the indirect force depends on the fluid density, so the potential of mean force at high density can be quite different from the interatomic pair potential. For our model fluid, the variation of w共r兲 with N is presented in Fig. 3. At N = 3500,

w共r兲 has a deep attractive well at the first-peak position of g共r兲; this attractive well is attributed to the compactness in

the fluid at such a high density. As N is decreased, the depth

of this attractive well gradually reduces, while a shoulder near␴appears. At some critical N, corresponding to a criti-cal density, the roles of this attractive well and the shoulder in w共r兲 switch with each other. As N is decreased further, the new attractive well in w共r兲 continuously moves out, with its depth becoming deeper first and then attenuated to␧ at final. In the WCA theory,11 the effective␴HS of the model is

determined by the solution of the following integral equation:

冕

0 r₀ yHS共r兲exp共−␤␾0共r兲兲r2dr =

冕

␴HS r₀ yHS共r兲r2dr. 共9兲

This integral equation implies equal compressibility between the reference fluid with the repulsive core␾₀共r兲 and the fluid of HSs with diameter␴HSat the same temperature and

den-sity. In the Lado-WCA method,14 the effective HS diameter is given by the solution of a similar equation, by replacing

y_HS共r兲 in both integrands in Eq. 共9兲 with⳵y_HS共r兲/⳵␴_HS. To solve these equations efficiently, we have replaced yHS共r兲

and ⳵yHS共r兲/⳵␴HS with their analytical expressions in the

hard-fluid 共HF兲 approximation,37,38 which can be simply evaluated in terms of the packing fraction␩ of the HS fluid. By using the discrete data of␾0共r兲 and ␩ evaluated at

the simulated density, we have solved these equations to ob-tain the WCA and Lado-WCA values of HS diameter, and their density variations are shown in Fig.4. For each method, our results are consistent with the general understanding that the effective HS diameter of a fluid is decreased with in-creasing the fluid density.25At each density, the HS diam-eters obtained by the two methods are close, with the WCA value larger. At high densities, both WCA and Lado-WCA

␴HS values are larger than the main-peak position of g共r兲,

which is located in the exponential region of the repulsive core. But, the HS diameters estimated at low densities be-comes smaller than the main-peak position of g共r兲, which has a sudden jump to a larger distance in the repulsive core with a LJ-type behavior around ␳ⴱ= 2. With the results shown in Fig. 2 for some particle numbers, we have exam-ined the approximated g共r兲 given in Eq. 共6兲 for our model, with y_HS共r兲 evaluated in the HF approximation at the Lado-WCA value of HS diameter. For our model at high densities, Eq. 共6兲 is a good approximation of g共r兲 for distances less than the main-peak position, while at low densities the va-lidity of the approximation extends over the main peak of

g共r兲.

IV. THERMODYNAMIC PROPERTIES

Through the pressure equation with g共r兲 generated by our simulations as an input, we have calculated the com-pressibility factor Z⬅␤P共␳兲/␳ of our model and the results are shown by the symbols in Fig.5. Theoretically, Z共␳兲 as a

function of density can be expanded into a series with the leading term to be one.

Z共␳兲 = 1 + b2␳+ b3␳2+ b4␳3+ b5␳4+ ¯ . 共10兲

Here, biis the ith Virial coefficient, which is associated with

␾共r兲 via some diagrams.2

In the second Virial approximation,

Z共␳兲 is approximated to be 1+b2␳and the calculation of b2is FIG. 4. 共Color online兲 Effective HS diameter␴HSas a function of density.

␴HSis estimated by various theories: WCA共triangles兲, Lado-WCA 共dia-monds兲, MCRS 共circles兲, and E-MCRS 共squares兲. The crosses stand for the main-peak position of g共r兲. The symbols of each kind are connected by a line to guide the eye.

(7)

straightforward. Shown in the inset of Fig. 5 is the second Virial approximation of our model, which serves as a check for the accuracy of our calculations at low densities. In order to perform numerical integration in Eq. 共8兲, we fit Z共␳兲 in Fig.5with the series in Eq.共10兲truncated beyond the fourth order of␳ and treat the Virial coefficients as fitting param-eters. To make the fitting accurate enough, the data of Z共␳兲 for␳ⴱ larger and smaller than 1 are fitted separately and the values of the fitting parameters are given in TableI. Although the fitting function has an unnoticed discontinuity at ␳ⴱ= 1, this discontinuity has little effect on the numerical values of

A

˜ at higher densities. By using the fitting function of Z共␳兲, we have done the integration in Eq.共8兲and the calculated A˜ are indicated by the symbols in Fig.7.

In the MCRS and E-MCRS theories, the variational A˜ is minimized with␴HS as the only variational parameter.

Dur-ing the minimization process, gHS共r兲 in the integral in Eqs.

共3兲and共7兲is generated by the Verlet algorithm,39and yHS共r兲

in the second integral in Eq.共7兲is replaced by the analytical expression of the HF approximation37,38with the Lado-WCA

␴HSvalue. Shown in Fig.6are the variational curves of A˜ at

some particle numbers. For each curve, the values of A˜ and the effective HS diameter at the minimum are determined numerically. Since the value of Mon’s correction term de-creases with increasing␴HS, adding this correction term into

the variational function causes the minimum of the varia-tional curve shifting to a lower value of excess free energy and a larger value of effective HS diameter. Thus, for all thermodynamic states, the E-MCRS theory gives a better

prediction on the excess Helmholtz free energy than the MCRS theory does, and the value of the HS diameter ob-tained by the E-MCRS theory is larger than that by the MCRS theory. These results are expected to be true for all kinds of fluid.25

In order to investigate the effect of the approximation in Eq. 共6兲 on the values predicted by the E-MCRS theory, we also evaluate the variational A˜ with Mon’s correction term ⌬␩given in Eq.共5兲 with the simulated g共r兲 used. As shown

in Fig.6, the corresponding variational curve of A˜ generally shifts upward relative to the one obtained by the E-MCRS theory, except that at very low densities the curves of the two methods are found to be almost the same. Therefore, our results lead to the conclusion that the approximation in Eq.

共6兲 makes the prediction of the E-MCRS theory better. In Fig.4, the density variations of the effective HS di-ameters obtained by the MCRS and E-MCRS theories are compared to those obtained by the WCA and Lado-WCA theories. Similar as the results of the WCA and Lado-WCA theories, the ␴HS values estimated by the MCRS and

E-MCRS theories also increase with decreasing density. For

␳ⴱ_{larger than 2, the}_␴

HS values of the MCRS and E-MCRS

theories are limited by the main-peak position of g共r兲, with the density dependence of the E-MCRS value almost along the track of the peak position. For ␳ⴱ less than 2 around which the main-peak position of g共r兲 makes a jump to a larger value, the MCRS and E-MCRS values are apparently released from the restriction due to the main-peak position of

g共r兲 and ascend manifestly as␳ⴱ varies from 2 to 1. For␳ⴱ less than 1, the MCRS value becomes saturated roughly to 0.9␴ as density approaches to zero, but the E-MCRS value continuously increases at very low densities and even passes over the values of the WCA and Lado-WCA theories.

The comparison of the excess Helmholtz free energies estimated by various HS perturbation theories is shown in Fig.7. Apparently, the variational approaches do a better job than the WCA and Lado-WCA theories, which are only good for ␳ⴱless than 1. For all calculated densities, the predicted value of the E-MCRS theory is the one closest to the simu-lation data.

The excess entropy Sex _{of a realistic fluid is defined as}

the difference between the entropy of the fluid and that of an ideal gas at the same density and temperature. In thermody-namics,

Sex_{= −}_共_⳵_Aex_/_⳵_T_兲

N,V. 共11兲

S

˜ is defined as Sex_/Nk

B. In the MCRS theory, the excess

entropy S˜MCRSis equal to −A˜HS, which is the excess entropy

of a HS fluid with the HS diameter of the MCRS value.19In the E-MCRS theory, Mon’s correction term gives an extra contribution to S˜. In order to derive the extra contribution, we use the formula of⌬_␩in Eq.共5兲, rather than the approxi-mated one in Eq. 共7兲. After substituting g共r兲=exp共−␤w共r兲兲

into Eq.共5兲, S˜ in the E-MCRS theory can be expressed as

FIG. 5. 共Color online兲 Compressibility factor, Z=␤P/␳, as a function of density. The circles are calculated by the pressure equation with the simu-lated g共r兲 as an input. The solid line is the fitting function as described in the text with the fitting parameters given in TableI. The inset shows Z at low densities and the second Virial approximation共dashed line兲.

TABLE I. The fitting coefficients of Z共␳兲 for␳ⴱ⬍1 and␳ⴱ⬎1. The data are given in unit of␴.

Z共␳兲 b2 b3 b4 b5

␳ⴱ_⬍1 _0.995 _4.129 _−4.969 _4.670

␳ⴱ_⬎1 _−0.323 _5.228 _−0.523 _0.055

(8)

S ˜ EMCRS= − A˜HS+ 2␲␳

冕

0 ␴HS g共r兲关1 +␤w共r兲兴r2dr. 共12兲

Here, the second term on the right-hand side is due to Mon’s correction term. Then, after using the approximation in Eq.

共6兲 for g共r兲, S ˜ EMCRS= − A˜HS+ 2␲␳

冕

0 ␴HS yHS共r兲exp共−␤␾共r兲兲 ⫻关1 +␤␾共r兲 − ln yHS共r兲兴r2dr. 共13兲

To evaluate S˜EMCRS, the HS diameter is required in three

places in the above equation: the packing fraction in A˜HS, the

upper limit of the integral in the second term, and yHS共r兲,

which is approximated to be the analytical expression of the HF approximation.37,38Care must be taken for that the value of ␴HS in the first and second places should be the one

ob-tained by the minimization of variational A˜ in the E-MCRS theory but the Lado-WCA value of␴_HSshould be used in the third place. Also, we have calculated Eq.共12兲in terms of the simulated g共r兲 and w共r兲, and no significant difference is

found between the results calculated by Eqs.共12兲 and共13兲. Thus, for our model, the two equations give almost identical results.

On the other hand, S˜ can be evaluated by the difference between the excess internal energy U˜ and the A˜ obtained by Eq.共8兲. By using␾共r兲 and the simulated g共r兲, the calculation of U˜ is straightforward2 and the result of S˜ is shown by the symbols in Fig.8. The value of S˜ at␳ⴱ= 3.305 is estimated to be −4.10, which is reasonably comparable to the experimen-tal value of −4.62 of the excess entropy of liquid Ga at a higher temperature.19 At low densities, S˜ decreases almost linearly with density; this decrease arises from the excluded volume effect due to the repulsive core of␾共r兲. Intriguingly, in the intermediate range of density, roughly from␳ⴱ= 1.4 to 1.8, the density curve of S˜ has an anomalous behavior, which has a positive slope of a small value. Beyond this anomalous region, S˜ decreases again in a slower rate with increasing density. Observed also in the two-scale ramp potential,40the

FIG. 6. 共Color online兲 Variational A˜ as a function of

␴HSat the indicated particle number. In each panel, the solid and dashed lines are the E-MCRS and MCRS A˜ , respectively. The dot-dashed line is the variational A˜ with Mon’s correction term given in Eq. 共5兲. For N = 500, the solid and dot-dashed lines are almost the same.

FIG. 7.共Color online兲 The excess Helmholtz free energy A˜=Aex_/Nk

BT

pre-dicted by the E-MCRS共solid line兲, MCRS 共dash line兲, WCA 共dotted line兲, and Lado-WCA共dot-dashed line兲 theories and the results calculated by ther-modynamic integration in Eq.共8兲共symbols兲.

FIG. 8. 共Color online兲 Density dependence of the excess entropy Sex_/Nk

B.

The solid and dashed lines are the predictions of the E-MCRS and MCRS theories, respectively. The dot-dashed line is the excess entropy of the HS fluid with the HS diameter estimated by the E-MCRS theory. The symbols are the simulated results as described in the text.

(9)

anomalous behavior of S˜ is associated with the density anomaly, where the density increases with increasing tem-perature at constant pressure.41

For our model, the comparison between S˜ calculated by simulation and those obtained by the MCRS and E-MCRS theories is shown in Fig.8; the three results of S˜ give similar density behavior. The predictions of the MCRS theory agree with the simulation results only at low densities and produce an excess entropy anomaly in a range shifting a little toward smaller density, but are deviated from the simulation data at densities beyond the anomalous region. The E-MCRS theory has a triumph over the MCRS theory for a perfect agreement with the simulation results in the entire range of density we have calculated. We also present the excess entropy of the HS fluid with the E-MCRS HS diameter, which shows a manifested anomaly with a larger positive slope. The abso-lute value of the excess entropy of this HS fluid is larger than that estimated by the E-MCRS theory. This implies that be-cause part of the phase space is inaccessible by the HS fluid, in which the free volume of each particle is reduced, the difference between the entropy of the HS fluid and that of the ideal gas is, therefore, enhanced. Mon’s correction makes a compensation for the excluded volume effect so that the pre-dictions of the E-MCRS theory for all densities agree well with the simulation results. As indicated in Fig.8, the larger the fluid density is, the more significant the compensation due to Mon’s correction term. Thus, Mon’s correction funda-mentally improves the prediction of the HS variational approach for the excess entropy of a fluid.

To address the reason why the anomaly in S˜ predicted by the MCRS and E-MCRS theories occurs in our model, we have plotted in Fig.9the density dependence of the packing fractions of the reference HS fluids in various perturbation theories. In the WCA or the Lado-WCA theories, the packing faction␩ increases monotonically with the increase in den-sity. At␳ⴱ= 3.3, which is close to that of realistic liquid Ga, the ␩ values estimated by the two theories are over 0.7, which is in the solid phase of a HS system. This indicates that because of the overestimation in the effective HS

diam-eter, the reference HS fluids in the two theories are too much deviated from our model. On the other hand, with the ␴HS

values estimated by the MCRS and E-MCRS theories, the density curve of the packing fraction apparently has a nega-tive slope in the intermediate region, roughly from␳ⴱ= 1.3 to 1.7, which corresponds to the substantial shrinkage in the size of the effective HS with increasing density, as shown in Fig.4. Thus, the packing fraction of the HS fluid is reduced with increasing density in this density region. This anomaly in the packing fraction is the origin for the excess entropy anomaly predicted by the two theories. Beyond this anoma-lous region in packing fraction, the size of the effective HS is limited by the exponentially increasing repulsion in the inner core of the pair potential so that the packing fraction in-creases again with density at a smaller rate. At␳ⴱ= 3.3, the packing fractions estimated by the MCRS and E-MCRS theories are, respectively, 0.473 and 0.5, which are in the fluid-solid coexistence region of the HS system, and the ref-erence HS system is physically reasonable.

V. CONCLUSIONS

In terms of various HS perturbation theories, we have investigated thermodynamic behaviors of a model fluid, with which the structure and dynamic properties of liquid Ga close to the triple point are reproduced.26,27 The interatomic pair potential of the model has two features: a ledge-shape repulsive core and a long-range oscillatory part. Somewhat like the two-scale ramp potential,31,32 the repulsive core of the pair potential varies continuously from an exponential-decay inner core, through a ramplike intermediate region, to a LJ-type outer core. The model has been simulated at con-stant NVE conditions from low to high densities; at the high-est density the simulated system corresponds to realistic liq-uid Ga at T = 323 K. The radial distribution function of the model has an interesting variation with number density: At high density, the main peak of the radial distribution function is located in the inner region of the repulsive core with an exponential decay. At some critical density, the position of the main peak makes a discontinuous jump to the outer core of the LJ-type region. As the density is further decreased, the main peak moves toward the minimum of the first attractive well of the pair potential.

We have used the WCA, Lado-WCA, MCRS, and E-MCRS methods to estimate the effective HS diameter of the model. With the HS diameter estimated by the WCA or Lado-WCA method, the packing fraction of the HS system monotonically increases with density and, therefore, the packing fractions at high densities are over the freezing point of the HS system, which makes the WCA and Lado-WCA theories breakdown for our model at high densities. In the MCRS and E-MCRS theories, due to the structural change resulted from the discontinuous jump in the main-peak posi-tion of the radial distribuposi-tion funcposi-tion, the estimated HS di-ameter is substantially reduced with increasing density, which renders the density curve of the packing fraction hav-ing a negative slope in a region of intermediate density. Be-yond this negative-slope region, the packing faction of the

FIG. 9.共Color online兲 Density dependence of packing fraction␩of the HS fluid with the E-MCRS共solid line兲, MCRS 共dashed line兲, WCA 共dotted line兲, and L-WCA 共dot-dashed line兲␴HSvalues shown in Fig.4.

(10)

HS reference system increases again with density and has a physically reasonable value at high densities.

By including Mon’s correction term to compensate the deficiency of the HS perturbation theory, the E-MCRS theory indeed gives better predictions than the MCRS theory does; the predictions of the E-MERS theory for the excess Hel-mohltz free energy and entropy agree well with the simula-tion results from low to high densities. The excess entropy of our model is found to have an anomalous behavior, which has a density function with a positive slope. This excess entropy anomaly also occurs in the two-scale ramp potential,39with which the waterlike thermodynamic proper-ties are investigated recently.30,31 The excess entropy anomaly in our model can be described by the E-MCRS and MCRS theories, and the anomaly is associated with the sub-stantial reduction of the effective HS diameter estimated by the two variational theories. Our results indicate that the E-MCRS theory is superior to other first-order HS perturba-tion theories for investigating thermodynamics of core-softened potentials.

ACKNOWLEDGMENTS

T.M.W. would like to acknowledge financial support from the National Science Council of Taiwan, R. O. C. under Grant No. NSC 95-2112-M-009-027-MY2.

1_{C. Chandler, J. D. Weeks, and H. C. Andersen,}_Science _{220, 787}_共1983兲. 2_{J. P. Hensen and I. R. McDonald, Theory of Simple Liquids}_共Academic,

New York, 1986兲.

3_{J. A. Barker and D. Henderson,}_{Rev. Mod. Phys.} _{48, 587}_共1976兲. 4_{D. Ben-Amotz and G. Stell,}_{J. Phys. Chem. B} _{108, 6877}_共2004兲. 5_{F. O. Raineri, G. Stell, and D. Ben-Amotz,}_{J. Phys.: Condens. Matter}_16,

S4887共2004兲.

6_{S. Zhou,}_{Phys. Rev. E} _{74, 031119}_共2006兲. 7_{S. Zhou,}_{J. Chem. Phys.} _{125, 144518}_共2006兲. 8_{S. Zhou,}_{J. Phys. Chem. B} _{111, 10736}_共2007兲. 9_{A. B. Adib,}_{Phys. Rev. E} _{75, 061204}_共2007兲.

10_{D. Chandler and J. D. Weeks,}_{Phys. Rev. Lett.} _{25, 149}_共1970兲.

11_{J. D. Weeks, D. Chandler, and H. C. Andersen,}_{J. Chem. Phys.} _{54, 5237} 共1971兲.

12_{H. C. Andersen, J. D. Weeks, and D. Chandler,}_{Phys. Rev. A} _{4, 1597} 共1971兲.

13_{H. C. Andersen, D. Chandler, and J. D. Weeks,}_{Adv. Chem. Phys.} _34, 105共1976兲.

14_{F. Lado,}_{Mol. Phys.} _{52, 871}_共1984兲.

15_{G. A. Mansoori and F. B. Canfield,}_{J. Chem. Phys.} _{51, 4958}_共1969兲. 16_{J. Rasaiah and G. Stell,}_{Mol. Phys.} _{18, 249}_共1970兲.

17_{D. Stroud and N. W. Ashcroft,}_{Phys. Rev. B} _{5, 371}_共1972兲. 18_{J. Hafner,}_{Phys. Rev. A} _{16, 351}_共1977兲.

19_{R. Kumaravadivel and R. Evans,}_{J. Phys. C} _{9, 3877}_共1976兲.

20_{J. Hafner, From Hamiltonians to Phase Diagrams}_{共Springer-Verlag,} Ber-lin, 1987兲.

21_{K. K. Mon,}_{J. Chem. Phys.} _{112, 3245}_共2000兲. 22_{K. K. Mon,}_{J. Chem. Phys.} _{115, 4766}_共2001兲. 23_{K. K. Mon,}_{J. Chem. Phys.} _{116, 9392}_共2002兲. 24_{K. K. Mon,}_{Phys. Rev. E} _{63, 061203}_共2001兲.

25_{D. Ben-Amotz and G. Stell,}_{J. Chem. Phys.} _{120, 4844}_共2004兲. 26_{S. F. Tsay and S. Wang,}_{Phys. Rev. B} _{50, 108}_共1994兲.

27_{K. H. Tsai, T. M. Wu, S. F. Tsay, and T. J. Yang,}_{J. Phys.: Condens.} Matter 19, 205141共2007兲.

28_{E. A. Jagla,}_{J. Chem. Phys.} _{111, 8980}_共1999兲.

29_{J. R. Errington and P. G. Debenedetti,}_Nature_共London兲 _{409, 318}_共2001兲. 30_{Z. Yan, S. V. Buldyrev, N. Giovambattista, and H. E. Stanley,}_{Phys. Rev.}

Lett. 95, 130604共2005兲.

31_{Z. Yan, S. V. Buldyrev, N. Giovambattista, P. G. Debenedetti, and H. E.} Stanley,Phys. Rev. E 73, 051204共2006兲.

32_{N. F. Carnahan and K. E. Starling,}_{J. Chem. Phys.} _{51, 635}_共1969兲. 33_{T. M. Wu, S. F. Tsay, S. L. Chang, and W. J. Ma,} _{Phys. Rev. B} _64,

064204共2001兲.

34_{T. M. Wu, W. J. Ma, S. L. Chang, and S. F. Tsay, Physica B}_{共Amsterdam兲}

316–317, 606共2002兲.

35_{A. Rahman,}_{Phys. Rev.} _{136, A405}_共1964兲.

36_{A. B. de Oliveira, P. A. Netz, T. Colla, and M. C. Barbosa,}_{J. Chem.} Phys. 125, 124503共2006兲.

37_{L. E. S. de Souza and D. Ben-Amotz,}_{Mol. Phys.} _{78, 137}_共1993兲. 38_{D. Ben-Amotz and G. Stell,}_{J. Chem. Phys.} _{119, 10777} _{共2003兲; 120,}

4994共2004兲.

39_{L. Verlet and J. J. Weiss,}_{Phys. Rev. A} _{5, 939}_共1972兲.

40_{J. R. Errington, T. M. Truskett, and J. Mittal,}_{J. Chem. Phys.}_{125, 244502} 共2006兲.

41_{R. M. Lynden-Bell and P. G. Debenedetti,}_{J. Phys. Chem. B} _{109, 6527} 共2005兲.