Retension of the bridge diagrams in the pair-distribution function's integral equation

Retension of the bridge diagrams in the pair-distribution function's integral equation

Physica A 181 (1992) 297-308 North-Holland I [~1 Retention of the bridge diagrams in the pair-distribution function's integral equation M.-E. B o u...

547KB Sizes 0 Downloads 10 Views

Physica A 181 (1992) 297-308 North-Holland

I

[~1

Retention of the bridge diagrams in the pair-distribution function's integral equation M.-E. B o u d h - H i r The University of Illinois at Chicago, CEB, 810 South Clinton Street, Chicago, 1L 60607, USA

Received 9 January 1991 Revised manuscript received 27 August 1991

A convenient theoretical method for calculating the structure properties of a homogeneous fluid, retaining the bridge diagrams, is developed. Using the diagrammatic expansion theory, the exact pair-distribution function's integral equation is reduced to a hypernetted chain-like equation in which the interaction between particles is replaced by an effective pair potential depending on the bridge function. The bridge diagrams are classified according to their number of field points. The contributions of these different diagrams are analyzed and an approximation for the bridge function is proposed.

I. Introduction

Considerable theoretical efforts have been expended in order to gain knowledge on the structure of homogeneous fluids. Several approximated integral equations such as the mean spherical approximation (MSA) [1], the hypernetted chain (HNC) [2-7] and Percus-Yevick (PY) [8] approximations as well as the HNC and PY derivatives have been proposed for the pairdistribution function (PDF), g(1,2). These integral equations are now well established and their weaknesses are known. It has been shown that the HNC approximation is generally more accurate than the PY approximation for systems in which the particles interact via long-ranged potentials [9], while PY approximation leads to better results with short-ranged repulsive interactions [10]. In fact, a fine comparison of HNC and PY approximation results for hard sphere systems shows the following [11]: (i) At low density (i.e., 0 0 -3 =- 0.4) the PY approximation gives results in excellent agreement with those obtained from computer simulations. (ii) At fairly high density (i.e., p0-3= 0.9) the accuracy of these approximations' results depends on the distance, r~2, between the two particles 1 and 2. Nevertheless, PY approximation always leads to 0378-4371/92/$05.00 © 1992- Elsevier Science Publishers B.V. All rights reserved

298

M.-E. Boudh-Hir / An approximation Jbr the hrid~e [kmction

results in better agreement with simulations for small values of r~2 (i.e., #'12 slightly larger than (r, where ~r denotes the hard spherc diametcr of the particles). Hence, PY is the appropriate approximation for these systems since their thermodynamic properties are related to the PDF contact value. Concerning the MSA approximation, perhaps its only advantage is that it may be solved analytically. Considering combinations of HNC find PY approximations, attempts have been made in ordcr to find a self-consistent approximation giving results in better agreement with those obtained by numerical simulations [12]. The coefficients that this method generates have to be calculated in such a way that guarantees the thermodynamic consistency. These coefticients have been assumed to be dependent only on the fluid density, p, and chosen such that the same virial coefficients fire obtained using the virial or the compressibility equations of state [12]. The first coefficients have been calculated and the improvement is considerablc. To extend successfully the HNC approximation, it has been suggested that, instead of neglecting it, the bridge function should be replaced by its value for an assumed known system [13]. It should however be noted that although the bridge function is a short-ranged one for homogeneous systems (indeed, for large values of the distance between particles, it decreases at least as the square of the pair potential) it is not a universal function. In fact, it may be: (i) A function depending only on the particle particle distance. That is the case of systems in which the interactions between particles depend only on their mutual separations. (it) An anisotropic function, which depends on the distance between the particles and on their orientations. This for instance is the case of polar fluids. The most interesting property of the bridge function is the fact that. for a dipolar fluid in front of a hard wall, it becomes a long-ranged function (i.e.. it decreases at infinity as the particle-particle potential) [14]. This point provides the proof that the bridge diagrams do not have a universal contribution. Hence, the suitable reference system to study a given fluid will depend on the nature of the interactions between its particles. Since the only system for which the bridge function is known with sufficient accuracy is the hard sphere fluid, the knowledge of this function, for this single case is not very useful because of its non-universality. For instance, when the above mentioned approximation is introduced in the treatment of polar fluids, the bridge function which is an anisotropic function is replaced by an isotropic one. Therefore, the only anisotropic components that the bridge function will generate are obtained through the convolution (i.c., given by the O r n s t e i n - Z e r n i k e equation). Hence, if this contribution is expanded in spherical harmonics, it will not cover all the possible symmetries.

M.-E. Boudh-Hir / An approximationfor the bridgefunction

299

2. Theoretical developments We consider a system of N identical particles interacting via an additive pairwise potential, v(i, j). The two-particle density, p(1, 2), is defined by p(1,2) = -~ ' Z

N(NNwl)zN

fexp(-t3,,jv(i,

j ) ) d { N - 2 } = p2g(1 2)

N>~0

in which the PDF, g(1,2), is given by the graphical expansion: g(1,2) = exp[-/3v(1,2)] exp(sum of all the distinct connected diagrams free of: articulation points, consisting of two non-directly connected root points 1 and 2, p field points and f bonds, such that the two root points do not form an articulation pair) = exp[-/3v(1,2)] e x p ( & + ~ +

~:J~+ ~7~+ ~ + . . . ) .

(2)

Here the white and the black circles denote the root and the p field points, respectively, and the lines represent the f bonds. The Mayer functions, f, are associated to the pair-potential via the following equation:

f(i, j) = exp[-/3v(i, j)] - 1.

(3)

Two non-directly connected points are defined as a pair of points between which there is no direct path without intermediate points. An articulation point is defined such that upon its removal from a connected diagram, the diagram splits into two or more components. Of the resulting components, at least one will contain no root point. An articulation pair is defined such that upon its removal from a connected diagram, the diagram splits into two or more components free of root point. The diagrams in eq. (2) can be classified into two groups according to their structure: (i) Bridge or elementary diagrams, they are free of nodal points. We denote their contribution by B(1, 2). (ii) Nodal or series diagrams, they have at least one nodal point. Their contribution, N(1, 2), is therefore given by a convolution. It is easy to see that it satisfies

N(1, 2) = h ( l , 2 ) -

c(1,2) = p f c(1, 3) h(3,2) d3,

(4)

M . - E . B o u d h - t t i r / A n approximation ~ r the hridL, e fimcti(m

300

where c(1,2) and h(1,2) are the direct and the total correlation functions, respectively. A nodal point is defined such that upon its removal from a connected diagram, the diagram splits into two componcnts each of them having one root point. Using the two functions B(1,2) and N(I,2) the PDF, g(1,2), can be expressed by g(1,2) = exp[-/3v(1,2) + B(1,2) + N(1,2)].

(5)

Eqs. (4) and (5) together with

g(1,2)-h(1,2)+ 1,

(6)

lead to the following expression for the direct correlation function, c( 1,2):

c(1,2)-f(1,2)exp[N(l,2)]+f(1,2){exp[B(l,2)] + {exp[B(1,2)]-1} exp[N(l,2)] + ,~, pl

1}exp[N(l,2)]

N"(

n!1,2)

I

(7)

At this point it is useful to recall and compare the HNC and PY approximations. This may allow us to gain some insight on how to implement the treatment of the contribution of the bridge diagrams.

2.1. HNC and PY approximations The HNC approximation assumes that the bridge function, B(I,2), is identically zero. Eq. (7) becomes c~N~(l,2)

f(l,2)exp[XHN"(l,2)]+ ~ [NHY"(I,2)] '' J~ I

(~)

11 !

It follows that the HNC approximation becomes more accurate for particle systems of which the bridge function is a short-ranged function having a small magnitude. The PY approximation is given by c e ~ ( 1 , 2 ) - f ( 1 , 2 ) explNe~(l, 2) I +./'(1,2) {explB"~(1, 2)1

1} exp[NF'V(I,2)].

(9)

One may verify that the functions NeV(1,2) and Bv~(I, 2) satisfy the condition

M.-E. Boudh-Hir / An approximation for the bridge function

exp[BPV(1, 2) + NPv(1, 2)] - [ 1 + NeV(1, 2)] = 0,

301

(10)

which can be written in the form BPV(1,2) = - ~ n>l

[-NPV(I' 2)l"

(11)

F/

Concerning these approximations it should be noted that: First, the nodal diagrams' contributions, NHNC(1,2) and NPv(1,2), in eqs. (8) and (9), respectively, are different because they satisfy eq. (4). Second, the first term in eq. (7), f(1, 2) exp[N(1,2)], which constitutes the main contribution to the direct correlation function is retained by both the approximations. A little difference exists, however. It is due to the fact that the nodal contributions, as mentioned above, are not identical. This difference is, nevertheless, of higher order in pair potential. Indeed, the lowest order is the same and is equal to -/3v(1, 2). Third, the fact that the third term in eq. (7) is neglected in both of these approximations, should have the same consequences on HNC and PY results. Fourth, in homogeneous systems the contribution of the bridge diagrams is always a short-ranged function. This is the consequence of their high connectivity. In contrast, the contribution of the nodal diagrams has the same range as the pair potential except in the case of Coulombic systems because of the electroneutrality condition. Considering the first term in the right hand side of eq. (7) as the main part of c(1, 2) and the last three terms as a correction, it follows that: (i) In the case of long-range potentials, the contribution of the last term is the most important. (ii) For short-range potentials the two last terms, which become of the same order, are the most important. Consequently, the accuracy of the HNC and PY approximations with respect to the nature of the interactions between the particles of the system can really not be attributed to the difference in numbers of diagrams retained or to a fortuitous cancellation of errors. Indeed, regardless of the kind of particleparticle potential, the term neglected by HNC, exp[-/3v(1, 2)] e x p [ B ( 1 , 2 ) 1] exp[N(1, 2)], is really negligible inside the hard core because of the Boltzmann factor and goes rapidly to zero outside because of the bridge function. This term is essentially non-zero for distances between the particles slightly larger than o-. Its contribution to the total correlation function through Ornstein-Zernike equation which is, of course, neglected by HNC is small at low density and large distances. This term affects the PDF essentially at small values of the distance, r12, between the two particles 1 and 2. Therefore, it is

302

M.-E. Boudh-ttir / An approximation .for the bridge .flmction

not expected to have good results for the PDF contact value using HNC. PY approximation assumes that B(1,2) and N(1,2) satisfy eq. ( 11 ). This is clearly not realistic for long-range potentials. However, for short-ranged repulsive interactions, both of these functions decrease rapidly to zero with respect to the distance, rte, and the approximation becomes more reasonable. We now turn back to eq. (2) and consider a bridge diagram of this expansion. We note that the bridge structure of this diagram is not changed by replacing an f bond by a bridge subdiagram. A bridge subdiagram is such that, when isolated from the main diagram, it leads to a bridge diagram. We therefore define the function f * as ,l"~'(i, j ) = exp[B(i, j ) ] -

(12)

1,

and expand the PDF in terms of functions f and f*. The PDF's diagrams obtained are free of bridge diagrams and subdiagrams. We have

g ( l , 2 ) = e x p [ - / 3 v ( l , 2) + B(I, 2) 1 exp(sum of all the distinct connected diagrams free of: articulation points and bridge articulation pairs, consisting of: two non-directly connected root points 1 and 2, p field points, and f a n d / o r f * bonds such that the two root points do not form an articulation pair) cxpl-/~v(l, 2) + R(l, 2)1

(13) Here the solid and the dashed lines represent the f and the .f.... bonds, respectively. A bridge articulation pair is defined such that upon its removal from a connected diagram, the diagram splits into two or more components. Of the resulting components, at least one will be a bridge diagram. The expansion (13) is precisely the HNC approximation of the PDF of a system of particles interacting via the effective potential, v*(i, j ) , defined by v * ( i , j ) = v(i, ./) - k T B ( i ,

j) .

(14)

The P D F , g(l, 2), can be written rigorously in the form g ( l , 2) = e x p [ - / 3 v * ( l , 2) + N*(I, 2)],

(15)

M.-E. Boudh-Hir / An approximation for the bridge function

303

where N*(1,2) is the contribution of the nodal diagrams for a system of particles interacting via the effective pair potential, v*(i, j). Hence, it is clear that the exact PDF's integral equation can be reduced to an HNC-like equation provided that the particle-particle interaction is replaced by an effective pair potential which itself depends on the bridge function according to eq. (14). The eternal problem is how to calculate this function which plays here a central role. Unfortunately, this calculation is not possible at the present time and, consequently, approximating methods become unavoidable.

2.2. Approximation of the bridge function As a consequence of the relatively simple structure of the nodal diagrams, the function N(1,2) has a compact and simple form. Indeed, the fact that eq. (4) cannot be solved is not related to the equation itself but is a consequence of the closure problem (i.e., if the direct and the total correlation functions were known, the convolution would be easy to perform). Due to the high connectivity of the bridge diagrams, such a simple form for the function B(1,2) cannot be obtained. This function can however be expressed in terms of residual diagrams of n-particle distribution functions with n >~4, but this does not simplify the calculation of B(1, 2) and therefore it does not make any difference regarding the determination of the PDF. It follows that the only remaining alternative is to find a criterion according to which the bridge diagrams can be classified and their contribution analyzed, and then, on the basis of this analysis, find an approximation for the bridge function. For simplicity, we will first proceed to the topological reduction which eliminates the Mayer functions, f, in favour of the total correlation functions, h, according to the fourth and sixth lemma of Morita and Hiroike [15]. The diagrams obtained are therefore free of field points having less than three h bonds. This topological reduction has the following advantages: (i) It does not modify the structure of the diagrams, that is, it replaces a series of bridge diagrams by a single bridge diagram (and a series of nodal diagrams by a single nodal diagram). As a consequence, the number of diagrams of each of these two classes is reduced which should simplify their analysis. (ii) It replaces the Mayer function by the total correlation function whose variations are, in general, more moderated. This should accelerate the convergence of the developments. It is true that the total correlation function introduced is an unknown function but, for our purpose, the knowledge of the general behaviour of this function is sufficient (i.e., the details are not necessary). Consider the PDF's bridge diagrams characterized by the same number, n,

304

M.-E. Boudh-Hir / An approximation for the bridge ]}mction

of field points. It can be noted that the simplest bridge diagrams of this series can be obtained as follows: First, the n field points are connected to each other in such a way to form a chain. Second, the two field points in the ends of the chain are connected to both of the root points. Third, one of the n - 2 remaining field points is connected to one of the root points, the others are connected (i) to each other or (ii) to any of the root points such that each of the root point has a minimum of two h bonds and each of the field points has, at least, three h bonds. The contribution of all the bridge diagrams having n field points is simply obtained from these basic diagrams by: (i) connecting by g bonds any non-connected pair, or (ii) replacing by g bonds all the h bonds whose removal left the bridge structure of the diagram unchanged. Let us, for simplicity, consider a simple fluid characterized by the parameter (r (which stands for the hard sphere diameter for hard sphere fluids or the collision parameter, otherwise) and note that: (i) In a given bridge diagram, each field point has, at least, three h bonds. The function h(r,~) is - 1 for distances smaller than o-. It presents a main peak. whose magnitude depends on the density of the fluid, at r~, slightly larger than (r and then vanishes after few oscillations. (ii) The functions g and h have opposite effects. Indeed, the function 9(r,i) is really not negligible only for distances rii larger than {r while h(r,) may have an important contribution only for r,~ smaller or slightly larger than (r. It follows that the main contribution to the bridge function comes from the configurations in which: - The field points having g and h bonds move inside small spheres of radius ~O'. - The field points free of g bonds move inside small spheres of radius l = A(r. The parameter l stands for the average range of the total correlation function h(r~i ). This parameter is of the order of {r and is, because of the above mentioned properties of h(r,i), a decreasing function of the density of the fluid. It follows that the contribution, b ( l , 2; n), of a bridge characterized by n field points, k of which have g and h bonds, can be estimated by b( l, 2; n) ~ r/'a 3°' k)~3aA(l,2:n).

(16)

In this equation, r/ denotes the packing fraction of the fluid defined by

vl = ], qvl(r' ,

(17)

and A ( 1 , 2 ; n) is a function of magnitude of the order of unity, depending on

M.-E. Boudh-Hir / An approximation for the bridge function

305

the distance rij and the parameters n and p. The sign of b(1, 2; n) is essentially governed by the number of h bonds that the diagram has. Indeed, each h(ri~) contributes by a factor of essentially - 1 to this function. The function magnitude converges rapidly to zero. In fact, the number k is larger than or equal to 2 (except in the case n = 2 where it is zero). Consequently, only the first few bridge diagrams will contribute effectively to B(1,2). Here, this function is approximated by the contribution of the simplest bridge diagram. We have B(1,2)~

½p2 f

h(1,3) h(1,4) h(3,4) h(2, 3) h(2,4) d 3 d 4 .

(18)

The most important contributions to this integral come from the domains: D 1 (in which the particles 3 and 4 are in the neighbourhood of one of the particles 1 and 2) and D 2 (where each of the particles 1 and 2 has in its neighbourhood one of the two particles 3 and 4). The contributions, i~(1,2) and i2(1, 2), of these domains are estimated by i,(1, 2) ~ - o~(p, T) rt2h2(1, 2),

(19a)

i2(1, 2) ~- co(p, T) "q2h3(1, 2),

(19b)

where the coefficients a(p, T) and co(p, T) depend on the density of the fluid and its temperature (except in the case of hard-sphere fluids where they become functions of the density only). The magnitudes of these coefficients are expected to be of the order of unity, increasing functions of temperature and decreasing functions of density. The coefficients ~(p, T) and co(p, T) can be related through the approximate relations

T)= co(p,

T),

(20a)

T),

(20b)

in which the function I(p, T) depends essentially on the connectivity of the diagram and, therefore, it is the same for a(p, T) and co(p, T). The coefficient ~"comes from the total correlation functions between particles in two different domains. (These total correlation functions are approximated by h(1, 2). Those between particles in the same domains D 1 or D 2 are replaced by their average value, i.e., -1.) As a consequence of the oscillating form of the total correlation function, ~" should be a decreasing function of the density of the fluid. Using eqs. (20a) and (20b), a(p, T) and co(p, T) can be expressed by co(p, T ) ~ ~'2a(p, T ) .

(21)

306

M.-E. Boudh-Hir / An approximation lbr the bridge fhnction

Inserting eqs. (19a) and (19b) in (18), the bridge function takes the form B ( 1 , 2 ) = ~/21- c~(p, T) h e ( l , 2) + oo(p, T ) h3(l,

2)1.

(22)

which t o g e t h e r with eq. (14) leads to the effective pair potential u*{i, j ) = v{i, j) + kTr/~[-c~{ p, T ) h ' - ( 1 , 2 ) + w { p , T) h ~(1,2)1

= v*(i, ./: c~, co).

(23)

T h e P D F of the system is now a function of the p a r a m e t e r s (t and {o aiven by g ( 1 , 2 : c~, ~o) = e x p l - ~ v * ( 1 , 2 :

{r, {o)+ N * ( I , 2: ~, ~o}].

{24)

This function can bc o b t a i n e d by solving the familiar equations (4), (8) and (24) defining the H N C a p p r o x i m a t i o n t o g e t h e r with eq, (23) giving the pair potential. T o have an idea about the bridge function's a p p r o x i m a t i o n , given by eq. (22), the coefficients {~(p, 7") and £ arc roughly estimated and the calculations of B( l, 2) are p e r f o r m e d for hard-sphere systems in both cases of low and high density. In both cases this function is, as expected, negative (sec tables I and I1). At low density, B ( 1 , 2 ) is reduced to - r / - ' E : ( l , 2) with E_,(1,2) being the c o n t r i b u t i o n of first bridge diagram. This result is available in the literature [161. T h e c o m p a r i s o n of our results with the exact ones [16] shows a qualitative a g r e e m e n t (table i). At high density, such a result is not available, at least to o u r knowledge. T h e r e f o r e , we have to confine ourselves by examining the general b e h a v i o u r of the bridge function which turns out to be satisfactory (table I1). Table 1 Bridge function {}1 a hard-sphere lluid at p(r: = (I.4. ~ = 3. w = 2.3, ~"=: 0.87. ,.:,r

h(r)[lll

1.0 1.1 1.2 1.3 1.4

{).77 0.56 {).36 0.23 0.1 I (1.11(} 0.04 0.07 I).09 0.07 0.03

1.5

1.6 1.7 1.8 1.9 2.0

I~{,),',F'

tl.727 0.537 I).281 0. 131 (t.(133 0.000 0.1}05 0.014 0.026 I).(114 {}.{}03

E:{r) [l~q 11.088 0.494 11.331 I).225 O. 125 (I.tl62 0.025 11.013 0.1)00 0.000 {}.00{1

M.-E. Boudh-Hir / An approximation for the bridge function

307

Table II Bridge function of a hard-sphere fluid at po-3=0.9. a =0.75, ~o =0.19, ~r = 0.51.

r/o

h(r) [111

B(r)

1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4

3.26 1.77 0.04 0.06 -0.26 -0.40 -0.39 -0.26 -0.11 0.06 0.26 0.29 0.14 0.03 -0.04

- 1.010 -0.611 0.001 -0.012 -0.025 -0.062 -0.058 -0.025 -0.004 -0.001 -0.025 -0.027 - 0.007 -0.000 -0.001

A more rigorous way to determine the coefficients u(p, T) and w(p, T) and, therefore, a systematic calculation of B(1,2) and g(1, 2) should consist in imposing to the PDF a supplementary constraint that requires the thermodynamic consistency. For instance, these coefficients are such that the virial and the compressibility equations lead to the same pressure. This constraint can be stated by the following equations:

2 2 i r~g(r; , [3P = p - ~'rrp o~, w) Orv(r ) d r ,

(25a)

0

(~OoP)r,v -- 1 - 4~rp f r2c(r; ~, ~o) d r ,

(25b)

1)

where the symbol 0,~ denotes the derivative with respect to x. Consequently to the behaviour of the total correlation function, it can be noted, regarding the contribution of the bridge function to the effective potential, v*(1, 2), that: (i) Inside the core, the contribution of the bridge diagrams is a bonded function. Hence, the effective potential, v*(i, j), has the same behaviour as the pair potential, v(i, j), in this region. (ii) For isotropic interactions, v(i, j) = v(ri~ ), this contribution is repulsive for distances between particles slightly larger than o-; and alternatively attractive and repulsive for the larger values of the interparticle distances. Due to the

31)8

M.-E. Boudh-Hir / An approximation/br the bridge/unction

fact t h a t t h e m a g n i t u d e of the t o t a l c o r r e l a t i o n f u n c t i o n goes to z e r o for r~, l a r g e r t h a n (r, the effect of the p a i r - p o t e n t i a l , v(i, j ) , d o m i n a t e d a n d , t h e r e f o r e , this c o r r e c t i o n will h a v e no d i r e c t c o n s e q u e n c e in this region ( i . e . . the m o d i f i c a t i o n o f the P D F in this r e g i o n arises t h r o u g h the c o n v o l u t i o n ) . (iii) F o r a n i s o t r o p i c p a i r p o t e n t i a l s , v(i, ]) = v(r~i, ,Q~, .Q~), with ~(2, d e n o t i n g t h e o r i e n t a t i o n s of the m o l e c u l e s i, this c o r r e c t i o n has an i s o t r o p i c a n d a n i s o t r o p i c c o m p o n e n t s . It will t h e r e f o r e have a n o n - z e r o c o n t r i b u t i o n to lhe o r i e n t a t i o n a l p a r t of the P D F . A s in the p r e v i o u s case, this c o r r e c t i v e t e r m is d o m i n a t e d by the p a i r p o t e n t i a l in the region r~i < (r: and its i s o t r o p i c p a r t is a r e p u l s i v e p o t e n t i a l for r~i slightly l a r g e r than ~r.

3. Concluding remarks T h e p u r p o s e of this p a p e r is to retain the c o n t r i b u t i o n of the b r i d g e d i a g r a m s to t h e P D F . By m e a n s of the d i a g r a m m a t i c e x p a n s i o n t h e o r y , the e x a c t P D F ' s i n t e g r a l e q u a t i o n is r e d u c e d to an H N C - l i k c e q u a t i o n in which the i n t c r a c t i o n b e t w e e n p a r t i c l e s is r e p l a c e d by an effective p a i r p o t e n t i a l d e p e n d i n g on the b r i d g e function. T h e b r i d g e d i a g r a m s are classified a c c o r d i n g to t h e i r n u m b e r of field p o i n t s . T h e c o n t r i b u t i o n s of the d i f f e r e n t d i a g r a m s are a n a l y z e d and an a p p r o x i m a t i o n for the b r i d g e function is p r o p o s e d .

References [ I I J. I,. l,cbowitz and ,I.K. I'crcus, t'hys. Rcv. 144 (1966) 251. 121 J.M.J. Wm Lccuwcm J. Grocnvcld and J. l)e Boer. Physica 25 (I959) 792. 131 E. Mecrom J. Math. Phys. I 11961)) 192. [41 M.S. Green, J. Chem. Phys. 33 (19611) 1403. [5] T. Morita and K. Hiroikc, Progr. Thcor. Phys. 23 (1960} 111t)3. [6] L. gerlct. Nuovo Cimcnto 18 (1960) 77. [7] G.S. Rushbrookc, Physica 28 (1960) 259. [8] J.K. Porous and G..I. Ycvick, Phys. Roy. 110 (lg5S) 1. [9] J.-P. Hansen and I.R. McDonald, Phys. Rcv. A 11 (1975) 2111. [1t)] A.A. Broyles, S.U. Chung and tl.L. Sahlin, J. Chem. Phys. 37 (1962} 2462. [I I] R.O. Watts and l.J. McGee, Liquid State Chemical Physics (Wiley, New York, 197fi) p. 137. [121 ,1.S. Rowlinson, Mol. Phys. 9 (1965) 217. [13] F. Lado, Phys. Rcv. A 8 (1973) 2548. [14] M.-E. Boudh-Hir, Physica A 158 (1989) 619. [15] T. Morita and K. Hiroikc, Progr. Thcor. Phys. 25 (1961) 537. [16] F.H. Rcc, R.N. Keclcr and S.L. McCarty, J. Chem. Phys. 44 (1966) 3407.