The calculation of magnetisabilities using current density functional theory

28 October

1994

CHEMICAL PHYSICS LETTERS ELSEVIER

Chemical

Physics Letters 229 ( 1994) 225-232

The calculation of magnetisabilities using current density functional theory Aaron M. Lee, Susan M. Colwell, Nicholas C. Handy University Chemical Laboratory, Lensjeld Road, Cambridge CB2 IEW UK Received 3 August

1994; in final form 24 August

1994

Abstract

The theory for the calculation of magnetisabilities using current density functional theory, which follows from the original theory of Vignale, Rasolt and Geldart, has been implemented. We present an initial application of this theory to the set of small molecules HZ, HF, Nz. CO, H20, and NH3.

1. Introduction

The application of density functional methods to molecular properties is of interest because experience with calculations of geometries and vibrational frequencies indicates that density functional theory (DFT) may have the potential to provide answers that are more accurate than those given by Hartree-Fock self consistent field (SCF) calculations, but at similar cost. SCF calculations of electric and magnetic properties are known to be inaccurate, but in general correlated methods, particularly for frequency-dependent properties, are complicated and expensive. There have been many calculations of electric properties using finite field DFT techniques, and recently some using coupled Kohn-Sham theory [ 1,2]. Although there have been calculations of magnetic properties using DFI (e.g. Malkin et al. [ 31) these have not been based on a full coupled Kohn-Sham treatment. The difficulty, which is widely acknowledged, is that with standard DFT using an exchange-correlation functional which depends only on the density, the corresponding exchange-correlation part of the mag-

0009.2614/94/$07.00 @ 1994 Elsevier Science B.V. All rights reserved rcnrn,?nn .lI, “ln”\,T.,x,T_ _

netic Hessian vanishes. Therefore it is necessary to use a density functional theory which introduces the current density and specifically introduce an exchange-correlation functional which depends on both densities. The basic theory for this current density functional theory (CDFT) has been presented by Vignale, Rasolt and Geldart [4], and they have given the specific form of the Kohn-Sham equations. In a previous paper [5] we presented an explicit expression for the magnetic Hessian within this theory using the local density approximation. We found that a Hellmann-Feynman identity holds, and hence the expressions for second-order properties are straightforward. In this Letter, we discuss the computational implementation of these expressions and present some initial results on a sequence of small molecules. The total magnetisability is not directly measurable, although the paramagnetic part can be obtained from measurements of the rotational g factor and therefore we compare our results with previous calculations rather than with experiment; in particular we refer to a recent series of papers by Cybulski and Bishop [ 69] which use analytic derivatives of the second-order

226

A.M.

Lee et al. /Chemical

Physics

Moller-Plesset (MP2) energy. Other correlated calculations of magnetisabilities have used multiconfiguration self-consistent field (MCSCF) methods (in particular CASSCF) [lo] and configuration interaction methods [ 111 or the second-order polarisation propagator (SOPPA) method [ 121. One problem with all methods is that the answers are gauge-origin dependent with finite basis sets. Various schemes have been used to overcome this problem for SCF based methods, namely LORG (localised orbital/localised origin) [ 131, GIAO (gauge invariant atomic orbitals) [ 141 and IGLO (individual gauge for localised orbitals) [ 151. These methods have been combined with correlated methods to give, for example, MC-IGLO [ 161. All of the correlated calculations are expensive, and the different methods do not always agree even as to the sign of the correlation corrections. This suggests that the perturbation theory methods may not have converged and should ideally be taken to higher order, further complicating these techniques. Density functional theory may provide a much simpler way to include correlation effects.

Letters

229

(1994)

where o(r), J[p], and Ts[p] have the usual definitions and E,, [p, v] is the exchange-correlation energy. The density and current density are now given by

and

j,(r)

(4) The generalised Hohenberg-Kohn theorem [ 41 states that 8”~ is minimised when p and j, take the correct values for the potentials A and Vext, and this minimum value is the true ground state energy for the system. Writing the equation for the non-interacting orbitals, +p, in the form

&(-ifiV 2. Current

density functional

theory (CDFT)

= The theory for this section was presented in full in Ref. [ 51, but we outline it here for clarity. Within CDFT [ 41 the exchange-correlation energy is a function, not only of the density p(r), but also of the paramagnetic current density j,(r) It is convenient to express the problem in terms of p and of the gauge invariant vector v(r) = V x [j,(r) /p( r) 1. Therefore the CDFT energy of a system of N (= 2M) electrons moving in an external electromagnetic field (E, B), i.e. under external potentials 4, A is

&"A[P,Vl + e

Jp(r) (

=&oh~l+ j,(r).Adr,

(1)

s where E = -04

+ eA,E)‘+

kdr)

q$

cp+p

(3

9

one finds expressions

for V& and A,E in terms of

K,(r) =

and eA,,(r)

= LO

x

P

Hence, given Ex,[ p, v] Eq. (5) can, in principle, be solved to give the orbitals, and hence the density, current and energy. In standard Kohn-Sham [ 171 language

and B = V x A. Vext = -e4,

EoIp,vl = K:,[p, jpl +

+

225-232

JLPI + &[P,

J

VI 3

p(r)

and

E = 2(ilh,ffli) -e

dr

(2)

J

+ 2(iiljj)

j,(r)-&,dr,

where we define

+ E,,[p,vl

(8)

A.M. Lee et al. /Chemical

h,n = ho - g

(A-V

+ V-A)

+ $A2

+ v,,,

+ - E

(A,, . V + V-A,,)

= h + h,,.

(9)

3. Coupled-perturbed (CPCKS)

227

Physics Letters 229 (1994) 225-232

s

gtdlW12dc

we find expressions sians [5]

for the electric and magnetic

Hes-

- ei)Soi,b,i + 4(ailbj)

Htni,bj = (8

current K-S theory

(15)

(16)

In Ref. [ 51 we showed that t

10)

We have used the standard repeated suffix notation, with i,j, . . . being occupied orbitals, a, b, . . . being virtual orbitals and p, q, . . . being any molecular orbitals. Note that in Eq. (10) it is dh/d3* not dh,n/d3* i.e. a Hellmann-Feynman identity holds for CDFT. Therefore first- and second-order electric and magnetic properties are given by expressions which are formally the same as in the SCF case. The Kohn-Sham matrix is FE

= tplh&)

+2(&i)

= ens,,.

+ (plv,&)

-4e/2g(Vx$).(Vx%)

The derivative of a molecular orbital with respect to the applied field, 3*, can be written as

’ = -g[4iV4q Jpq

- V+T,$q]

= jf.

(18)

HI and H2 are real, and so the same deductions can be made as in the RHF case. For purely electric fields U is real and one recovers the usual CPKS equations. For purely magnetic fields U is pure imaginary, hence the magnetisability is given by + 4Im( U,i)‘N~,

(19)

where iNL = (ilf$lB=o/a) A 2+$‘(r x V&),dr,

= -g

(12)

(20)

J’ and

where we define ur:, = o;Jl IF=O’

(13)

Differentiating the condition Fzs = 0 with respect to 3A, and evaluating at 3 = 0 gives the coupledperturbed current Kohn-Sham equations for the U matrices,

(21)

4. Computational

+ i Im t Ubj) AH2ai,b,j = 0.

(17)

where

X,+ = -2QF (11)

dr,

(14)

Using the expression given by Vignale, Rasolt and Geldart [ 41 for the exchange-correlation energy in the weak field case,

implementation

The DFT calculations require a suitable form for the exchange-correlation energy functional in a uniform applied magnetic field, E!$* [p, Within the local density approximation, the uniform electron gas has represented a starting point for such functionals. From Eq. ( 15) we can write ,?$F* [ p, v] as the sum of

v]

228

A.M. Lee et al./Chemical

Physics Letters 229 (1994) 22S-232

two components: EiF* [ p]~, which is the usual KohnSham exchange-co~elation functional in the absence of an applied magnetic field, and a v-dependent part, which according to Vignale, Rasolt and Geldart can be written as Jg(p)lv(r)12dr (the form for g(p) is given below). We will write this last functional as VRG. For E,, [ p] 0 we use the Dirac electron gas exchange formula [ 191, with correlation included explicitly by adding the Vosko-Wilk-Nusair [ 201 correlation functional (VWN) which is a fit to the correlation energy of the electron gas. Vignale, Rasolt and Geldart [25] give the following form for g(p) :

(22) where XL is the Landau orbital magnetic susceptibility, ,& = - ( e2kF/ 12&72) the Landau non-interacting orbital susceptibility, and kF(r) = [ 3~~2p(r)] ‘i3 the Fermi momentum. The ratio XL/X: is tabulated (Table 1 of Vignale, Rasolt and Geldart [25] ) as a functionofr,(r) = [47rp(r)/3]-1/3,0nr, E [O,lO].We fitted functions of the form exp( -aor,) x (al +a2rs + ug-f + . . .) to the tabulated data. In this work we use a0 = 0.042, at = I, a - 0.028, and a, = 0 for n > 3. The rms error m our t%was N 1.5 x 10P3, and it was found that variations in the fitted form used altered the components of the paramagnetic magnetisability tensor by an amount of this magnitude. This form is useful as it is everywhere positive ( rs > 0)) continuous, and has the correct value at rs = 0 (p --t co). We note that Vignale, Rasolt and Geldart [25] give an exact form for this function in the limit r, + 0, but that this form is unlikely to be useful for non-metallic densities. For the molecules examined here rs will usually be larger than 0.1. We wili denote the Dirac + VWN + VRG combination of functionals as LDA, and the Dirac + VRG combination as LDAx. In the calculation of molecular properties the Dirac exchange term alone is known not to give reliable results. Nonetheless, we have used both the LDA and LDAx functionals in the present work to illustrate the effect of including the VWN correlation term. For the numerical integrations needed for the density functional calculations, we used Euler-Maclaurin radial quadrature and Lebedev angular quadrature [ 181, and proceeded as detailed in previous work

[ I]. The size of the systems examined made precise quadrature possible; all figures calculated were stable to at least 10e4 atomic units. We confirmed the implementation of our formulae by comparing the total susceptibility tensor to the numerical second derivative of the energy calculated as a finite difference of energies at different applied magnetic field strengths. The CDFT energy for a closedshell system in the presence of an external electromagnetic field has been given, Eq. ( 1). The solution of the self-consistent equations now involves the construction and manipulation of complex Hermitian matrices, but the procedure is nonetheless straightforward. We represented the complex eigenfunctions in the presence of a non-zero magnetic field, as an expansion with complex coefficients Cai in a real atomic orbital basis of h?, functions,

(23)

@i(r) = &*is*(r). a=1 By defining the complex density matrix D as

(24) i=l

the density was expressed in terms of the real part DR ofthedensitymatrix,p(r) = C>i+, D&~~(r)qp(r), and the paramagnetic current density in terms of the

imaginary part D’, j,(r)

= (@‘tfz) ~~Pzl D$w(r)

xVqa( t). The energy and the Kahn-Sham matrix in the presence of a field were then written in terms of the atomic orbital basis. Agreement to 4 and 5 decimal places was achieved between our finite field and analytic calculations.

5. Calculations In Tables l-3, the results of our calculations of the magnetisabilities of the hydrogen, hydrogen fluoride, nitrogen, carbon monoxide, water and ammonia molecules are presented. Results are given in atomic units ( e2a$‘m) _For each molecule we give the nonzero components of the paramagnetic and diamagnetic magnetisability tensors, xr and xd. For the diatomics we also give the isotropic magnetisability, X = (2x1 + ,yll>/3, and the anisotropic magnetis-

A.M. Lee et al./Chemical Table I The magnetisability

of the (a) hydrogen

Quantity

(a)

xi (cm.)

-0.789

x$ (c.m.)

-0.929

Ax(c.m.)

0.021

Ax(H)

xi (cm) XT (c.m.)

x: (c.m.)

F(c.m.) Ax(c.m.)

Ax(F)

-0.769

-0.841

-0.806

-0.777

-0.767

-0.762

-0.909

-0.978

-0.942

-0.899

-0.894

-0.893

0.026

0.025 -0.915

-0.843

0.110

-0.915

0.024

-0.844

0.110

-2.114

-2.275

-2.221

-2.115

-2.193

-2.374

-2.541

-2.488

-2.373

-2.454

0.163 -2.235

0.113

0.113 -2.237

0.116

0.150

0.111

0.110

0.129 -2.362

0.111

-2.189

-2.189

0.110

0.130

-2.268

0.110

-2.295

0.104

0.148

0.148

-2.295

-2.361

-2.188

0.116

0.156

0.137

-2.188

-0.831

0.104

-2.436

0.113 -2.269 0.113

0.110

and (b) carbon monoxide

Basis A

molecules

in atomic units

Basis B LDA

SCF

Literature LDAx

LDA

SCF

a MI?

x;t(c.m.)

-3.870

-3.800

-3.861

-3.939

-3.850

-3.859

-3.809

xd,(c.m.)

-7.953

-7.874

-7.859

-7.995

-7.901

-7.857

-7.873

xP,(c.m.)

5.814

5.820

5.887

5.971

5.998

5.912

5.749

X( cm.)

-2.716

-2.634

-2.602

-2,782

-2.552

-2.583

-2.686

Ax(c.m.)

-1.731

- 1.746

-1.888

-1.914

- 1.948

-1.914

- 1.685

X(N)

-3.000

-2.915

-2.633

-2.684

-2.577

-2.583

-2.694

Ax(N)

- 1.304

-1.327

-1.841

-1.881

-1.910

-1.910

- 1.667

xi;CC) xd,(C) xp,(C)

-3.905

-3.825

-3.799

-3.963

-3.865

-3.796

-3.838

-13.321

-13.245

-13.332

-13.358

-13.267

-13.333

-13.233

10.331

X(C)

-3.295

10.378

11.207

-3.187

-2.683

11.142

11.189

11.311

1

0.104

-0.837

0.100

-2.160

0.164

-0.83

0.103

-2.472

-2.269

0.027

-0.836

0.100

-0.880

0.111

exact h

MF2

0.022

-0.880

0.111

0.190

LDAx

a Ref. 191.

SCF

-2.194

of the (a) nitrogen

Quantity

(b)

LDA

-0.895

0.194

a

h Ref. 1231.

Table 2 The magnetisability

(a)

in atomic units

LDAx

0.118

-2.271

X(F)

SCF

-0.848

0.119

229

Literature

0.022

-0.869

-0.919

X(H)

fluoride molecules

Basis B LDA

xy(cm.)

a Ref. 191.

Basis A LDAx

z(c.m.)

(b)

and (b) hydrogen

Physics Letters 229 (1994) 225-232

10.991

-2.799

-2.674

-2.613

-2.774

Ax(C)

-0.915

-0.957

- 1.674

- 1.746

-1.787

- 1.774

-1.596

F(c.m_)

-2.805

-2.662

-2.624

-2.743

-2.625

-2.607

-2.754

Ax(c.m.)

- 1.650

- 1.744

-1.762

- 1X30

-1.860

-1.784

- 1.626

230

A.M. Lee et d/Chemical

Table 3 The magnetisability

(a)

of (a) water and (b) ammonia

x:: (0)

Physics Letters 229 (1994) 225.232

in ntamic units

LDAx

LDA

LDAx

LDA

SCF

ME2

-3.324

-3.270

-3.419

-3.339

-3.215

-3.301

x~,fo~

-3.192

-3.137

-3.309

-3.224

-3.089

-3.187

X~)J(~~

-3.605

-3S52

-3.694

-3.619

-3.4%

-3.577

xf, (0)

0.265

0.266

0.214

u.270

0.270

0.260

x!&(O)

0.189

0.188

0.204

0.188

0.192

0.195

xyy(O)

0.539

0.538

0.570

0.546

0.546

0.544

X(O)

-3.0443

-2,989

-3.145

-3.059

-2.930

-3.022

zic.m.1

-3.035

-2‘985

-3.145

-3.058

-2.930

-3.02

ability Bx = 211 - xl; for the remaining molecules we give the isotropic mag~~t~sabiIity~ defined as X = $x~~+x~~+,Y~~) (wherexnl_L=x~,+x~,).The gauge for the calculation is given in parentheses following the quantity being calculated. Our initial calculations were performed with Dunning’s [ 21) DZP and TZZP bases. We present some TZ2P (basis A) results in the tables. For our larger calculations we chose to compare our results with those given in the papers I [S] and III [S] of Cybulski and Bishop’s recent series and we have used similar or identical basis sets to theirs (although our code uses Cartesian and not spherical Gaussians). Cybulski and Bishop construct their basis sets from the (10s) for H, and (13~8~) for C, N, 0, F uncontracted sets of van Duijneveldt [ 221. The large basis sets used for Hz, HF, H20 and NH3, are termed basis B. For Hz this is the uncontracted (1 ls9p6d2f) basis for If from paper III, and for the other molecules, the ( lOs6p id) / [ MpIsd] contraction for H, and ( 16s9p5d2f) / [ 12sQ4d2fj contractions for C, N, 0, and F of Cybulski and Bishop from paper IX [7]. Basis C was constructed from basis B with the f functions and most diffuse d function removed, i.e.

I

( 16s9p4d) /I 12sSp3d], and was used for Nz and CO. We compare our results to those of Cybulski and Bishop using the same or similar basis. For Hz, and HzO, and NH3 we reproduced exactly the SCF results of Cyhulski and Bishop (as given in the appropriate paper). For HF, CO and Nz, Cybulski and Bishop [9} report larger basis calculations, and therefore we give our SCF values as some form of comp~ison between both, We remark that large basis sets are certainly necessary for the accurate and gauge invariant calculation of magnetic properties. All our calculations were performed at the geometries used in Ref. [ 71.

Far hydrogen we used the bond length THH = 1.4 no, and for hydrogen fluoride the experimental bond length of ?&F = 1.7328 a~_ In Table 1 we see that the DFT calculations vary sign~~ca~tly between the two bases. The calculated values for X are much too negative with the exchange-only functional LDAx (x(c.m.) = -0.915 for Hz, and y(c.m.) = -2.361 for I-IF), with LDA improving on this (xfc.m.) = -0.880 for HZ, and x(c.m.) = -2.295 for HF). For

A.M. Lee et al./Chemical

Physics Letters 229 (1994) 225-232

Hz, with the LDA functional, the components of x in the centre of mass (c.m.) gauge using the smaller basis agree more closely than those of the larger basis with the exact values of Rychlewski and Raynes [ 231, which is fortuitous. We give x and Ax at both the cm. and H or F gauges showing that, for the larger basis, the results are less gauge dependent. 5.2. The nitrogen and carbon monoxide molecules For nitrogen and carbon monoxide, we used the experimental bond lengths of rNN = 2.07432 ao, and loo = 2.132 aa, respectively. For Nz, we obtained a value of X(c.m.) = -2.552 for LDA and X(c.m.) = -2.782 with LDAx. Our SCF value is z(c.m.) = -2.602 with this basis, whereas Cybulski and Bishop calculate z = -2.583 (SCF) and X = -2.685 (MP2) [ 91. For CO, we obtained a value of x(C) = -2.674 with LDA and x(C) = -2.799 for LDAx. Our SCF value is x(C) = -2.683 with the larger basis, while Cybulski and Bishop calculate x(C) = -2.607 (SCF) and x(C) = -2.754 (MP2). We expect our DFT results to change further with basis set, perhaps considerably. For both molecules, MP2 [9] and MC-IGLO give a diamagnetic (i.e. negative) correlation correction (7 = -2.578 (SCF-IGLO) and X = -2.651 (MCIGLO) for N2, and 7 = -2.611 (SCF-IGLO) and X = -2.725 (MC-IGLO) [ 161 for CO), whereas SOPPA gives a paramagnetic (i.e. positive) correction (X = -2.597 (RPA), 7 = -2.505 (SOPPA) for Nz, and 7 = -2.608 (RPA), X = -2.584 (SOPPA) [24] for CO). Our LDA results which explicitly include correlation via the VWN term, are more paramagnetic than our LDAx ones. 5.3. The water and ammonia molecules For water we used rou = 1.808846 ao, LHOH = 104.52”, and for ammonia, NH3 , rNH = 1.91316 ao, LHNH = 106.67”. The results for both molecules using the larger bases are encouraging; the LDA results being closer to the MP2 than SCF values. For H20, we find X(O) = -3.059 for LDA, compared with the values of Cybulski and Bishop [ 91 of X( 0) = -2.930for SCF, and X(O) = -3.022 for MP2. The LDA value in the c.m. gauge was -3.058 showing little variation with gauge. For NHs, X(N) = -3.781 for LDA com-

231

pared with the values of Cybulski and Bishop [9] of F(N) = -3.647 for SCF, and x(N) = -3.737 for MP2. Again little variation with gauge is shown.

6. Conclusions We have demonstrated that density functional calculations of magnetisabilities are feasible. It is clear that the density functional numbers are not yet in competition with those of the conventional correlated methods. It is also clear that reliable results for DFT will only be obtained with large basis sets, and that misleading conclusions may well be drawn from small basis set results. We observe the SCF and DFT results to be broadly similar, although the latter appear more sensitive to basis set effects. These are essentially preliminary calculations and we intend to investigate systematically the variation of results with basis set and ways of overcoming gauge dependence. As the structure of the orbital equations is the same as for Hartree-Fock, it should be possible to adapt the techniques (LORG [ 131, GIAO [ 141 and IGLO [ 151) developed in that field. The form of g( p) given explicitly by Vignale, Rasolt and Geldart [4] is appropriate only in the high density limit. We will investigate forms for EkF* [ p, v] appropriate for nonmetallic densities. The theory as presented here is only appropriate for local functionals; we will investigate its application to non-local functionals, i.e. those that depend also on the derivatives of the density. Nuclear shielding constants and hence chemical shifts are also determined by expressions formally similar to the Hartree-Fock case, consequently our program will calculate these as well. We expect the variation in these quantities to be larger than for the magnetisabilities, and so provide a more sensitive test of basis sets and functionals. An important aspect of Kohn-Sham density functional calculations of magnetisabilities is that twoelectron integrals are required only in the calculation of the energy and the Kohn-Sham matrix. From Eq. ( 17) it is seen that the DFT magnetic Hessian involves only one-electron integrals, in contrast with HartreeFock theory and methods based on it. The coupling term in the DFT magnetic Hessian is straightforward to calculate, and hence CPCKS calculations, although more expensive than uncoupled approximations (see

232

A.M. Lee et al./Chemical

Physics Letters 229 (1994) 225-232

Ref. 31)) are less costly than expected. By using direct methods whereby the required two-electron integrals are evaluated as needed in the construction of the orbitals and eigenvalues, the DFT calculation of magnetisabilities (and related properties such as chemical shifts) is feasible for systems of considerable size. Our aim remains to calculate frequency-dependent properties, and we are now in a position to do so.

Acknowledgement The authors wish to acknowledge useful discussions with Dr. Roger Amos. AML gratefully acknowledges the award of a Commonwealth Scholarship from the Commonwealth Scholarship Commission of the United Kingdom.

References [ 11 S.M.Colwell,

C.W. Murray, NC. Handy and R.D. Amos, Chem. Phys. Letters 210 (1993) 261; A.M. Lee and SM. Colwell, J. Chem. Phys., submitted for publication. [2] G. Senatore and K.R. Subbaswamy, Phys. Rev. A 34 (1986) 3619. [3] V.G. Malkin, O.L. Malkina and D.R. Salahub, Chem. Phys. Letters 204 (1993) 80, 87; V.G. Malkin, O.L. Malkina, M.E. Casida and D.R. Salahub, J. Am. Chem. Sot. 116 (1994) 5898. [4] G. Vignale, M. Rasolt and D.J.W. Geldart, Advan. Quantum Chem. 21 (1990) 235; G. Vignale and M. Rasolt, Phys. Rev. Letters 59 (1987) 2360; Phys. Rev. B 37 (1988) 10685.

[5] SM. Colwell and N.C. Handy, Chem. Phys. Letters 217 (1994) 271. [61 SM. Cybulski and D.M. Bishop, Mol. Phys. 76 (1992) 1289. 171 SM. Cybulski and D.M. Bishop, J. Chem. Phys. 98 ( 1993) 8057. [81 D.M. Bishop and S.M. Cybulski, Mol. Phys. 80 ( 1993) 199. [91 S.M. Cybulski and D.M. Bishop, J. Chem. Phys. 100 (1994) 2019. 1101 G.T. Dabom and N.C. Handy, Mol. Phys. 49 (1983) 1277. 1111 G.T. Dabom and N.C. Handy, Chem. Phys. Letters 81 ( 1981) 201. [I21 S.P.A. Sauer, V. Spikko and J. Oddershede, Chem. Phys. 153 (1991) 189. [I31 Aa. E. Hansen and T.D. Bouman, J. Chem. Phys. 82 (1985) 5035. [I41 K. Ruud, T. Helgaker, K.L. Bak, P Jorgensen and H.J.Aa. Jensen, J. Chem. Phys. 99 (1993) 3847. [I51 M. Schindler and W. Kutzelnigg, J. Chem Phys. 76 (1982) 1919. 1161 Ch. van Wiillen and W. Kutzelnigg, Chem. Phys. Letters 205 (1993) 563. [I71 W. Kohn and L.J. Sham, Phys. Rev. 140 (1965) Al 133. 1181 V.I. Lebedev, Sibirsk. Mat. Zh. 18 ( 1975) 132; Zh. Vychisl. Mat. Mat. Fiz. 15 (1975) 48: 16 (1976) 293; Russian Acad. Sci. Dokl. Math. 45 (1992) 593. [I91 P.A.M. Dirac, Cambridge Phil. Sot. 26 (1930) 376. 1201 S.J. Vosko, L. Wilk and M. Nusair, Can. J. Phys. 58 (1980) 1200. [21] T.H. Dunning, J. Chem. Phys. 53 (1970) 2823. [22] F.B. van Duijneveldt, IBM Res. Rept. ( 1971) 945. 1231 J. Rychlewski and W.T. Raynes, Mol. Phys. 41 (1980) 843. 1241 SPA. Sauer, 1. Paidarova and J. Oddershede, Mol. Phys. 81 (1994) 87. [25] G. Vignale, M. Rasolt and D.J.W. Geldart, Phys. Rev. B 37 (1988) 2502.