- Email: [email protected]
Perturbation theory in terms of electron density
11 November 1996 PHYSICS LETTERS A
Physics Letters A 222 ( 1996) 3 1S-323
ELSEVIER
Perturbation theory in terms of electron density Manoj K. Harbola, Amp Banerjee Laser Programme, Centre for Advanced Technology, Indore 452013, India
Received 25 June 1996; revised manuscript received 3 September 1996; accepted for publication 6 September 1996 Communicated by B. Fricke
Abstract When an inhomogeneous electron gas is subjected to a perturbation, its energy and density both change by small amounts. We calculate the changes in the energy explicitly in terms of the density changes within the density-functional theory of many-electron systems. We also derive the equations for the induced densities, and using these show that a density correct up to order n in terms of the perturbation parameter is sufficient to give energy changes up to order (2n + 1) . As a corollary, the even-order energy changes E(“‘) are variational with respect to the density changes p(“). The equations for the induced densities also follow from this corollary. The even-order corollary also gives a variational method of calculating the induced densities. The theory is demonstrated by applying it to calculate the polarizability and hyperpolarizability of the hydrogen, helium and neon atoms.
When a system of N electrons throughout the paper)
Ho
=&o: +uo(ri)l i=l
+
described
$c (izj)
by the Hamiltonian
(atomic
units,
lel = m, = Fi = 1, are used
-i ‘ij
is subjected to an applied perturbative field Au(l), the changes in its properties are also small and as such can be described by an expansion in powers of the perturbative parameter A. Thus the stationary state wavefunction !P and the ground-state energy E for the Hamiltonian HO + Au(‘) are written as [l] !P = To +
c
AW’)
(2)
i=l
and E = E. +
F A'E('),
(3)
i=l
where PO and EO are the wavefunction and the energy, respectively, corresponding to the Hamiltonian I$ (1). The corrections P(‘),P(2), . . . to the wavefunction and EC’), EC*), . . . to the energy are calculated 0375~9601/%/$12.00 Copyright 0 1996 Published by Elsevier Science B.V. All rights reserved. PII SO375-9601(96)00675-S
of by
316
M.K. Ha&la.
A. Banerjee/Physics
Lerters A 222 (1996) 315-323
solving the perturbation equations [ 11 to the appropriate order. Interestingly, a wavefunction correct only up to order n is required [ 21 to calculate the energy up to order (2n + 1) . This is the (2n + 1) theorem, The theorem is valid not only for the exact solutions of the Schriidinger equation but also [ 31 for approximate methods such as Hartree-Fock (see for example Ref. [ 41) theory. The purpose of this Letter is to develop perturbation theory in terms of the density within density-functional theory (DFT) [ 51 of many-electron systems. Thus we derive: (i) the energies to various orders in terms of the induced densities to appropriate orders, and (ii) the equations for the induced densities to each order, Using these equations, we show that the density correct to order n in terms of the perturbation parameter is sufficient to give energy to order (2n + 1) ( (2n + 1) theorem), as expected [ 61 from the variational nature of the energy functional with respect to the electronic density. As a consequence of this, the even-order corollary of the theorem follows [ 61: if the density is known exactly to order (n - 1), the even order energy change E(2n) is variational with respect to the density change p (‘). This is then shown to lead variationally to the equations for induced densities to each order. Further, this also gives a variational method of determining the induced densities. We demonstrate the theory developed with the examples of the hydrogen, helium and neon atoms. Before presenting our work, we point out that the theory developed below is distinct from the well established perturbative scheme [ 7,8] within the Kohn-Sham approach [ 51 of DFT. The Kohn-Sham method is an orbitalbased approach and as such its methodology [ 7,8] is similar to that of Hartree-Fock theory [ 31. On the other hand we do everything in terms of the density, obviating thereby the need to invoke single-particle orbitals. In DFT the ground-state energy of a system of electrons is expressed as a functional E[p] of the density p(r). This functional is made up of two components,
UpI = Npl +
s
uext(r>dr>drv
(4)
where F[ p] is a universal functional of p(r) and comprises the kinetic and interaction energy of the electrons; bxt is the external potential electrons are moving in. The energy is shown [9] to be variational with respect to p(r), and therefore the ground-state density is obtained by solving the following Euler equation,
SE[PI --p=o.
(5)
Q(r)
Once the ground-state density is known, the corresponding energy is given by the functional E[ p] . In Eq. (5) p, identified as the chemical potential [lo] (also see the discussion in Ref. [5] ), is the Lagrange multiplier to ensure that the total number of electrons is conserved. Now for the system with Hamiltonian Ho (Eq. ( 1) ), with ground state energy Eo, density pa(r) and the chemical potential ,UQ,Eqs. (4) and (5) can be written as
Eobol = Rpol +
J
uo(r>po(r>
dr,
(6)
and
SFbol G(r)
+ uo(r) = piI,
where SF [ po] /6p( r) means that the functional derivative SF/Sp is evaluated at p(r) will be followed throughout the paper. When the perturbation Au(‘) is applied, the density changes to po + Ap with
J
Ap( r) dr = 0,
(7) = po. The same notation
(8)
M.K. Harbola, A. Banerjee/Physics
and the chemical
Etpo+bl
potential
Letters A 222 (19%) 315-323
317
to m + A/L. From Eq. (4), we have
=F[po+Ap]
+
which can be Taylor expanded
(uo+~~(l’)(Po+Ap)dr,
(9)
J
to get
+S~opodr+/L’aapd’+~/~“)podr+~/~“)A~dr. Using Eqs. (6) -( 8) this can be rewritten as
E[pO+ Ap] = Eo[Po]+ g ; n=2
+ A Similarly
J
u”‘Podr+
A
.
1 Spo~~;r~~MISp(T . .Ahp(r,) n fph)Ap(r2). ...
dr, dr2..
. dr,
J
(11)
#)Apdr.
from Eq. (5), /-LO+ Ap can be written as
@‘[PO + API + f@(r) + /G(r) /LQLO+A/X= G(r)
=%‘+%
n=l
.
S("+')F[po] Adm> drl Wrl) ...Q(r,) J~p~r~ap~r~~
dr, +ua(r)
+ Au”)(r). (12)
Since Ap(r) Ap(r)
and Ap can also be expanded perturbatively,
we write
= xA’p(‘)(r),
(13)
and Ap
=
c
A’,&‘,
(14)
i
with
J#j(r) dr=O
(15)
for each i. In the following we first demonstrate with the examples of E(l), . . . , EC51 that the induced density only up to order two is required to calculate energy up to the fifth order. We then give the general proof for the (2n + 1) theorem. From Eqs. (11) and (12) expressions for EC’), . . . ,Ec5) and ~(~),p(~) are E”’ =
Ec2’ =
J J
u(‘)(rl)pO(rl)
dl)(rl)p(‘)(rl)
drl,
dr, + 12
(16)
J
~c2)mol
G(nhMr2)
p(1)(r,)p(1)(r2) dr, dr2,
(17)
M.K. Ha&la.
318
Ec3’=
u”)(r,)p’2)(r,)
s
S’3’F[pOl
+t
Q(r1
p(1)(r,)p(2)(r2) drt dr2
)G(r2)
p(‘)(r,)p”)(r2)p(‘)(r3)
dr, dr2,
(18)
Sd~l)~~(~2)&'(~3)
s
s
Letters A 222 (1996) 315-323
~(2’FbOl
dr, +
s
Ec4)=
A. Banerjee/Physics
up
dr, + L 2
~‘2’Fbol
+
~c2mkd s
dr, dr2
p’2)(q)p(2)(r2)
&J(Q)&d~2)
p(‘)(r1)p(~)(r2) drt dr2
&J(r1h%(r2)
s
~(3)Fbol
+; J I
+2;i Ec5’=
+(rl
p(‘)(r’)p(1)(r2)p(2)(r3)
s
6(4’F[Pol
p(‘)(r’)p(‘)(r2)p(l)(r3)p(‘)(r4)
6p(rl)Sp(r2)6p(r3)Sp(r4)
u”‘(r’)pc4)(r’)
SC2’F[pol
dr’ +
J
.I
pc1)(r’)pc4)(r2)
J J ~'3'Fbol +; J +;J J ~(3’Fbol
+;
drl drzdrgdr4,
(19)
dr’dr2
~p(r’)G(r2)
Sc2’Fbol Sp~rl~Gp~r2~pc2)(n )d3’(rd
+
dr’drzdq
)&(r2)+(r3)
drl dr2
p(‘)(r’)p(‘)(r2)p(3)(r3)
dr’ drzdr3
p(2)(r’)p’2’(r2)p(‘)(r3)
dr’dr;!dq
&(rl)Mr2Vp(r3)
Wrlb%(r2>Mr3)
~(4'F[pol
p(‘)(r’)p(1)(r2)p(1)(r3)p(2)(r4)
dr’ drzdrsdr4
Sp(rl)Sp(r2)6p(r3)Sp(r4)
+m
~‘5’F[pol
1
Sp(rl
)sp(rz)sp(r3>sp(r4)sp(r5)
XP (‘)(r’)p(1)(r2)p(1)(r3)p(1)(r4)p(1)(r5)
dr’drzdq
dr4dr5
(20)
and
p(I) = u(‘)(q) #I
+
J
a2F[pol p(‘)(r1) dn, Jjp(r)~p(rl>
a2F[pol
=
J
ap(r
pc2)(r’)
dr’ + i
J
(21) ~3F[~ol
pt1)(r’)p(‘)(r2)
dr’dr2.
(22)
G(r)MrI)&(r2)
Having been derived from the Euler equation (Q. (5) ) for the density, Eqs. (21) and (22) are the equations for the induced densities p(l) and pC2), respectively. Later in this paper we also derive them variationally from the even order corollary of the (2n + 1) theorem. Using Eqs. ( 15)) (21) and (22) expressions for E”) , . . . can be simplified (E(l) and Ec2) remain the same) and written as E’3’
= t
I
St3’F[ PO] ~p(n)~p(r2)@P(r3)
p(‘)(r’)p(‘)(r2)p(‘)(q)
drl dr2dr3,
(23)
M.K. Ha&la,
E(4)
2
E’5’
1
319
drt dr2
p(2)(r,)p(2)(r2) Q(u)Mr2)
s
J J
~(3’wcll
+(rl
fz
Letrers A 222 (1996) 315-323
6(2’F[po]
= I
+;
A. Banerjee/Physics
drr dr2drg
p’1’(rl)p”‘(r2)p’2’(~3)
)+(r2)&‘(r3)
~(4mPol
drr dr2dr3drq,
~~‘~(~1)~~‘~(~2)~~‘~(~3)p(l~(~4)
‘%(rl)+(r2)+(r3)&‘(r4)
J +i J S'4'F[Pol J
(24)
~c3m%l
= 1 2
S.(rl)Gp(r2)6p(r3)
d2)(rl
)d2b2)d1)(r3)
drl
dr2
dr3
drr drzdrsdrq
p’11(~,)p”‘(~2)pf1f(~3)pf2’(~4)
&drl)+(r2)@(r3)&‘(r4)
+m
@5’mol
1
+(n
)~~(r2)@(~3)~A~4)~~(~5)
(25)
Note that for the applied perturbation considered here (first order in A), energy corrections from Ec3) onwards are calculated from the universal functional F [ p] alone. It is evident from these expressions that p(t) alone gives an energy up to Ec3) and p(*) up to E (5) . This indicates the existence of the (2n + 1) theorem. We now prove it in general. By writing out the full expression for E (2n+1) from Eq. ( 1 I), and making use of the expression for p (n) from Eq. (12), we show that all the terms in the expression for Ec2”+*) involving induced densities of the order greater than n integrate to zero. From Eq. ( 11) the general expression for Ec2”+‘) is
,$2n+‘)
=
J X
u(1)(r,)p(2”)(rl) drt
J
Similarly
kl
@+‘Fbol
G(r1
=k 2 j=l
X
J
2 ili2...ik+l=l
pyq)p(iqr2)
. . . Mn+1)
)@(f2)
the expression
p(n)
+ 2
S(ir + i2 +
. . . +ik+l - (2n+ (k+
1))
l)!
. . . p(“+‘) (rk+l > drt dr;! . . . drk+, .
(26)
for p (n) (n > 1 )can be written from Eq. ( 12) and is given as
&it
+htl;.+ij-n)
ili2...ij=l
@+'Fbol aP(rl)aP(r2)
. . . aP(rj+l)
p(i’)(r2)p(i2)(T3)
. . . pciJ)(rj+l)
dr2dr3..
. drj+r.
(27)
As discussed above, Eq. (27) is the equation for the induced density p (“). In Eq. (26) the largest value of k for which the expression will have induced densities of order (n + I) and above is n; all the lower values of k would involve higher order induced densities. We therefore partition the k sum as follows,
320
M.K. Ha&da.
pnf’)
=
u(1)(r,)p(2”)(rl) dr,
J X
J
+ 2 5 k=l iliz...it+l=l
~k+lFbol
~P(~l
)Sp(r2)
Letters A 222 (19%) 315-323
6(i, + i2 + . . .+&+I
p(‘qr*)p(iqr2)
’ . .8P(Tk+, 8(i,
+F
(k+
1))
- (2n+
l)!
. . .pci’+‘)(rk+,)
dr, dr2..
. drk+,
. . .#i”+‘)(rk+,)
dr, dr2..
. drk+,.
)
+&+...+ik+l
-(2n+l))
2 k=(n+l)
X
A. Banerjee/Physics
J
(k+
i,iz...it+l=l
~k+lmol
&‘(rl
)@(r2)
. . . &'('.k+l
)
l)!
#&l) (r, )p(i2) (Q)
The third term in the above equation (28) involves pCn, and lower order induced densities. We now separate the highest order induced densities p(*“, and P(~‘+~, (m = 1, . . . , (n - 1) ) in the first two terms and write Eq. (28) as
+5
J
drlp(2”-“t)
[@*+” (m+‘) c c (k+
k=l
llt=l
X
(rl)
l)6(il +“;;;l;;“‘l))
iliz...ik=l
J~~~~~~~~~2~!po~,,,+,~ di”(r2)di2)(rd..
+ (terms involving
.#‘)(rk+d
ptn, and lower order induced densities).
dQ..
. drk+,]
(29)
Notice that once p(2”-“‘) is taken out, the maximum that k or il, i2,. . . can go up to is (m -t1) because of the delta function. The factor (k + 1) arises because p(2”-m) can be chosen in that many different ways. The expressions in the square brackets above, however, are the expressions for I_c(‘, (IQ. (21)) and P(~+,) (Eq. (27) ), respectively. Thus by the normalization condition (Eq. ( 15)) the integral s dr,p(2”-m) ( r,)p(m+‘) is zero. The highest order of the induced density in the terms left is therefore n. One can similarly show that the next order correction to the energy, i.e., E (2n+2), does require the induced density p(“+‘,. This is because for this energy correction one of the terms involving p (‘+l) is quadratic (the rest of the terms with p(“+*) are linear in it and have products of induced densities of order n or less as factors). This term comes from the k = 1 term in the sum and has the factor 1/2k! (= $) which is not the right factor ( l/k! = 1) for pCn+‘, to be pulled out of it. Thus these terms do not integrate to zero. This completes the proof of the (2n + 1) theorem. As mentioned in the beginning, the results above follow [6] from the variational nature of the energy functional. From Eq. (5) it is evident that any number conserving changes in the density would give an energy change which is quadratic in Ap. Thus for a trial density p’, there exists a number K such that [ 6,111
0 6 E[dl
- E[pl G Klb - ~11~7
(30)
where p is the exact density. Therefore, if the trial density is correct up to order n, the error in the energy is of the order (2n + 2). Thus for this density the energy is determined correctly up to order (2n + 1). From Eq. (30), it is also clear [ 61 that E(2n+2) is variational with respect to a trial p(“+‘) if the density is known exactly up to order n. It is also evident from the discussion above that the even order energy E(2n+2) has one term which is quadratic along with others which are linear in p(‘+‘). The variational nature of the even order energy terms leads to the Euler equations (27) for the induced densities p(“) with pCn) as the Lagrange multiplier to satisfy Eq. ( 15). Thus by minimizing Ef2) with respect to p(‘,, Eq. (21) for p(l) is obtained; ,u(‘, is the corresponding Lagrange multiplier. Knowledge of p(l) gives
M.K. Ha&la,
A. Banerjee / Physics Letters A 222 (19%) 315-323
321
energy correct up to order Ec3). Now p(l) is kept fixed and Eq. (22) for p’*’ is obtained by minimizing Ec4) with respect to p(*), with PC*) as Lagrange multiplier. In this manner one can go on to obtain higher order corrections both to the density as well as the energy. These equations are DFT counterparts of the more familiar equations 1l] for wavefunction corrections to various orders. As pointed out before, the theory above is different from the orbital-based perturbation scheme within the Kohn-Sham DFT methods. The difference between the two approaches is quite evident when the expressions for the energy changes to various orders in the two schemes are compared. As shown above, in the present method energy corrections from EC3) onwards are given by the universal functional F[ p] alone. In the KohnSham method, on the other hand, energy changes to all orders involve [ 7,8 ] the full Hamiltonian including the applied perturbation. In the past we have calculated [ 121 polarizabilities (energy up to order EC*)) of atoms and metal clusters variationally by employing the expression for EC’) given by Eq. ( 17). Bartolotti [ 131 has also mentioned the (2n + 1) theorem in terms of the density and has calculated polarizabilities of atoms using Eq. ( 17). However, it is only now that density based perturbation theory has been fully developed. This makes it possible to calculate higher order response functions of a system in terms of its density. As a demonstration we apply the theory presented above to calculate variationally the hyperpolarizabilities (energies up to the fourth order) for the hydrogen, helium and neon atoms subjected to a static electric field. Induced densities for the hydrogen atom in a static electric field are known [ 14,151 exactly, and are given as an
p(‘)
=
(2r + r*) cosBp0,
p ~*‘~[($~*+~~~+~~~)C0s~e+(~r~+~r3-~)1p0.
(31)
Taking a cue from Eq. (27), the variational form for the induced densities we choose is p”‘(r)
=Al(r)cosBpo(r),
pC2)(r) = [AZ(~) +43(r)
Al(r)=(aIr+b,r2+c,r3+dIr4+e,r5), 43(r)
=
(a3r
+
by* + c3r3 + d3r4 + e$),
cos2Blpo(r)
AZ(r) = (u2r +
b2r2
+
(32)
c2r3 + d2r4 + e2r5),
(33)
with al,. . . , ei, a*, . . . , e2 and ag, . . . , es being the variational parameters. The universal functional F[ p] for the hydrogen atom represents the kinetic energy of the atom, and is given exactly by the van Weizsacker functional [ 161 Tw[Pl = $
I
IVp(r) p(r)
I2 dr
(34)
For the helium atom we perform our calculation within the Hartree-Fock theory so that the universal functional F[ p] has two terms additional to the kinetic energy (given by Eq. (21) ) . These are the Hartree energy, EHbl
=
IS
$
(35)
and the exchange energy E, [p] , which for the helium-like atoms is given exactly as - $EH [ p] . The ground-state density of helium is obtained from its analytic Hat-tree-Fock wavefunction 1171. As described above, variational calculations to determine p(‘) and p’*) are performed in a sequential manner [ 151. First EC2’ is minimized to obtain p(l). Then keeping p(l) fixed p(*) is varied to minimize Ec4). This then is the density based counterpart of the similar wavefunctional [ 181 (also see Ref. [3]) or Kohn-Sham orbital method [8,15]. When this procedure is applied to the problem [ 141 of the hydrogen atom in a static field, it leads to the correct densities, as given by Rqs. (31), and the exact energies [ 141 EC*) = -$ a.u. and
M.K. Harbola, A. Banerjee/Physics
322
Letiers A 222 (1996) 315-323
E(4)
= -F a.u. for it (E(l) = Ec3’ = Ec5) = 0 by symmetry). Notice that Ec4) is calculated from the kinetic energy functional alone. Thus all the information about the applied perturbation is carried to EC41 by the induced densities. This, as discussed above, is in contrast to the Kohn-Sham method [ 8,151 where the expectation value with respect to the induced orbitals involves the applied perturbation also. For the helium atom the second order energy is -0.6611 a.u. This gives the polarizability to be 1.322 a.u. which is equal to the Hartree-Fock theory result [ 15,191. The fourth order energy is -1.5008 a.u. which gives the hyperpolarizability to be 36.02. This is close to but slightly lower (because of the approximate form chosen for the induced densities) than the more accurate Hartree-Fock value of 37.77 a.u. [ 151. In the calculations above, the kinetic energy and the exchange energy are treated exactly. Thus it is not clear how accurate a density-based calculation would be if these energies are to be approximated. For this reason, we have also performed the calculations within exchange-only for the neon atom with the kinetic energy given by Eq. (34) and an exchange energy by the local density approximation [5] as 113 E,LDA[P]
= -t
a 0
p4’3 (r)
dr.
(36)
s
The reasons for choosing the von Weizsacker functional over the other kinetic energy functionals such as the Thomas-Fermi [ 51 or the gradient-expansion approximation [ 51 functionals are discussed elsewhere [ 121 in the context of polarizability calculations. More discussion will follow in our subsequent work. When HartreeFock densities [ 171 are employed to calculate the response properties of neon, they give cy = 3.1 a.u. and y = 74.15 a.u. which compare well with the coupled Hartree-Fock values [ 191 of 2.34 a.u and 71.9 a.u., respectively. This, along with the calculation for hydrogen and helium, demonstrates the applicability of the density-based perturbation theory, the (2n + 1) theorem and its even-order corollary to orders higher than two. These calculations also show that although in principle exact knowledge of induced densities to order (n - 1) is required to determine p(“) variationally, in practice a sequential variational procedure also leads to stable and accurate results. For a further discussion in terms of Kohn-Sham orbitals, we refer the reader to Ref. [ 151. To conclude, we have developed a density based perturbation theory and shown that in the perturbative regime, ground-state density accurate to order n is sufficient to determine energy up to order (2n + 1) . Further, the even-order energy changes EC*“) are variational with respect to the nth order density change p(“). This then gives the method of calculating the induced density, and leads to equations that determine them. These results would make perturbative calculations within DFT quite easy as compared to orbital based methods, as demonstrated in Refs. [ 12,131 for second-order calculations. We have demonstrated the applicability of the theory by employing it to calculate linear and nonlinear response properties of hydrogen, helium and neon atoms. Work in the direction of application of the theory to other atoms is in progress and will be reported in the future. We thank Dr. Kailash Rustagi for his constructive
comments.
References [ I ] L.I. Schiff, Quantum mechanics (McGraw-Hill, New York, 1969). [2] A. Dalgarno and A.L. Stewart, Proc. R. Sot. A 238 (1956) 269; PO. Lcwdin, J. Molec. Spectr. 13 ( 1964) 326. [3] PW. Langhoff, J.D. Lyons and RI? Hurst, Phys. Rev. 148 (1966) 18. [4] CF. Fischer, The Hartree-Fock method for atoms (Wiley, New York, 1977). [5] R.G. Parr and W. Yang, Density-functional theory of atoms and molecules (Oxford Univ. Press, Oxford, R.M. Dreizler and E.K.U. Gross, Density-functional theory (Springer, Berlin, 1990); N.H. March, Electron density theory of atoms and molecules (Academic Press, New York, 1992). [6] X. Gonze, Phys. Rev. A 52 (1995) 1086.
1989);
M.K. Ha&la,
A. Banerjee/
Physics Letters A 222 (1996) 315-323
323
[ 71 G. Mahan and K.R. Subbaswamy, Local density theory of polarizability (Plenum, New York, 1990), and references therein. [8] X. Gonze and J.-P. Vigneron, Phys. Rev. B 39 (1989) 13120; X. Gonze, Phys. Rev. A 52 (1995) 1096. [9] P Hohenbetg and W. Kohn, Phys. Rev. 136 (1964) B864. [ 101 R.G. Parr, R.A. Donnelly, M. Levy and WE. Palke, J. Chem. Phys. 68 (1978) 3801. [ 111 1.M. Gelfand and S.V. Fomin, Calculus of variations (Prentice-Hall, Englewood Cliffs, NJ, 1963) [ 121 M.K. Harbola, Phys. Rev. A 48 ( 1993) 2696; Chem. Phys. Lett. 217 (1994) 461; lnt. J. Quantum Chem. 51 (1994) 201; Solid State Commun. 98 ( 1996) 629. [ 131 L.J. Bartolotti, lnd. Acad. Sci. (Chem. Sci.) 106 (1994) 103. [ 141 CA. Coulson, Proc. R. Sot. Edinburgh A 61 (1941) 20; G.S. Sewell, Proc. Cambridge Philos. Sot. 45 (1949) 678. [ 151 M.K. Harbola and A. Bane+, Phys. Rev. A 54 ( 1996) 283. [ 161 CF. von Weizsacker, 2. Phys. 96 (1935) 431. [ 171 E. Clementi and C. Roetti, At. Data Nut. Data Tables 14 ( 1974) 179. [ 18) M.N. Grasso, K.T. Chung and RI? Hurst, Phys. Rev. 167 (1968) 1. [ 191 J. Rice, P.R. Taylor, T.J. Lee and J. Almolf, J. Chem. Phys. 94 ( 1991) 4972.
Physics Letters A 222 ( 1996) 3 1S-323
ELSEVIER
Perturbation theory in terms of electron density Manoj K. Harbola, Amp Banerjee Laser Programme, Centre for Advanced Technology, Indore 452013, India
Received 25 June 1996; revised manuscript received 3 September 1996; accepted for publication 6 September 1996 Communicated by B. Fricke
Abstract When an inhomogeneous electron gas is subjected to a perturbation, its energy and density both change by small amounts. We calculate the changes in the energy explicitly in terms of the density changes within the density-functional theory of many-electron systems. We also derive the equations for the induced densities, and using these show that a density correct up to order n in terms of the perturbation parameter is sufficient to give energy changes up to order (2n + 1) . As a corollary, the even-order energy changes E(“‘) are variational with respect to the density changes p(“). The equations for the induced densities also follow from this corollary. The even-order corollary also gives a variational method of calculating the induced densities. The theory is demonstrated by applying it to calculate the polarizability and hyperpolarizability of the hydrogen, helium and neon atoms.
When a system of N electrons throughout the paper)
Ho
=&o: +uo(ri)l i=l
+
described
$c (izj)
by the Hamiltonian
(atomic
units,
lel = m, = Fi = 1, are used
-i ‘ij
is subjected to an applied perturbative field Au(l), the changes in its properties are also small and as such can be described by an expansion in powers of the perturbative parameter A. Thus the stationary state wavefunction !P and the ground-state energy E for the Hamiltonian HO + Au(‘) are written as [l] !P = To +
c
AW’)
(2)
i=l
and E = E. +
F A'E('),
(3)
i=l
where PO and EO are the wavefunction and the energy, respectively, corresponding to the Hamiltonian I$ (1). The corrections P(‘),P(2), . . . to the wavefunction and EC’), EC*), . . . to the energy are calculated 0375~9601/%/$12.00 Copyright 0 1996 Published by Elsevier Science B.V. All rights reserved. PII SO375-9601(96)00675-S
of by
316
M.K. Ha&la.
A. Banerjee/Physics
Lerters A 222 (1996) 315-323
solving the perturbation equations [ 11 to the appropriate order. Interestingly, a wavefunction correct only up to order n is required [ 21 to calculate the energy up to order (2n + 1) . This is the (2n + 1) theorem, The theorem is valid not only for the exact solutions of the Schriidinger equation but also [ 31 for approximate methods such as Hartree-Fock (see for example Ref. [ 41) theory. The purpose of this Letter is to develop perturbation theory in terms of the density within density-functional theory (DFT) [ 51 of many-electron systems. Thus we derive: (i) the energies to various orders in terms of the induced densities to appropriate orders, and (ii) the equations for the induced densities to each order, Using these equations, we show that the density correct to order n in terms of the perturbation parameter is sufficient to give energy to order (2n + 1) ( (2n + 1) theorem), as expected [ 61 from the variational nature of the energy functional with respect to the electronic density. As a consequence of this, the even-order corollary of the theorem follows [ 61: if the density is known exactly to order (n - 1), the even order energy change E(2n) is variational with respect to the density change p (‘). This is then shown to lead variationally to the equations for induced densities to each order. Further, this also gives a variational method of determining the induced densities. We demonstrate the theory developed with the examples of the hydrogen, helium and neon atoms. Before presenting our work, we point out that the theory developed below is distinct from the well established perturbative scheme [ 7,8] within the Kohn-Sham approach [ 51 of DFT. The Kohn-Sham method is an orbitalbased approach and as such its methodology [ 7,8] is similar to that of Hartree-Fock theory [ 31. On the other hand we do everything in terms of the density, obviating thereby the need to invoke single-particle orbitals. In DFT the ground-state energy of a system of electrons is expressed as a functional E[p] of the density p(r). This functional is made up of two components,
UpI = Npl +
s
uext(r>dr>drv
(4)
where F[ p] is a universal functional of p(r) and comprises the kinetic and interaction energy of the electrons; bxt is the external potential electrons are moving in. The energy is shown [9] to be variational with respect to p(r), and therefore the ground-state density is obtained by solving the following Euler equation,
SE[PI --p=o.
(5)
Q(r)
Once the ground-state density is known, the corresponding energy is given by the functional E[ p] . In Eq. (5) p, identified as the chemical potential [lo] (also see the discussion in Ref. [5] ), is the Lagrange multiplier to ensure that the total number of electrons is conserved. Now for the system with Hamiltonian Ho (Eq. ( 1) ), with ground state energy Eo, density pa(r) and the chemical potential ,UQ,Eqs. (4) and (5) can be written as
Eobol = Rpol +
J
uo(r>po(r>
dr,
(6)
and
SFbol G(r)
+ uo(r) = piI,
where SF [ po] /6p( r) means that the functional derivative SF/Sp is evaluated at p(r) will be followed throughout the paper. When the perturbation Au(‘) is applied, the density changes to po + Ap with
J
Ap( r) dr = 0,
(7) = po. The same notation
(8)
M.K. Harbola, A. Banerjee/Physics
and the chemical
Etpo+bl
potential
Letters A 222 (19%) 315-323
317
to m + A/L. From Eq. (4), we have
=F[po+Ap]
+
which can be Taylor expanded
(uo+~~(l’)(Po+Ap)dr,
(9)
J
to get
+S~opodr+/L’aapd’+~/~“)podr+~/~“)A~dr. Using Eqs. (6) -( 8) this can be rewritten as
E[pO+ Ap] = Eo[Po]+ g ; n=2
+ A Similarly
J
u”‘Podr+
A
.
1 Spo~~;r~~MISp(T . .Ahp(r,) n fph)Ap(r2). ...
dr, dr2..
. dr,
J
(11)
#)Apdr.
from Eq. (5), /-LO+ Ap can be written as
@‘[PO + API + f@(r) + /G(r) /LQLO+A/X= G(r)
=%‘+%
n=l
.
S("+')F[po] Adm> drl Wrl) ...Q(r,) J~p~r~ap~r~~
dr, +ua(r)
+ Au”)(r). (12)
Since Ap(r) Ap(r)
and Ap can also be expanded perturbatively,
we write
= xA’p(‘)(r),
(13)
and Ap
=
c
A’,&‘,
(14)
i
with
J#j(r) dr=O
(15)
for each i. In the following we first demonstrate with the examples of E(l), . . . , EC51 that the induced density only up to order two is required to calculate energy up to the fifth order. We then give the general proof for the (2n + 1) theorem. From Eqs. (11) and (12) expressions for EC’), . . . ,Ec5) and ~(~),p(~) are E”’ =
Ec2’ =
J J
u(‘)(rl)pO(rl)
dl)(rl)p(‘)(rl)
drl,
dr, + 12
(16)
J
~c2)mol
G(nhMr2)
p(1)(r,)p(1)(r2) dr, dr2,
(17)
M.K. Ha&la.
318
Ec3’=
u”)(r,)p’2)(r,)
s
S’3’F[pOl
+t
Q(r1
p(1)(r,)p(2)(r2) drt dr2
)G(r2)
p(‘)(r,)p”)(r2)p(‘)(r3)
dr, dr2,
(18)
Sd~l)~~(~2)&'(~3)
s
s
Letters A 222 (1996) 315-323
~(2’FbOl
dr, +
s
Ec4)=
A. Banerjee/Physics
up
dr, + L 2
~‘2’Fbol
+
~c2mkd s
dr, dr2
p’2)(q)p(2)(r2)
&J(Q)&d~2)
p(‘)(r1)p(~)(r2) drt dr2
&J(r1h%(r2)
s
~(3)Fbol
+; J I
+2;i Ec5’=
+(rl
p(‘)(r’)p(1)(r2)p(2)(r3)
s
6(4’F[Pol
p(‘)(r’)p(‘)(r2)p(l)(r3)p(‘)(r4)
6p(rl)Sp(r2)6p(r3)Sp(r4)
u”‘(r’)pc4)(r’)
SC2’F[pol
dr’ +
J
.I
pc1)(r’)pc4)(r2)
J J ~'3'Fbol +; J +;J J ~(3’Fbol
+;
drl drzdrgdr4,
(19)
dr’dr2
~p(r’)G(r2)
Sc2’Fbol Sp~rl~Gp~r2~pc2)(n )d3’(rd
+
dr’drzdq
)&(r2)+(r3)
drl dr2
p(‘)(r’)p(‘)(r2)p(3)(r3)
dr’ drzdr3
p(2)(r’)p’2’(r2)p(‘)(r3)
dr’dr;!dq
&(rl)Mr2Vp(r3)
Wrlb%(r2>Mr3)
~(4'F[pol
p(‘)(r’)p(1)(r2)p(1)(r3)p(2)(r4)
dr’ drzdrsdr4
Sp(rl)Sp(r2)6p(r3)Sp(r4)
+m
~‘5’F[pol
1
Sp(rl
)sp(rz)sp(r3>sp(r4)sp(r5)
XP (‘)(r’)p(1)(r2)p(1)(r3)p(1)(r4)p(1)(r5)
dr’drzdq
dr4dr5
(20)
and
p(I) = u(‘)(q) #I
+
J
a2F[pol p(‘)(r1) dn, Jjp(r)~p(rl>
a2F[pol
=
J
ap(r
pc2)(r’)
dr’ + i
J
(21) ~3F[~ol
pt1)(r’)p(‘)(r2)
dr’dr2.
(22)
G(r)MrI)&(r2)
Having been derived from the Euler equation (Q. (5) ) for the density, Eqs. (21) and (22) are the equations for the induced densities p(l) and pC2), respectively. Later in this paper we also derive them variationally from the even order corollary of the (2n + 1) theorem. Using Eqs. ( 15)) (21) and (22) expressions for E”) , . . . can be simplified (E(l) and Ec2) remain the same) and written as E’3’
= t
I
St3’F[ PO] ~p(n)~p(r2)@P(r3)
p(‘)(r’)p(‘)(r2)p(‘)(q)
drl dr2dr3,
(23)
M.K. Ha&la,
E(4)
2
E’5’
1
319
drt dr2
p(2)(r,)p(2)(r2) Q(u)Mr2)
s
J J
~(3’wcll
+(rl
fz
Letrers A 222 (1996) 315-323
6(2’F[po]
= I
+;
A. Banerjee/Physics
drr dr2drg
p’1’(rl)p”‘(r2)p’2’(~3)
)+(r2)&‘(r3)
~(4mPol
drr dr2dr3drq,
~~‘~(~1)~~‘~(~2)~~‘~(~3)p(l~(~4)
‘%(rl)+(r2)+(r3)&‘(r4)
J +i J S'4'F[Pol J
(24)
~c3m%l
= 1 2
S.(rl)Gp(r2)6p(r3)
d2)(rl
)d2b2)d1)(r3)
drl
dr2
dr3
drr drzdrsdrq
p’11(~,)p”‘(~2)pf1f(~3)pf2’(~4)
&drl)+(r2)@(r3)&‘(r4)
+m
@5’mol
1
+(n
)~~(r2)@(~3)~A~4)~~(~5)
(25)
Note that for the applied perturbation considered here (first order in A), energy corrections from Ec3) onwards are calculated from the universal functional F [ p] alone. It is evident from these expressions that p(t) alone gives an energy up to Ec3) and p(*) up to E (5) . This indicates the existence of the (2n + 1) theorem. We now prove it in general. By writing out the full expression for E (2n+1) from Eq. ( 1 I), and making use of the expression for p (n) from Eq. (12), we show that all the terms in the expression for Ec2”+*) involving induced densities of the order greater than n integrate to zero. From Eq. ( 11) the general expression for Ec2”+‘) is
,$2n+‘)
=
J X
u(1)(r,)p(2”)(rl) drt
J
Similarly
kl
@+‘Fbol
G(r1
=k 2 j=l
X
J
2 ili2...ik+l=l
pyq)p(iqr2)
. . . Mn+1)
)@(f2)
the expression
p(n)
+ 2
S(ir + i2 +
. . . +ik+l - (2n+ (k+
1))
l)!
. . . p(“+‘) (rk+l > drt dr;! . . . drk+, .
(26)
for p (n) (n > 1 )can be written from Eq. ( 12) and is given as
&it
+htl;.+ij-n)
ili2...ij=l
@+'Fbol aP(rl)aP(r2)
. . . aP(rj+l)
p(i’)(r2)p(i2)(T3)
. . . pciJ)(rj+l)
dr2dr3..
. drj+r.
(27)
As discussed above, Eq. (27) is the equation for the induced density p (“). In Eq. (26) the largest value of k for which the expression will have induced densities of order (n + I) and above is n; all the lower values of k would involve higher order induced densities. We therefore partition the k sum as follows,
320
M.K. Ha&da.
pnf’)
=
u(1)(r,)p(2”)(rl) dr,
J X
J
+ 2 5 k=l iliz...it+l=l
~k+lFbol
~P(~l
)Sp(r2)
Letters A 222 (19%) 315-323
6(i, + i2 + . . .+&+I
p(‘qr*)p(iqr2)
’ . .8P(Tk+, 8(i,
+F
(k+
1))
- (2n+
l)!
. . .pci’+‘)(rk+,)
dr, dr2..
. drk+,
. . .#i”+‘)(rk+,)
dr, dr2..
. drk+,.
)
+&+...+ik+l
-(2n+l))
2 k=(n+l)
X
A. Banerjee/Physics
J
(k+
i,iz...it+l=l
~k+lmol
&‘(rl
)@(r2)
. . . &'('.k+l
)
l)!
#&l) (r, )p(i2) (Q)
The third term in the above equation (28) involves pCn, and lower order induced densities. We now separate the highest order induced densities p(*“, and P(~‘+~, (m = 1, . . . , (n - 1) ) in the first two terms and write Eq. (28) as
+5
J
drlp(2”-“t)
[@*+” (m+‘) c c (k+
k=l
llt=l
X
(rl)
l)6(il +“;;;l;;“‘l))
iliz...ik=l
J~~~~~~~~~2~!po~,,,+,~ di”(r2)di2)(rd..
+ (terms involving
.#‘)(rk+d
ptn, and lower order induced densities).
dQ..
. drk+,]
(29)
Notice that once p(2”-“‘) is taken out, the maximum that k or il, i2,. . . can go up to is (m -t1) because of the delta function. The factor (k + 1) arises because p(2”-m) can be chosen in that many different ways. The expressions in the square brackets above, however, are the expressions for I_c(‘, (IQ. (21)) and P(~+,) (Eq. (27) ), respectively. Thus by the normalization condition (Eq. ( 15)) the integral s dr,p(2”-m) ( r,)p(m+‘) is zero. The highest order of the induced density in the terms left is therefore n. One can similarly show that the next order correction to the energy, i.e., E (2n+2), does require the induced density p(“+‘,. This is because for this energy correction one of the terms involving p (‘+l) is quadratic (the rest of the terms with p(“+*) are linear in it and have products of induced densities of order n or less as factors). This term comes from the k = 1 term in the sum and has the factor 1/2k! (= $) which is not the right factor ( l/k! = 1) for pCn+‘, to be pulled out of it. Thus these terms do not integrate to zero. This completes the proof of the (2n + 1) theorem. As mentioned in the beginning, the results above follow [6] from the variational nature of the energy functional. From Eq. (5) it is evident that any number conserving changes in the density would give an energy change which is quadratic in Ap. Thus for a trial density p’, there exists a number K such that [ 6,111
0 6 E[dl
- E[pl G Klb - ~11~7
(30)
where p is the exact density. Therefore, if the trial density is correct up to order n, the error in the energy is of the order (2n + 2). Thus for this density the energy is determined correctly up to order (2n + 1). From Eq. (30), it is also clear [ 61 that E(2n+2) is variational with respect to a trial p(“+‘) if the density is known exactly up to order n. It is also evident from the discussion above that the even order energy E(2n+2) has one term which is quadratic along with others which are linear in p(‘+‘). The variational nature of the even order energy terms leads to the Euler equations (27) for the induced densities p(“) with pCn) as the Lagrange multiplier to satisfy Eq. ( 15). Thus by minimizing Ef2) with respect to p(‘,, Eq. (21) for p(l) is obtained; ,u(‘, is the corresponding Lagrange multiplier. Knowledge of p(l) gives
M.K. Ha&la,
A. Banerjee / Physics Letters A 222 (19%) 315-323
321
energy correct up to order Ec3). Now p(l) is kept fixed and Eq. (22) for p’*’ is obtained by minimizing Ec4) with respect to p(*), with PC*) as Lagrange multiplier. In this manner one can go on to obtain higher order corrections both to the density as well as the energy. These equations are DFT counterparts of the more familiar equations 1l] for wavefunction corrections to various orders. As pointed out before, the theory above is different from the orbital-based perturbation scheme within the Kohn-Sham DFT methods. The difference between the two approaches is quite evident when the expressions for the energy changes to various orders in the two schemes are compared. As shown above, in the present method energy corrections from EC3) onwards are given by the universal functional F[ p] alone. In the KohnSham method, on the other hand, energy changes to all orders involve [ 7,8 ] the full Hamiltonian including the applied perturbation. In the past we have calculated [ 121 polarizabilities (energy up to order EC*)) of atoms and metal clusters variationally by employing the expression for EC’) given by Eq. ( 17). Bartolotti [ 131 has also mentioned the (2n + 1) theorem in terms of the density and has calculated polarizabilities of atoms using Eq. ( 17). However, it is only now that density based perturbation theory has been fully developed. This makes it possible to calculate higher order response functions of a system in terms of its density. As a demonstration we apply the theory presented above to calculate variationally the hyperpolarizabilities (energies up to the fourth order) for the hydrogen, helium and neon atoms subjected to a static electric field. Induced densities for the hydrogen atom in a static electric field are known [ 14,151 exactly, and are given as an
p(‘)
=
(2r + r*) cosBp0,
p ~*‘~[($~*+~~~+~~~)C0s~e+(~r~+~r3-~)1p0.
(31)
Taking a cue from Eq. (27), the variational form for the induced densities we choose is p”‘(r)
=Al(r)cosBpo(r),
pC2)(r) = [AZ(~) +43(r)
Al(r)=(aIr+b,r2+c,r3+dIr4+e,r5), 43(r)
=
(a3r
+
by* + c3r3 + d3r4 + e$),
cos2Blpo(r)
AZ(r) = (u2r +
b2r2
+
(32)
c2r3 + d2r4 + e2r5),
(33)
with al,. . . , ei, a*, . . . , e2 and ag, . . . , es being the variational parameters. The universal functional F[ p] for the hydrogen atom represents the kinetic energy of the atom, and is given exactly by the van Weizsacker functional [ 161 Tw[Pl = $
I
IVp(r) p(r)
I2 dr
(34)
For the helium atom we perform our calculation within the Hartree-Fock theory so that the universal functional F[ p] has two terms additional to the kinetic energy (given by Eq. (21) ) . These are the Hartree energy, EHbl
=
IS
$
(35)
and the exchange energy E, [p] , which for the helium-like atoms is given exactly as - $EH [ p] . The ground-state density of helium is obtained from its analytic Hat-tree-Fock wavefunction 1171. As described above, variational calculations to determine p(‘) and p’*) are performed in a sequential manner [ 151. First EC2’ is minimized to obtain p(l). Then keeping p(l) fixed p(*) is varied to minimize Ec4). This then is the density based counterpart of the similar wavefunctional [ 181 (also see Ref. [3]) or Kohn-Sham orbital method [8,15]. When this procedure is applied to the problem [ 141 of the hydrogen atom in a static field, it leads to the correct densities, as given by Rqs. (31), and the exact energies [ 141 EC*) = -$ a.u. and
M.K. Harbola, A. Banerjee/Physics
322
Letiers A 222 (1996) 315-323
E(4)
= -F a.u. for it (E(l) = Ec3’ = Ec5) = 0 by symmetry). Notice that Ec4) is calculated from the kinetic energy functional alone. Thus all the information about the applied perturbation is carried to EC41 by the induced densities. This, as discussed above, is in contrast to the Kohn-Sham method [ 8,151 where the expectation value with respect to the induced orbitals involves the applied perturbation also. For the helium atom the second order energy is -0.6611 a.u. This gives the polarizability to be 1.322 a.u. which is equal to the Hartree-Fock theory result [ 15,191. The fourth order energy is -1.5008 a.u. which gives the hyperpolarizability to be 36.02. This is close to but slightly lower (because of the approximate form chosen for the induced densities) than the more accurate Hartree-Fock value of 37.77 a.u. [ 151. In the calculations above, the kinetic energy and the exchange energy are treated exactly. Thus it is not clear how accurate a density-based calculation would be if these energies are to be approximated. For this reason, we have also performed the calculations within exchange-only for the neon atom with the kinetic energy given by Eq. (34) and an exchange energy by the local density approximation [5] as 113 E,LDA[P]
= -t
a 0
p4’3 (r)
dr.
(36)
s
The reasons for choosing the von Weizsacker functional over the other kinetic energy functionals such as the Thomas-Fermi [ 51 or the gradient-expansion approximation [ 51 functionals are discussed elsewhere [ 121 in the context of polarizability calculations. More discussion will follow in our subsequent work. When HartreeFock densities [ 171 are employed to calculate the response properties of neon, they give cy = 3.1 a.u. and y = 74.15 a.u. which compare well with the coupled Hartree-Fock values [ 191 of 2.34 a.u and 71.9 a.u., respectively. This, along with the calculation for hydrogen and helium, demonstrates the applicability of the density-based perturbation theory, the (2n + 1) theorem and its even-order corollary to orders higher than two. These calculations also show that although in principle exact knowledge of induced densities to order (n - 1) is required to determine p(“) variationally, in practice a sequential variational procedure also leads to stable and accurate results. For a further discussion in terms of Kohn-Sham orbitals, we refer the reader to Ref. [ 151. To conclude, we have developed a density based perturbation theory and shown that in the perturbative regime, ground-state density accurate to order n is sufficient to determine energy up to order (2n + 1) . Further, the even-order energy changes EC*“) are variational with respect to the nth order density change p(“). This then gives the method of calculating the induced density, and leads to equations that determine them. These results would make perturbative calculations within DFT quite easy as compared to orbital based methods, as demonstrated in Refs. [ 12,131 for second-order calculations. We have demonstrated the applicability of the theory by employing it to calculate linear and nonlinear response properties of hydrogen, helium and neon atoms. Work in the direction of application of the theory to other atoms is in progress and will be reported in the future. We thank Dr. Kailash Rustagi for his constructive
comments.
References [ I ] L.I. Schiff, Quantum mechanics (McGraw-Hill, New York, 1969). [2] A. Dalgarno and A.L. Stewart, Proc. R. Sot. A 238 (1956) 269; PO. Lcwdin, J. Molec. Spectr. 13 ( 1964) 326. [3] PW. Langhoff, J.D. Lyons and RI? Hurst, Phys. Rev. 148 (1966) 18. [4] CF. Fischer, The Hartree-Fock method for atoms (Wiley, New York, 1977). [5] R.G. Parr and W. Yang, Density-functional theory of atoms and molecules (Oxford Univ. Press, Oxford, R.M. Dreizler and E.K.U. Gross, Density-functional theory (Springer, Berlin, 1990); N.H. March, Electron density theory of atoms and molecules (Academic Press, New York, 1992). [6] X. Gonze, Phys. Rev. A 52 (1995) 1086.
1989);
M.K. Ha&la,
A. Banerjee/
Physics Letters A 222 (1996) 315-323
323
[ 71 G. Mahan and K.R. Subbaswamy, Local density theory of polarizability (Plenum, New York, 1990), and references therein. [8] X. Gonze and J.-P. Vigneron, Phys. Rev. B 39 (1989) 13120; X. Gonze, Phys. Rev. A 52 (1995) 1096. [9] P Hohenbetg and W. Kohn, Phys. Rev. 136 (1964) B864. [ 101 R.G. Parr, R.A. Donnelly, M. Levy and WE. Palke, J. Chem. Phys. 68 (1978) 3801. [ 111 1.M. Gelfand and S.V. Fomin, Calculus of variations (Prentice-Hall, Englewood Cliffs, NJ, 1963) [ 121 M.K. Harbola, Phys. Rev. A 48 ( 1993) 2696; Chem. Phys. Lett. 217 (1994) 461; lnt. J. Quantum Chem. 51 (1994) 201; Solid State Commun. 98 ( 1996) 629. [ 131 L.J. Bartolotti, lnd. Acad. Sci. (Chem. Sci.) 106 (1994) 103. [ 141 CA. Coulson, Proc. R. Sot. Edinburgh A 61 (1941) 20; G.S. Sewell, Proc. Cambridge Philos. Sot. 45 (1949) 678. [ 151 M.K. Harbola and A. Bane+, Phys. Rev. A 54 ( 1996) 283. [ 161 CF. von Weizsacker, 2. Phys. 96 (1935) 431. [ 171 E. Clementi and C. Roetti, At. Data Nut. Data Tables 14 ( 1974) 179. [ 18) M.N. Grasso, K.T. Chung and RI? Hurst, Phys. Rev. 167 (1968) 1. [ 191 J. Rice, P.R. Taylor, T.J. Lee and J. Almolf, J. Chem. Phys. 94 ( 1991) 4972.