Structure and thermodynamic properties of Sutherland fluids from computer simulation and the Tang–Lu integral equation theory

Fluid Phase Equilibria 253 (2007) 67–73

Structure and thermodynamic properties of Sutherland fluids from computer simulation and the Tang–Lu integral equation theory A. D´ıez a , J. Largo b , J.R. Solana a,∗ a

Departamento de F´ısica Aplicada, Universidad de Cantabria, 39005 Santander, Spain b Dipartimento di Fisica, Universit` a di Roma La Sapienza, I-00185 Roma, Italy

Received 24 November 2006; received in revised form 15 January 2007; accepted 16 January 2007 Available online 21 January 2007

Abstract Monte Carlo simulations in the NVT ensemble have been performed for fluids of hard spheres plus attractive inverse power tails, to obtain the radial distribution function, the equation of state and the excess energy. The results obtained from the Monte Carlo simulations are compared with the predictions of the Tang–Lu [Y. Tang, B.C.-Y. Lu, Mol. Phys. 90 (1997) 215–224] integral equation theory. It is found that the agreement between theory and simulation is quite good for the radial distribution function and excellent for the equation of state and the excess energy. © 2007 Elsevier B.V. All rights reserved. Keywords: Monte Carlo simulation; Structure; Equation of state; Sutherland potential; Integral equation theory

1. Introduction The Sutherland potential, consists of a hard-sphere repulsive core plus an attractive inverse power tail, that is ⎧ ⎨ ∞, r≤σ  σ γ u(r) = (1) ⎩ − , r>σ r in which σ is the diameter of the particles, − is the maximum potential depth, and γ is a parameter that determines the effective range of the potential. Relatively little attention has been paid to the study of systems with particles interacting by means of a potential model of the form (1) both from computer simulation [1–4] and from theory [4–7]. However, this potential model might be useful because of its resemblance with other intermolecular potentials. Thus, the Sutherland potential resembles the Lennard–Jones potential, for γ = 6, and therefore would be useful for studying simple fluids, and the Mie potential, widely used for modeling complex molecular fluids [8–10], for other values of γ. The Girifalco potential for C60 is also well reproduced by a potential of the form (1) with γ ≈ 12 [11]. Quite recently the Sutherland potential has ∗

Corresponding author. Tel.: +34 942201447; fax: +34 942201402. E-mail address: [email protected] (J.R. Solana).

0378-3812/$ – see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.fluid.2007.01.011

been used to modelize different kind of interactions in complex fluids and mixtures [12–14]. An advantage of a hard-core potential, such as the Sutherland potential, with respect to the Lennard–Jones or Mie potentials, is that it is particularly suitable for the application of the Tang–Lu theory [15,16]. In the present work we will test the performance of this theory against computer simulations for the structural and thermodynamic properties of fluids with potentials of the form (1). To this end, we have performed computer simulations for the radial distribution function (r.d.f.) of this kind of fluids, which are described in the next section. These data together with those for the excess energy and the compressibility factor previously reported [4] are used to test the accuracy of the analytical expression derived from the Tang and Lu [15,16] theory, based on the mean spherical approximation and a perturbative expansion. One of the advantages of this theory over other perturbative theories, such as the Barker–Henderson [17] or the Weeks–Chandler–Abdersen [18] theories, is that it provides analytical expressions for the direct correlation function c(r) for different kinds of intermolecular potentials [19,20]. This function is particularly suitable for the study of inhomogeneous fluids. The derivation of the expression for the r.d.f. from the Tang and Lu theory is outlined in Section 3, and the calculated values of the r.d.f., the excess energy, and the equation of state are compared in Section 4 with computer simulations and discussed.

68

A. D´ıez et al. / Fluid Phase Equilibria 253 (2007) 67–73

• Hard spheres in the generalized mean spherical approximation (GMSA):

2. Computer simulations We have performed NVT Monte Carlo computer simulations for fluids with Sutherland potentials of the form (1) with γ = 6, 12, 18, and 36, for several reduced temperatures T ∗ = kB T/, and reduced densities ρ∗ = ρσ 3 covering most of the density range of the fluid phase. The cut-off distance of the potential was fixed at rc = 3σ. The systems consisted in N = 500 particles. For each state considered, the system was equilibrated for 2 × 104 cycles, each of them consisting of an attempted move per particle. Acceptance ratio was fixed at around 50%. The r.d.f. g(r), the excess energy U E over an ideal gas, and the compressibility factor Z = pV/NkB T were determined from averaging over the next 5 × 104 cycles. The two latter quantities were reported in a previous paper [4]. The values of the r.d.f. that will be used here are listed in Table A.1.

h(r) = h0 (r) + εh1 (r) + ε2 h2 (r) + · · ·

(2)

c(r) = c0 (r) + εc1 (r) + ε2 c2 (r) + · · ·

(3)

where subscript 0 refers to the hard-sphere reference system, and subscripts i refer to the ith perturbative contribution. Combining the preceding expressions with the Ornstein–Zernike equation:  h(r) = c(r) + ρ c(r )h(r − r  ) dr (4) these authors developed a general procedure, within the framework of the mean spherical approximation, to obtain the radial distribution function for fluids with a hard-core potential as a sum of contributions in the form: g(x) =

∞ 

gi (x)

(5)

i=0

where x = r/σ is the reduced distance. The procedure lead to analytical expressions for the zero- and first-order contributions in expression (5) for a number of hard-core potentials. These expressions were reported in Ref. [16] where the interested reader is addressed for further details. Here we will give the explicit expressions of g0 (x) and g1 (x) for several cases that are relevant here for convenience as well as for the fact that in Ref. [16] there are several misprints that are not easily detected and we have corrected here. • Hard spheres in the Percus–Yevick (PY) approximation: xg0 (x) =

∞  n=0

(−12η)n C(1, n + 1, n + 1, x − n − 1)

∞

+

η2 (1 − η)  (1 + n)(−12η)n 2 n=0

× D(6, n, n + 2, z0 , x − n − 1)

(7)

• Hard core Lennard–Jones: The potential function in this case is defined in the form: ⎧ r < σeff ⎨ ∞,  σ 12  σ 6

(8) u(r) = ⎩ 4ε , r > σeff − r r 

σeff =

Tang and Lu [15,16] express the total h(r) and direct c(r) correlation functions as a series expansion in terms of a suitable parameter  in the form:

(−12η)n C(1, n + 1, n + 1, x − n − 1)

n=0

with

3. Tang–Lu integral equation theory for Sutherland fluids

∞ 

xg0 (x) =

σ

[1 − e−βu(r) ] dr

(9)

0

For this potential model we have [16] xg1 (x) = 4βε(1 − η)4  ∞ ((σ/σeff )6 /4!)t 4 − ((σ/σeff )12 /10!)t 10 −t × dt e Q20 (t) 0 ×

∞ 

(1 + n)(−12η)n D(6, n, n + 2, t, x − n − 1)

n=0

(10)

• From the preceding expression one readily obtains that corresponding to the Sutherland potential:  ∞ γ−2 ∞ t /(γ − 2)! −t  xg1 (x) = βε(1 − η)4 dt (1 + n) e Q20 (t) 0 n=0 × (−12η)n D(6, n, n + 2, t, x − n − 1)

(11)

In expressions (6), (7), (10), and (11) n is the number of the shell, x is the reduced distance (x = r/σ in Eqs. (6), (7), and (11), and x = r/σeff in Eq. (10)), and Q0 (t) =

S(t) + 12ηL(t) e−t (1 − η)2 t 3

(12)

S(t) = (1 − η)2 t 3 + 6η(1 − η)t 2 + 18η2 t − 12η(1 + 2η) (13)  η L(t) = 1 + t + 1 + 2η (14) 2 where η = (π/6)ρσ 3 is the packing fraction for a hard-sphere system consisting of particles with diameter σ. Moreover C(n1 , n2 , n3 , x)

(6)

=

n3 2   B(n1 , n2 , n3 , i, α) α=0 i=1

(i − 1)!

xi−1 etα x H(x),

n1 + n2 < 3n3 ;

A. D´ıez et al. / Fluid Phase Equilibria 253 (2007) 67–73

C(n1 , n2 , n3 , x) =

(1 + η)n2 δ(x) + the preceding expression, (1 − η)2n3

n1 + n2 = 3n3

The terms in the sum (22) are related to those in the sum corresponding to the energy

(15)

where δ(x) is the Dirac δ and H(x) is the Heaviside step function,

69

U=

∞ 

Ui

(23)

i=0

n −i−k1 3 −i n3 1 (−1)n3 −i−k1 (n3 − 1 + k2 )!(2n3 − 1 − i − k1 − k2 )! A(n1 , n2 , k1 , α) B(n1 , n2 , n3 , i, α) = 2n3 n3 +k2 2 (1 − η) (t − t ) (tα − tγ )2n3 −i−k1 −k2 k1 !k2 !(n3 − i − k1 − k2 )!((n3 − 1)!) α β k1 =0 k2 =0 (16) n2 

A(n1 , n2 , k1 , α) =

i=max(k1 −n1 ,0)



× 1+

η i 2

which in turn can be obtained from the energy equation

n2 !(i + n1 )! i!(n2 − i)!(i + n1 − k1 )!

(1 + 2η)n2 −i tαn1 +i−k1

UE = 2πρ N (17)

D(n1 , n2 , n3 , z, r) n3 2   (−1)i B(n1 , n2 , n3 , i, α) α=0 i=1

(tα + z)i

⎛

× ⎝e−zr − etα r

i−1  (−1)j (tα + z)j r j j=0

j!

⎞

U2 = 2πρ∗ Nε

(1 + η/2)n2 −zr e H(r) + the preceding expression, (1 − η)2n3 (18)

(19)

Within this scheme, the compressibility factor would be given by Z=

∞ 

Zi

(20)

i=0

where the ith contribution Zi can be obtained from the general relationship

pi V ∂(Fi /NkB T ) Zi ≡ (21) =ρ NkB T ∂ρ T where Fi is the ith term in the sum corresponding to the Helmholtz free energy F=

∞  i=0

Fi



∞

g1(x)u∗ (x)x2 dx

(26)

0

U E = U1 + U2

This completes the expressions of the contributions g0 (x) and g1 (x) in the sum (5) and the complete r.d.f. up to first order would be given by g(x) = g0 (x) + g1 (x)

(24)

respectively, where u∗ (x) = u(x)/N. To this order, the excess energy will be given by

D(n1 , n2 , n3 , z, r)

n1 + n2 = 3n3

g(r)u(r)r 2 dr

0

and

⎠ H(r),

n1 + n2 < 3n3 ;

=

∞

by means of introducing in this equation the expression (5) for g(x). In practice, because we only have available analytical expressions for g0 (x) and g1 (x) we can obtain only the firstand second-order terms in the sum (23). They are given by  ∞ U1 ∗ g0(x)u∗ (x)x2 dx (25) = 2πρ Nε 0

and

=



(22)

(27)

The first- and second-order terms in the sum (22) can be obtained from the relationships F1 = U1

(30)

and F2 = 2U2

(31)

From them, and the general expression (21), we can readily obtain the corresponding terms Z1 and Z2 in the sum (20). To second order, the compressibility factor would be given by Z = Z0 + Z1 + Z2

(32)

where Z0 is the compressibility factor of the reference hardsphere fluid. 4. Results and discussion We have obtained the r.d.f. of Sutherland fluids with γ = 6, 12, 18, and 36 by means of the Tang–Lu procedure using expression (19) with g0 (x) given by the GMSA result (7) and g1 (x) given by Eq. (11). Results are compared in Figs. 1 and 2 for the reduced density ρ∗ = 0.9 and two different temperatures

70

A. D´ıez et al. / Fluid Phase Equilibria 253 (2007) 67–73

Fig. 1. Radial distribution function g(x) of Sutherland fluids with γ = 6, 12, 18, and 36, from down upwards, at reduced density ρ∗ = 0.9 and T ∗ = 0.6 in the first case and T ∗ = 0.5 in the remaining cases. Circles: MC simulations. Curves: Tang–Lu theory with the GMSA solution for g0 (x). For clarity, each curve, and the corresponding set of data, has been displaced upwards by a unit with respect to that immediately below.

with the simulation data from this work. We can see that the theory is in qualitative agreement with simulations, but underestimates the simulation data near contact distance. The situation is similar for other densities and temperatures, although the disagreement near contact increases as temperature and density decrease. This behaviour is due to the fact that for high values of γ and low temperatures the short effective range of the potential causes the particles to stick together, which results in a certain degree of clustering [1]. In spite of this, the agreement between theory for the excess energy U E , as determined from Eqs. (23)–(26), and computer simulation is very good, as shown in Figs. 3–6 , although some deviations appear at very low temperatures for short-ranged potentials, namely γ ≥ 12. A similar picture emerges from Figs. 7–10 for the compressibility factor Z determined from Eqs. (20) and (21), together with (30) and (31). Again the agreement between theory and

Fig. 2. As in Fig. 1 but for the fact that T ∗ = 1.0 for all values of γ.

Fig. 3. Excess internal energy −U E /N of a Sutherland fluid with γ = 6 as a function of the reduced density ρ∗ . Points: simulation data from Ref. [4] for T ∗ = 0.6, 0.8, 1.0, 1.5, 2.0, and 3.0, respectively, from the bottom upwards. Curves: Tang–Lu theory with the GMSA solution for g0 (x). For clarity, each curve, and the corresponding set of data, has been displaced upwards by a unit with respect to that immediately below.

Fig. 4. As in Fig. 3 for γ = 12 and T ∗ = 0.5, 0.7, 1.0, 1.5, 2.0, and 3.0.

Fig. 5. As in Fig. 3 for γ = 18 and T ∗ = 0.5, 0.7, 1.0, 1.5, 2.0, and 3.0.

A. D´ıez et al. / Fluid Phase Equilibria 253 (2007) 67–73

Fig. 6. As in Fig. 3 for γ = 36 and T ∗ = 0.5, 0.7, 1.0, 1.5, 2.0, and 3.0.

Fig. 7. Compressibility factor Z = pV/NkB T of a Sutherland fluid with γ = 6 as a function of the reduced density ρ∗ . Points: simulation data from Ref. [4] for T ∗ = 0.6, 0.8, 1.0, 1.5, 2.0, and 3.0, respectively, from the bottom upwards. Curves: Tang–Lu theory with the GMSA solution for g0 (x).

Fig. 8. As in Fig. 7 for γ = 12 and T ∗ = 0.5, 0.7, 1.0, 1.5, 2.0, and 3.0.

71

Fig. 9. As in Fig. 7 for γ = 18 and T ∗ = 0.5, 0.7, 1.0, 1.5, 2.0, and 3.0.

computer simulations is excellent except for very short-ranged potentials, γ ≥ 18 and very low temperatures, T ∗ ≤ 0.7. It seems worthwhile to compare the performance of the perturbative first-order Tang–Lu theory with the celebrated Barker–Henderson second-order perturbation theory. The application of the second of these theories to Sutherland fluids was analyzed in a previous paper [4]. The interested reader may address to this paper for details. Therefore, we have determined percent deviations D(%) =  the mean absolute sim )|/Asim , where A is either U E or Z, (100/n) ni=1 |(Acalc − A i i i between the calculated quantity by means of these two theories and the corresponding simulation data. From the secondorder Barker–Henderson perturbation theory we have obtained D(%) = 14.76 and 12.15 for the excess energy U E using the macroscopic compressibility (mc) and the local compressibility (lc) approximations, respectively, and D(%) = 4.86 (mc) and 4.67 (lc) for the compressibility factor Z. From the Tang–Lu theory we have obtained D(%) = 11.69 for U E and D(%) = 3.74 for Z. Therefore, the first-order Tang–Lu theory is, on the whole, more accurate than the Barker–Henderson second-order theory for both quantities. One might argue that these deviations are

Fig. 10. As in Fig. 7 for γ = 36 and T ∗ = 0.5, 0.7, 1.0, 1.5, 2.0, and 3.0.

72

A. D´ıez et al. / Fluid Phase Equilibria 253 (2007) 67–73

considerably high in both theories, especially for U E . However, it is to be noted that, for the cases considered here, the excess energy is usually a small quantity, frequently close to zero, so that a small absolute deviation in the calculated values with respect to simulations ones result in a large relative deviation. A similar situation arises for Z at low temperatures. The advantage of the TL theory over the BH one is more clearly seen when we analyze the radial distribution function g(r). Although in its original formulation [17] the BH perturbation theory did not provide an explicit expression for the radial distribution function, it may be obtained from the relationship [16]:

List of symbols c(r) direct correlation function F free energy g(r) radial distribution function h(r) total correlation function H(r) Heaviside step function u(r) intermolecular potential reduced temperature T∗ UE excess internal energy x = r/σ reduced distance Z compressibility factor

∂F 1 = Nρg(r) ∂u(r) 2

Greek letters  potential depth γ exponent in the Sutherland potential η packing fraction ρ∗ numerical density ρ in reduced units σ diameter of the particles

(33)

From the second-order BH perturbation theory [17] in the mc approximation this gives [16]  

∂ρ gmc (r) = g0 (r) 1 − u(r) (34) ∂p 0 and in the lc approximation   ∂[ρg0 (r)] glc (r) = g0 (r) − u(r) ∂p 0

(35)

where subscript 0 refers to the hard-sphere reference system. Expressions (34) and (35) are temperature-independent. Therefore, contrarily to the Tang–Lu theory, and to the simulation results, the Barker–Henderson theory does not provide a temperature-dependent radial distribution function. In summary, the Tang–Lu procedure provides an analytical expression for the radial distribution function of Sutherland fluids that is qualitatively right and leads to accurate results for the excess energy and the compressibility factor of these fluids for wide ranges of densities and temperatures. The theory compares favourably with the Barker–Henderson perturbation theory for these fluids.

Subscripts 0 hard-sphere reference system contribution 1 first-order contribution 2 second-order contribution Acknowledgement We are grateful to the Spanish Ministerio de Ciencia y Tecnolog´ıa (MCYT) for the financial support under grant no. BFM2003-001903. Appendix A Table A.1 contains the simulation data for the radial distribution function of Sutherland fluids obtained in the present work and used in Figs. 1 and 2. The estimated statistical uncertainty is approximately 0.3% for T ∗ = 0.5 and 0.6 and 0.2% for T ∗ = 1.

Table A.1 Monte Carlo simulation data for the radial distribution function of Sutherland fluids at reduced number density ρ∗ = 0.9 x

1.00 1.01 1.02 1.03 1.04 1.05 1.07 1.10 1.15 1.20 1.25 1.30 1.35 1.40 1.45 1.50

g(x) γ = 6, T ∗ = 0.6

γ = 12, T ∗ = 0.5

γ = 18, T ∗ = 0.5

γ = 36, T ∗ = 0.5

γ = 6, T ∗ = 1

γ = 12, T ∗ = 1

γ = 18, T ∗ = 1

γ = 36, T ∗ = 1

6.840 5.980 5.264 4.661 4.150 3.721 3.037 2.326 1.613 1.215 0.970 0.820 0.727 0.669 0.639 0.630

10.113 7.875 6.265 5.101 4.246 3.603 2.722 1.949 1.327 1.032 0.872 0.779 0.724 0.688 0.675 0.678

12.379 8.750 6.448 4.963 3.972 3.288 2.427 1.752 1.264 1.030 0.896 0.808 0.755 0.720 0.697 0.694

16.707 9.408 5.939 4.201 3.266 2.710 2.088 1.654 1.307 1.100 0.965 0.860 0.790 0.746 0.719 0.710

6.084 5.443 4.890 4.412 3.998 3.635 3.047 2.397 1.707 1.290 1.029 0.860 0.751 0.683 0.645 0.630

7.280 6.176 5.301 4.606 4.051 3.597 2.908 2.227 1.578 1.217 0.993 0.849 0.755 0.694 0.660 0.651

8.191 6.602 5.443 4.596 3.959 3.467 2.772 2.124 1.533 1.206 1.002 0.864 0.771 0.713 0.678 0.664

9.661 6.891 5.259 4.265 3.621 3.176 2.600 2.085 1.583 1.262 1.047 0.896 0.791 0.720 0.679 0.663

A. D´ıez et al. / Fluid Phase Equilibria 253 (2007) 67–73

73

Table A.1 (Continued ) x

g(x) γ = 6, T ∗ = 0.6

1.55 1.60 1.65 1.70 1.75 1.80 1.85 1.90 1.95 2.00 2.05 2.10 2.15 2.20 2.25 2.30 2.35 2.40 2.45 2.50 2.55 2.60 2.65 2.70 2.75 2.80 2.85 2.90 2.95 3.00

0.646 0.682 0.742 0.831 0.943 1.053 1.129 1.180 1.227 1.298 1.287 1.189 1.088 1.001 0.936 0.893 0.867 0.857 0.862 0.881 0.910 0.948 0.988 1.028 1.061 1.082 1.089 1.086 1.076 1.059

γ = 12, T ∗ = 0.5 0.696 0.737 0.792 0.874 0.987 1.077 1.130 1.172 1.226 1.311 1.243 1.114 1.017 0.953 0.913 0.891 0.882 0.883 0.896 0.918 0.946 0.977 1.011 1.039 1.061 1.071 1.072 1.066 1.056 1.041

γ = 18, T ∗ = 0.5 0.707 0.736 0.789 0.871 0.980 1.060 1.118 1.168 1.230 1.313 1.217 1.096 1.014 0.963 0.929 0.906 0.897 0.896 0.902 0.919 0.943 0.972 1.003 1.031 1.051 1.065 1.069 1.064 1.054 1.038

γ = 36, T ∗ = 0.5 0.715 0.746 0.795 0.863 0.958 1.032 1.091 1.150 1.221 1.308 1.185 1.091 1.033 0.986 0.951 0.925 0.911 0.908 0.912 0.921 0.941 0.966 0.991 1.018 1.041 1.056 1.059 1.058 1.053 1.041

References [1] [2] [3] [4] [5] [6] [7] [8] [9]

D.M. Heyes, L.V. Woodcock, Mol. Phys. 59 (1986) 1369–1388. P.J. Camp, G.N. Patey, J. Chem. Phys. 114 (2001) 399–408. P.J. Camp, Phys. Rev. E 67 (2003) 0115031–0115038. A. D´ıez, J. Largo, J.R. Solana, J. Chem. Phys. 125 (2006) 0745091– 07450912. V.I. Kurochkin, Sov. Phys. Ledebedev Inst. Rep. (U.S.A.) 8 (1990) 1–2. J. Largo, J.R. Solana, Int. J. Thermophys. 21 (2000) 899–908. S. Jiuxun, Can. J. Phys. 83 (2005) 55–66. M. Edalat, S.S. Lan, F. Pang, G.A. Mansoori, Int. J. Thermophys. 1 (1980) 177–184. J. Cho, I.C. Sanchez, Macromolecules 31 (1998) 6650–6651.

γ = 6, T ∗ = 1

γ = 12, T ∗ = 1

γ = 18, T ∗ = 1

γ = 36, T ∗ = 1

0.638 0.668 0.722 0.803 0.910 1.021 1.107 1.172 1.225 1.291 1.293 1.212 1.116 1.032 0.961 0.911 0.876 0.860 0.858 0.871 0.896 0.929 0.970 1.011 1.047 1.075 1.088 1.090 1.084 1.070

0.661 0.694 0.745 0.824 0.926 1.027 1.107 1.168 1.220 1.291 1.277 1.182 1.088 1.012 0.950 0.907 0.881 0.869 0.873 0.887 0.910 0.942 0.979 1.018 1.050 1.072 1.083 1.081 1.073 1.059

0.670 0.699 0.748 0.821 0.919 1.017 1.094 1.161 1.223 1.297 1.268 1.171 1.082 1.012 0.958 0.916 0.891 0.879 0.880 0.891 0.913 0.944 0.977 1.011 1.042 1.063 1.074 1.076 1.072 1.060

0.669 0.695 0.741 0.813 0.905 1.000 1.083 1.151 1.219 1.299 1.261 1.177 1.098 1.026 0.970 0.928 0.897 0.881 0.879 0.888 0.908 0.936 0.970 1.003 1.036 1.059 1.072 1.078 1.074 1.063

[10] J.A.P. Coutinho, P.M. Vlamos, G.M. Kontogeorgis, Ind. Eng. Chem. Res. 39 (2000) 3076–3082. [11] D. Ben-Amotz, G. Stell, J. Chem. Phys. 119 (2003) 10777–10788. [12] N. Elvassore, J.M. Prausnitz, Fluid Phase Equilbr. 194–197 (2002) 567–577. [13] M. Murakami, J. Chem. Phys. 120 (2004) 6751–6755. [14] Q. Zhang, L.A. Archer, J. Chem. Phys. 121 (2004) 10814–10824. [15] Y. Tang, B.C.-Y. Lu, J. Chem. Phys. 99 (1993) 9828–9835. [16] Y. Tang, B.C.-Y. Lu, Mol. Phys. 90 (1997) 215–224. [17] J.A. Barker, D. Henderson, J. Chem. Phys. 47 (1967) 2856–2861. [18] J.D. Weeks, D. Chandler, H.C. Andersen, J. Chem. Phys. 54 (1971) 5237–5247. [19] Y. Tang, J. Chem. Phys. 118 (2003) 4140–4148. [20] Y. Tang, J. Wu, Phys. Rev. E 70 (2004) 0112011–0112018.