First order quantum correction to the free energy of simple liquids
S K D A T T A
Department of Physics, Jalpaiguri Govt. Engg. College, Jalpaiguri 735 102, India MS received 14 November 1984; revised 10 June t985
Abstract. Using the perturbation theory of Weeks et ai the first order quantum correction to the free energy of a simple fluid characterised by a double Yukawa potential function has been expressed in a simple closed analytical form which allows numerical calculation simply on a desk calculator.
Keywords. Q u a n t u m correction; double Yukawa potential; radial distribution function;
perturbation theory.
PACS No. 61-20
1. Introduction
The derivation o f equilibrium properties o f simple dense fluids in closed analytical form has gained much importance after the work of Verlet and Weis (1972a, b). Recently the thermodynamics a n d structure o f simple fluids in a fully analytic form have been derived (Datta 1983) where the perturbation theory o f Weeks et al (1971) (hereafter termed WCA) has been utilised and the simple fluids characterized by a double Yukawa (bY) potential function (suggested by Foiles and Ashcroft 1981) with parameters adjusted to fit closely a Lennard-Jones (u) potential function. The DV potential is given by
t~)Dy (1 r) = ~ {exp [ -- a(r/a -- 1)] -- exp [ -- b(r/tr - 1)] }, (1) where E = 2.0199e, a = 14.735, b = 2.6793, and e and tr are the parameters o f the LJ potential.
In the present paper we derive a closed-form analytical expression for the first order quantum correction to the free energy o f a simple fluid characterised by the DV potential function utilising the WCA perturbation theory. It is well k n o w n that the h 2 correction to the free energy for an analytic potential is given by
h2p f o V2q~(r~Cl(r)r2 dr + . . . . (2)
flf = [3f'! + ( k a T ) 2 24rim
where m is the mass o f an atom and ~Z(r) is the radial distribution function (RDF) o f the classical system (i.e. h --* O, m --* ~ ) fl = I / k a T , kB the Boltzmann constant and T the absolute temperature.
335
In WCA perturbation theory the pair potential 4'(r) is divided as
4'(r) = 4'o(r) + 4'1(r), (3)
where 4'o(r) and 4'1 (r) are the reference (repulsive) part and the perturbation (attractive) part respectively o f the pair potential. These are given by
4'o(r) = 4'(r)+ e r ~< rm,
= 0 r > r,,, 4't (r) = - e r ~< r . ,
= 4'(r) r > r . , r m --~ 2 1 / 6 0 -.
(4)
(5)
We now proceed to calculate separately the contributions o f the reference part and the attractive part o f the potential to the integral in (2).
2. Calculation of the reference part contribution
To calculate the contribution o f the reference part to the integral in (2) we approximate the ROF of the actual fluid by that of the reference fluid system
i.e.
we assume~t(r) --* g0(r) which is valid at the high density range o f the fluid. With this choice the form of the integral in (2) may be rewritten as
lrcf = ~o V2d~o(r)o"(r)r2 dr ~ ~o V24'o(r)Oo(r)r2 dr,
(6) Again the WCA approximation for the RDF gO(/') is given bygo(r) = Yd(r)
exp [ -- fl4'o(r)], (7)where
Ya(r)
is the reduced pair correlation function o f a hard sphere fluid o f diameter d.It is defined as
Yd(r) = gd(r)
exp [fl4'a(r)], (8)gd(r) being the hard sphere g(r) and d is chosen so that
~ o Yd(r){exp [ - fl4'o(r)] - exp [ - fl4'~(r)] } dr = O. (9) We note that
d z 4,0 1 d dd~.
dr 2 exp [ - f14'o] = - -~-rr [ 6 4 ' o ( r ) ] f l - ~ t~o(r) (10) where the sharp peaked function
t~,o(r) ---- d {exp [ - fl4'otr)] } (11)
behaves like a delta function around r = d. The integral in (6) may then be separated as
lrcf = f o (d24'° +2d--ff-4'd° 2 r )Yd(r)exp [ -
2d d 2
fl f : (r/d) YAr/d)6+,o(r)dr--~ f : (r/d)~ YAr/d) d ['~¢o (r)]dr
- d 2 f ; (r/d)2 yJ(r/d)~r° 64~o (r)dr.
W e n o w e x p a n d
(r/d) Yd(r/d)
a n d(r/d) 2 Yd(r/d)
a r o u n d r = d as(r/d) Ya(r/d) = ~0 + oq (r/d -
1) +~-(r/d -
l) 2 + I I Oz
(12)
(13)
(r/d) 2 Ya(r/d) = tTo + ax (r/d - I) + ~ (r/d - I)2 + . . . (14) the Percus a n d Yevick 1958 (hereafter termed PV) approximation the
where
dn = J o { 1 - exp [ - fl~bo(r)] } dr. (19)
I n the a b o v e d e r i v a t i o n we have a s s u m e d t h a t all terms i n v o l v i n g 6", n > / 2 are negligibly small. T h e c o n t r i b u t i o n o f the last t e r m in (17) involving & is h o w e v e r extremely small.
where in
p a r a m e t e r s ao, ~ti, or2, ao, a l a n d a2 are given by
1 + q/2
1 -- 5q - 5q 2 -- 3q(2 -- 4q - 7q 2) ao - (1 - q ) 2 , 0tt -- (1 -- q)3 ' ~t2 = (1 -- q)+Oo - °to, al = ~tt +~to, a2 = 2ctx + a 2 . (15) W e again e x p a n d d4~o/dr in a T a y l o r series a r o u n d its value at r = d as
d4~o - gp'o(d) + (r - d)~'~(d) + ½ (r -
d) 2 ~b~'(d) + . . . . (16) d rwhere q~(d), q~g(d) a n d q~g'(d) are respectively the first, s e c o n d a n d t h e third derivatives o f 4>o(r) at r = d. Finally with the h a r d sphere d i a m e t e r d d e t e r m i n e d by the wcA criterion (9), using the e x p a n s i o n s o f
(r/d) Ya(r/d)
a n d(r/d) 2 Ya(r/d)
a n d utilising the properties o f delta f u n c t i o n it can be very easily s h o w n t h a t the integral in (12) reduces to,17,
where On (d) a n d
O'a(d )
are respectively the value a n d t h e derivative o f h a r d sphereon(r)
e v a l u a t e d at r slightly g r e a t e r t h a n d. 6 the function o f t e m p e r a t u r e a l o n e is defined as (Verlet a n d Weis 1972a).
r+:, ), fo°( , )"
= d.Jo \ d - I ,S+o(r)dr ~ ~ - 1 64~o(r)dr, (18)
In the av approximation the values o f on(d) and O'd(d) are given by (Lebowitz 1964)
1 +q/2 9q(1
+ q )#a(d) = tro = (1 - r / ) 2 and O'n(d) 2 ( 1 - q ) 3 (20) The value o f on(d) can however be correctly determined from the semi-empirical hard sphere equation o f state obtained by Carnahan and Starling (1969) while the values o f O'4(d), a~ and a2 can be accurately obtained from the semi-empirical hard sphere ROE derived by vw (1972a, b).
3. Calculation of the attractive contribution
The c o n t r i b u t i o n o f the attractive part o f the DY potential to the integral in (2) can be readily expressed in a closed form manner. Using the fact that Ya(r) = g~(r) for r > d, our desired integral reduces to [Hansen and M c D o n a l d 1976]
latt = E (r/o)ga(r) { a 2 exp [ - a(r/a -- 1)] -- b 2 exp [ - b(r/a - 1)] } dr
m
(21)
and in the P¥ approximation this can be expressed as
latt = E~rc2[ a2 exp [a] {G(ac) - F(a)} - - b 2 exp [b] {G(bc) - F(b)} ], (22) where G (z) is the Laplace transform o f r#~ v (r) already derived by Wertheim (1964) and is given by
zL(z) (23)
G(z) = 12r/{L(z) + S(z) exp [z] } '
L(z) = [(1 + q/2)z + (1 + 2n)]12q, (24)
S(z) = (1 - q)2z3 +6q(1 -rl)z 2 + 181'/27, - - 12q(1 + 2r/). (25) Again F(z) is given by
exp
~c ~ 2 - J exp [ -
zc]{x2 2xm 2 ~ - zcxm]] (26)
- ~ , - , + - ~ c + ~ ) e x p [
where x,, = rm/d, c = d/a, d(p, T) is the WCA diameter.
A = % - ~q + 2 ' B = ~t I - ~ 2 and C = ~t2/2.
(27)
However to get the correct value o f the integral in (21) we can utilise the vw expression for od(r) as given by Jedrijezek and Mansoori (1979).
Finally, the h z c o r r e c t i o n t e r m is given by h 2
(ks T) 2 241rm (Iref d- latt)- (28)
Using t h e wcA p e r t u r b a t i o n t h e o r y we are thus able to o b t a i n the first o r d e r q u a n t u m c o r r e c t i o n to the free e n e r g y o f the simple liquid in a fully analytic f o r m w h i c h a l l o w s numerical c a l c u l a t i o n s i m p l y o n a desk calculator.
M a n d e l (1972) calculated the integral in (2) using the LJ p o t e n t i a l a n d t h e WCA p r e s c r i p t i o n a n d f o u n d t h a t the WCA f o r m a l i s m w o r k e d surprisingly well. Since t h e earlier w o r k was d o n e f o r the LJ potential, utilising the present closed f o r m a n a l y t i c a l expression t h e values o f L a p l a c i a n ( V2q~ ) * (defined as
tr-- ~ ( V2~b = 47tp V2dp(r)ga(r)r2dr)
have been calculated using b o t h the o v a n d LJ potentials. All the c a l c u l a t i o n s in this p a p e r h a v e been d o n e using Pv values o f gd(r, r/). T h e calculated values a r e s h o w n in table 1 a l o n g with the results o f M a n d e l a n d the m o l e c u l a r d y n a m i c s (MO) results o f Verlet (1967). T h e c o n t r i b u t i o n o f the attractive part o f t h e LJ p o t e n t i a l t o t h e integral in (2) being difficult t o express analytically has been a p p r o x i m a t e d by the DV a t t r a c t i v e part c o n t r i b u t i o n . This r e p l a c e m e n t h o w e v e r entails very small error.
T h e wcA d i a m e t e r s have been calculated in the way suggested b y M e y e r et al (1980) using PV values o f Y given by
y = d In gd(r, r/) I 9r/(l +r/)
- ~ n r I , = , + o = - 2 ( 1 - ~/) (1 - v//2) (29)
Table 1. Values of the Laplacian ( V2~b )* as a function of density and temperature.
Present work Work of Mandel
p* T* DY-WCA LJ-WCA LJ-WCA LJ-MD
ff65 0-900 501 491 496 499
0-65 1.036 542 531 535 552
0-75 0-827 670 658 661 634
0-75 0-881 693 679 682 677
0-75 1.071 771 753 752 739
0-85 0-658 833 820 824 787
0-85 0-719 869 854 856 816
0-85 0-786 907 891 891 854
0-85 0-880 961 941 937 902
0-85 1-128 1093 1065 1055 1032
0-85 1-214 1136 1105 1078 1075
0-85 2.202 1585 1523 1486 1483
0-88 0-940 1095 1070 1062 1027
0-88 1.095 1184 1153 1140 1118
MD, molecular dynamics
The values o f 6 have been taken from the work o f v w (1972a, interpolated from table 12 contained therein). Again the same d and 6 (calculated for LJ potential) have been utilised also for the DV potential calculations.
An inspection o f table 1 shows that the present results for the u potential are almost exactly identical to that o f Mandel. Some deviations between the two sets o f results appear with the rise o f temperature. This is evident because with the rise o f temperature the contributions o f the terms involving
and o f higher order terms also become significant. However, most o f the deviations between the present result and that o f Mandel can be accounted for by including only the first order term in el. This term can be obtained in the following way. I f we introduce one more term in the vw expression for WCA diameter as
d = d n l + ~ z - ~ + el , (30)
0
then the term involving the first power o f el is given by
and this must be added to Ire r. Since the quantity ~ " ( d ) is o f o r d e r 103 while el is o f order 10- ~ and increases with temperature this e~ term will contribute significantly at the higher temperatures.
N o w coming to the case o f small disagreement between the LJ and DV results for ( V2~b >* we note that this happens due to the slight deviation between the values o f the higher order derivatives (~b~'(d), etc) o f the reference parts o f the two potentials.
Finally, considering the overall simplicity o f the m e t h o d we conclude that the present analytical route is the most easiest one for estimating the first order q u a n t u m correction to the free energy o f simple liquids.
References
Carnahan N F and Starling K E 1969 J. Chem. Phys. 51 635 Datta S K 1983 Pramana 20 313
Foiles S M and Ashcroft N W 1981 J. Chem. Phys. 75 3594
Hansen J P and McDonald I R 1976 Theory of simple liquids (New York: Academic Press) pp. 184 Jedrijezek C and Mansoori G A 1979 Acta Phys. Pol. A56 583
Lebowitz J L 1964 Phys. Rev. A!33 895 Mandel F 1972 J. Chem. Phys. 57 3929
Meyer A, Silbert M and Young W H 1980 Chem. Phys. 49 147 Percus J K and Yevick G J 1958 Phys. Rev. I t 0 l
Verier L and Weis J J 1972a Phys. Rev. A5 939 Verlet L and Weis J J 1972b Mol. Phys. 24 1013 Verlet L 1967 Phys. Rev. 159 98
Wertheim M S 1964 J. Math. Phys. 5 643
Weeks J D, Chandler D and Andersen H C 1971 J. Chem. Phys. 54 5237
