- Email: [email protected]
First-order correction to the three-body correlation function
Volume 183, number 5
6 September 199 I
CHEMICAL PHYSICS LETTERS
First-order correction to the three-body correlation function Ivo Nezbeda ’ and Gerhard Kahl Institut ftirTheoretische Physik, Technische Universiidt Wien, Wiedner Hauptstrasse 8- 10, A- 1040 Vienna, Austria Received 25 February 1991;
in final form 12 June 1991
The first-ordercorrection to the Kirkwood superposition approximation of the triplet correlation function, based on the density expansion due to Abe, has been calculated for the high-density Lennard-Jones fluid. While qualitative agreement with computer simulation data has been obtained, substantial disagreement with older literature data and conclusions based thereon have been found.
Although there has been much interest in the triplet correlation function g, both from the practical and theoretical points ov view, it is the great complexity of this field which has so far prevented an exhaustive treatment (for a review, see e.g. ref. [ 1 ] ). Nowadays, when the determination of the pair correlation function g of simple fluids has reached a high level of sophistication and computers have a reasonable performance rate to evaluate complex expressions, the time seems to be ripe for re-investigating the traditional approaches before adopting complex new methods via the density functional theory and direct correlation functions (see refs. [ 2,3 ] and references therein). The function g may be most conveniently written as
g3(r12,r13,b3)=g3(r, s, 0 =g(rk(sMWVc
4
1)
=ih(r, s,WV6 s, t) ,
(1)
where the subscript SA denotes the Kirkwood superposition approximation. Although the SA (i.e. R= 1)was for years the most commonly used approximation for g3, it was shown that it is substantially in error at high densities and, similarly, various
’ Permanent address: Institute of Chemical Process Foundations, Czechoslovak Academy of Sciences, 16502 Prague, Czechoslovakia. 0009-2614/91/S
03.50 0 1991 Elsevier
corrections based on semi-empirical expressions of R [1,4] also failed. Rigorously, may R be expressed as [ 51 R(r,s,f)=exp[pr(r,s,t)l =exp
F p”m(r,s,1)
( i=l
,
(2)
>
where p is the number-density and T(r, s, t) is an infinite sum of “simple 123 irreducible clusters” with f-bonds cf( r) = exp [ - /.?u(r) ] - 1 is the Mayer function and u(r) is the pair potential). I1 is possible to use also h-bonds, h =g- 1, instead of the f-bonds in the diagrams of (2)) which leads to a substantial reduction in the number of diagrams contributing to T(r, s, t), but then eq. (2) is no longer a density expansion. It is easy to see that for n& 3, the terms exceed practical applicability. Eq. (2 ) may now be seen from at least two points of view: (i) One may use eq. ( 1) as a closure in the BGY-hierarchy [ 6 ] and generate g and g3 in a selfconsistent way; this route was followed by Haymet et al. [ 71 at the time when no accurate methods for determining g( r) were available. (ii) With methods now available enabling one to determine g of a simple fluid with high accuracy [ 8 1, we can use eqs. (2) and ( 1) as defining equations for g,. Although the first route turned out to be not so bad, both approaches were practically abandoned at the beginning of the 1980’s: despite enormous effort, the results obtained did not reach expectations and the
Science Publishers B.V. All rights reserved.
337
Volume183.number 5
6 September1991
CHEMICALPHYSICSLETTERS
existing discrepancies between theory and experiment prevailed. Specifically, in 1982 Sane [ 91 showed that the rigorous first-order correction to SA,
R(r,s,r)IIR,(r,s,I)~exp[PT,(r,S,f)l,
9(r) r
3
2
with 1
is not even in “qualitative agreement” with Monte Carlo data. It is the purpose of this communication to show that Sane’s results are erroneous and, consequently, his conclusions unjustified. We consider the same system and arrangements as in paper [ 91: two state points of the Lennard-Jones fluid, (po3, k,T/e)= (0.75, 1.07) and (0.85, 0.719), each for two different isosceles triangles (r, s, s). To evaluate integral (4), we used two methods: (i) the well-known expansion in Legendre polynomials (used also by Sane) [LO]; and (ii) direct integration. For the latter, we used the Conroy method [ 111 with the new panel of points obtained recently by Nezbeda et al. [ 121 (details of the method will be given in a subsequent paper [ 131). From the dependence of the results on the number of integration points and on other parameters of integration, we claim that the error of the direct integration does not exceed 0.2% and both methods yield the same results within numerical accuracies. To check further the correctness of our numerical method, we also calculated r, using the{bonds in (4) for the same state point and arrangement considered by Rice and Young [ 141 and came to full agreement with their result. The value of ‘TVdepends evidently on h, and since this function used by Sane is not available, we use two different sets of the pair correlation function to investigate also the effect of the input h on TV: one obtained from the extended (self-consistent) MSA (HMSA) [ 151 which was shown to be very accurate [ 16 ] when compared with simulation data, and the other obtained from the HNC theory [6]. In fig. 1, we show both the HMSA and HNC g(r)s for the high-density state point. These functions differ in the heights and positions of the first maximum and minimum and should allow us to see the effect of the input on 7,. In figs. 2 and 3, we compare our results, obtained from the two different inputs using the Conroy method, with the simulation data [ 171 338
\ !
00 0
\\
_’
L 1
2
3
r/o. Fig. 1. The pair-correlation function g(r) of the LennardJones system at pc~~zO.85and kBT/c=0.719, calculated from the HMSA (full line) and HNC theory (broken line).
R(r,4 J) T/U = 2.025 I-
o.5 1.0
3.0
2.0 s/a
Fig. 2. Comparison of the approximation R(r, S, s) uR, (I, s, S) =exp [pr, (r, S, s)] for the LennardJones system atpa3=0.7S and kBT/c= I .07 with the simulation data [ 171 (0) for two different r-values. Eq. (4) with HMSA h (0), HNC h(...), and (--) ref. 191.
and those of Sane [ 91, It is seen that different hs produce only marginal differences in 7,. We have to bear this in mind when discussing the large difference between our and Sane’s data. Since Sane used parameterized molecular-dynamics data (i.e. an accurate input h which, in fact, should not differ from our HMSA input), and anyhow, different hs yield only minor differences in 7,. we conclude that the discrepancies obtained must be exclusively due to numerical inaccuracies.
Volume 183.number 5
CHEMICAL PHYSICS LETTERS
JQ, s,3)
r/u = 1.525
t
0
P/U = 1.025 0.5
i
-\ ‘_’ ’ ’ ’ ’ ! ’ ’ ’ ’ ! ’ ’ ’ ’ ! ’ ’ ’ ’ I I
0.5 1.0
2.0 s/o
3.0
6 September 1991
der correction term (3) provides a considerable improvement over the SA, and for large separations it practically removes the discrepancy between the SA and exact (simulation) g3. These findings put the rigorous expansion (2) in another light, and one may justifiably expect that the inclusion of the second-order term ~~ might push the applicability of eqs. (1) and (2) into quite a high-density range. The results, following this line, will be reported in detail in a subsequent paper [ 131. This work has been supported by the osterreichische Bundesministerium fiir Wissenschaft und Forschung and by the ijsterreichische Fonds zur Fijrderung der Wissenschaftlichen Forschung under Project Number P76 18-TEC.
Fig. 3. The same as for fig. 2 for pa’=O.ES and kJ/c=O.7 19.
References The discussion of the results must be made separately with respect to (i ) the simulation data and (ii) Sane’s data. Concerning (i), one may say that, in general, the first-order term captures the basic features of the full correction term R and thus provides an improvement over the SA. Significant discrepancies remain at separations, roughly, for S/O< 1.8, while at larger separations, the agreement is even nearly quantitative; this is especially true for the lowdensity state point. For the high-density state point, the discrepancy is much larger but nevertheless the trend of the correction is in the right direction. Concerning (ii), we see that at the lower density, our results push the correction even closer to the exact result and at the high density, they considerably remove the drastic qualitative disagreement. Moreover, Sane’s curves exhibit spurious peaks and wiggles. These do not show up in our calculations and may be only an artefact of his numerical method. To be consistent with the paper by Sane [9], we do not show results for other molecular arrangements which exhibit practically the same quality of agreement/disagreement with the simulation data as those shown in figs. 2 and 3. The following conclusions may thus be drawn from our results: (i) we may claim that the results for the first-order correction to g, presented by Sane are erroneous; (ii) the first-or-
[ 1] H.J. RavachC and R.D. Mountain, Progress in liquid physics, ed. C.A. Croxton (Wiley, New York, 1978) ch. 12.
[21 H. lyetomi and P. Vashishta, Phys. Rev. A 40 (1989) 305. 131 A.R. Denton and N.W. Ashcroft, Phys. Rev. A 39 ( 1989) 426. 141 H.J. Ravachtand R.D. Mountain, J. Chem. Phys. 57 ( 1972) 3987. [5] R. Abe, Progr. Theor. Phys. 21 (1959) 421. [6] J.-P. Hansen and I.R. McDonald, Theory of simple liquids (Academic Press, New York, 1986 ) . [71 A.D.J. Haymet, S.A.Rice and W.G. Madden, J. Chem. Phys. 74 (1981) 3033. [81 Y. Rosenfeld and N.W. Ashcroft, Phys. Rev. A 20 (1979) 1208; F. Lado, SM. Foiles and N.W. Ashcroft, Phys. Rev. A 28 (1983) 2374; F. Lado, Phys. Rev. A 8 ( 1973) 2548. [91R.N. Sane, Phys. Rev. A 25 ( 1982) 1779. [IO] J.A. Barker and J.J. Monaghan, J. Chem. Phys. 36 ( 1962) 2564. [ 1I] H.Conroy, J. Chem. Phys. 47 (1967) 5307. [ 121 I. Nezbeda, J. Kolafaand S. Labik,Czech.J. Phys.B 39 (1989) 65. [ 131 I. Nezbeda and G. Kahl, Mol. Phys., to be submitted. [ 141 S.A. Rice and D.A. Young, Discussions Faraday Sot. 43 (1967) 16. [ 151 J.-P. Hansen and G. Zerah, Phys. Letters A 108 ( 1985) 277; G. Zerah and J.-P. Hansen, J. Chem. Phys. 84 (1986) 2336. [ 16) J. Talbot, J.L. Lebowitz, EM. Waisman, D. Levesque and J.-J. Weis, J. Chem. Phys. 85 (1986) 2187. [ 171 H.J. RavachC, R.D. Mountain and W.R. Street& J. Cbem. Phys. 57 ( 1972) 4999.
339
6 September 199 I
CHEMICAL PHYSICS LETTERS
First-order correction to the three-body correlation function Ivo Nezbeda ’ and Gerhard Kahl Institut ftirTheoretische Physik, Technische Universiidt Wien, Wiedner Hauptstrasse 8- 10, A- 1040 Vienna, Austria Received 25 February 1991;
in final form 12 June 1991
The first-ordercorrection to the Kirkwood superposition approximation of the triplet correlation function, based on the density expansion due to Abe, has been calculated for the high-density Lennard-Jones fluid. While qualitative agreement with computer simulation data has been obtained, substantial disagreement with older literature data and conclusions based thereon have been found.
Although there has been much interest in the triplet correlation function g, both from the practical and theoretical points ov view, it is the great complexity of this field which has so far prevented an exhaustive treatment (for a review, see e.g. ref. [ 1 ] ). Nowadays, when the determination of the pair correlation function g of simple fluids has reached a high level of sophistication and computers have a reasonable performance rate to evaluate complex expressions, the time seems to be ripe for re-investigating the traditional approaches before adopting complex new methods via the density functional theory and direct correlation functions (see refs. [ 2,3 ] and references therein). The function g may be most conveniently written as
g3(r12,r13,b3)=g3(r, s, 0 =g(rk(sMWVc
4
1)
=ih(r, s,WV6 s, t) ,
(1)
where the subscript SA denotes the Kirkwood superposition approximation. Although the SA (i.e. R= 1)was for years the most commonly used approximation for g3, it was shown that it is substantially in error at high densities and, similarly, various
’ Permanent address: Institute of Chemical Process Foundations, Czechoslovak Academy of Sciences, 16502 Prague, Czechoslovakia. 0009-2614/91/S
03.50 0 1991 Elsevier
corrections based on semi-empirical expressions of R [1,4] also failed. Rigorously, may R be expressed as [ 51 R(r,s,f)=exp[pr(r,s,t)l =exp
F p”m(r,s,1)
( i=l
,
(2)
>
where p is the number-density and T(r, s, t) is an infinite sum of “simple 123 irreducible clusters” with f-bonds cf( r) = exp [ - /.?u(r) ] - 1 is the Mayer function and u(r) is the pair potential). I1 is possible to use also h-bonds, h =g- 1, instead of the f-bonds in the diagrams of (2)) which leads to a substantial reduction in the number of diagrams contributing to T(r, s, t), but then eq. (2) is no longer a density expansion. It is easy to see that for n& 3, the terms exceed practical applicability. Eq. (2 ) may now be seen from at least two points of view: (i) One may use eq. ( 1) as a closure in the BGY-hierarchy [ 6 ] and generate g and g3 in a selfconsistent way; this route was followed by Haymet et al. [ 71 at the time when no accurate methods for determining g( r) were available. (ii) With methods now available enabling one to determine g of a simple fluid with high accuracy [ 8 1, we can use eqs. (2) and ( 1) as defining equations for g,. Although the first route turned out to be not so bad, both approaches were practically abandoned at the beginning of the 1980’s: despite enormous effort, the results obtained did not reach expectations and the
Science Publishers B.V. All rights reserved.
337
Volume183.number 5
6 September1991
CHEMICALPHYSICSLETTERS
existing discrepancies between theory and experiment prevailed. Specifically, in 1982 Sane [ 91 showed that the rigorous first-order correction to SA,
R(r,s,r)IIR,(r,s,I)~exp[PT,(r,S,f)l,
9(r) r
3
2
with 1
is not even in “qualitative agreement” with Monte Carlo data. It is the purpose of this communication to show that Sane’s results are erroneous and, consequently, his conclusions unjustified. We consider the same system and arrangements as in paper [ 91: two state points of the Lennard-Jones fluid, (po3, k,T/e)= (0.75, 1.07) and (0.85, 0.719), each for two different isosceles triangles (r, s, s). To evaluate integral (4), we used two methods: (i) the well-known expansion in Legendre polynomials (used also by Sane) [LO]; and (ii) direct integration. For the latter, we used the Conroy method [ 111 with the new panel of points obtained recently by Nezbeda et al. [ 121 (details of the method will be given in a subsequent paper [ 131). From the dependence of the results on the number of integration points and on other parameters of integration, we claim that the error of the direct integration does not exceed 0.2% and both methods yield the same results within numerical accuracies. To check further the correctness of our numerical method, we also calculated r, using the{bonds in (4) for the same state point and arrangement considered by Rice and Young [ 141 and came to full agreement with their result. The value of ‘TVdepends evidently on h, and since this function used by Sane is not available, we use two different sets of the pair correlation function to investigate also the effect of the input h on TV: one obtained from the extended (self-consistent) MSA (HMSA) [ 151 which was shown to be very accurate [ 16 ] when compared with simulation data, and the other obtained from the HNC theory [6]. In fig. 1, we show both the HMSA and HNC g(r)s for the high-density state point. These functions differ in the heights and positions of the first maximum and minimum and should allow us to see the effect of the input on 7,. In figs. 2 and 3, we compare our results, obtained from the two different inputs using the Conroy method, with the simulation data [ 171 338
\ !
00 0
\\
_’
L 1
2
3
r/o. Fig. 1. The pair-correlation function g(r) of the LennardJones system at pc~~zO.85and kBT/c=0.719, calculated from the HMSA (full line) and HNC theory (broken line).
R(r,4 J) T/U = 2.025 I-
o.5 1.0
3.0
2.0 s/a
Fig. 2. Comparison of the approximation R(r, S, s) uR, (I, s, S) =exp [pr, (r, S, s)] for the LennardJones system atpa3=0.7S and kBT/c= I .07 with the simulation data [ 171 (0) for two different r-values. Eq. (4) with HMSA h (0), HNC h(...), and (--) ref. 191.
and those of Sane [ 91, It is seen that different hs produce only marginal differences in 7,. We have to bear this in mind when discussing the large difference between our and Sane’s data. Since Sane used parameterized molecular-dynamics data (i.e. an accurate input h which, in fact, should not differ from our HMSA input), and anyhow, different hs yield only minor differences in 7,. we conclude that the discrepancies obtained must be exclusively due to numerical inaccuracies.
Volume 183.number 5
CHEMICAL PHYSICS LETTERS
JQ, s,3)
r/u = 1.525
t
0
P/U = 1.025 0.5
i
-\ ‘_’ ’ ’ ’ ’ ! ’ ’ ’ ’ ! ’ ’ ’ ’ ! ’ ’ ’ ’ I I
0.5 1.0
2.0 s/o
3.0
6 September 1991
der correction term (3) provides a considerable improvement over the SA, and for large separations it practically removes the discrepancy between the SA and exact (simulation) g3. These findings put the rigorous expansion (2) in another light, and one may justifiably expect that the inclusion of the second-order term ~~ might push the applicability of eqs. (1) and (2) into quite a high-density range. The results, following this line, will be reported in detail in a subsequent paper [ 131. This work has been supported by the osterreichische Bundesministerium fiir Wissenschaft und Forschung and by the ijsterreichische Fonds zur Fijrderung der Wissenschaftlichen Forschung under Project Number P76 18-TEC.
Fig. 3. The same as for fig. 2 for pa’=O.ES and kJ/c=O.7 19.
References The discussion of the results must be made separately with respect to (i ) the simulation data and (ii) Sane’s data. Concerning (i), one may say that, in general, the first-order term captures the basic features of the full correction term R and thus provides an improvement over the SA. Significant discrepancies remain at separations, roughly, for S/O< 1.8, while at larger separations, the agreement is even nearly quantitative; this is especially true for the lowdensity state point. For the high-density state point, the discrepancy is much larger but nevertheless the trend of the correction is in the right direction. Concerning (ii), we see that at the lower density, our results push the correction even closer to the exact result and at the high density, they considerably remove the drastic qualitative disagreement. Moreover, Sane’s curves exhibit spurious peaks and wiggles. These do not show up in our calculations and may be only an artefact of his numerical method. To be consistent with the paper by Sane [9], we do not show results for other molecular arrangements which exhibit practically the same quality of agreement/disagreement with the simulation data as those shown in figs. 2 and 3. The following conclusions may thus be drawn from our results: (i) we may claim that the results for the first-order correction to g, presented by Sane are erroneous; (ii) the first-or-
[ 1] H.J. RavachC and R.D. Mountain, Progress in liquid physics, ed. C.A. Croxton (Wiley, New York, 1978) ch. 12.
[21 H. lyetomi and P. Vashishta, Phys. Rev. A 40 (1989) 305. 131 A.R. Denton and N.W. Ashcroft, Phys. Rev. A 39 ( 1989) 426. 141 H.J. Ravachtand R.D. Mountain, J. Chem. Phys. 57 ( 1972) 3987. [5] R. Abe, Progr. Theor. Phys. 21 (1959) 421. [6] J.-P. Hansen and I.R. McDonald, Theory of simple liquids (Academic Press, New York, 1986 ) . [71 A.D.J. Haymet, S.A.Rice and W.G. Madden, J. Chem. Phys. 74 (1981) 3033. [81 Y. Rosenfeld and N.W. Ashcroft, Phys. Rev. A 20 (1979) 1208; F. Lado, SM. Foiles and N.W. Ashcroft, Phys. Rev. A 28 (1983) 2374; F. Lado, Phys. Rev. A 8 ( 1973) 2548. [91R.N. Sane, Phys. Rev. A 25 ( 1982) 1779. [IO] J.A. Barker and J.J. Monaghan, J. Chem. Phys. 36 ( 1962) 2564. [ 1I] H.Conroy, J. Chem. Phys. 47 (1967) 5307. [ 121 I. Nezbeda, J. Kolafaand S. Labik,Czech.J. Phys.B 39 (1989) 65. [ 131 I. Nezbeda and G. Kahl, Mol. Phys., to be submitted. [ 141 S.A. Rice and D.A. Young, Discussions Faraday Sot. 43 (1967) 16. [ 151 J.-P. Hansen and G. Zerah, Phys. Letters A 108 ( 1985) 277; G. Zerah and J.-P. Hansen, J. Chem. Phys. 84 (1986) 2336. [ 16) J. Talbot, J.L. Lebowitz, EM. Waisman, D. Levesque and J.-J. Weis, J. Chem. Phys. 85 (1986) 2187. [ 171 H.J. RavachC, R.D. Mountain and W.R. Street& J. Cbem. Phys. 57 ( 1972) 4999.
339