Levels of self-consistency in the GW approximation
aa r X i v : . [ c ond - m a t . o t h e r] J un l. Levels of self-consistency in the GW approximation ∗ Adrian Stan,
1, 2, 3
Nils Erik Dahlen, and Robert van Leeuwen
1, 2 Department of Physics, Nanoscience Center, FIN 40014, University of Jyv¨askyl¨a, Jyv¨askyl¨a, Finland European Theoretical Spectroscopy Facility (ETSF) Rijksuniversiteit Groningen, Zernike Institute for Advanced Materials,Nijenborgh 4, 9747AG Groningen, The Netherlands. (Dated: October 30, 2018)We perform GW calculations on atoms and diatomic molecules at different levels of self-consistency and investigate the effects of self-consistency on total energies, ionization potentialsand on particle number conservation. We further propose a partially self-consistent GW scheme inwhich we keep the correlation part of the self-energy fixed within the self-consistency cycle. Thisapproximation is compared to the fully self-consistent GW results and to the GW and the G W approximations. Total energies, ionization potentials and two-electron removal energies obtainedwith our partially self-consistent GW approximation are in excellent agreement with fully self-consistent GW results while requiring only a fraction of the computational effort. We also find thatself-consistent and partially self-consistent schemes provide ionization energies of similar quality asthe G W values but yield better total energies and energy differences. PACS numbers: 31.15.-p,31.15.xm,71.15.-m
I. INTRODUCTION
Green function methods [1, 2] have been very suc-cesful in the description of various properties of many-electron systems, ranging from atoms and molecules tosolids [3, 4]. Within the Green function approach, theseproperties are completely determined by the self-energyoperator Σ, which incorporates all the effects of exchangeand correlation in a many-particle system [1]. One of themost widely used approximations to the self-energy is the GW approximation ( GW A) [5]. In the GW A, the self-energy operator has the simple form Σ = − GW , where G is the Green function that describes the propagationof particles and holes in the system, and W is the dy-namically screened interaction. This quantity describeshow the bare interaction v between electrons is modifieddue to the presence of the other electrons and appearsas a renormalized interaction in terms of Feynman dia-grams. In extended systems the screened interaction ismuch weaker than the bare interaction, and therefore itis much more natural to expand the self-energy in termsof the screened interaction than in terms of the bare in-teraction. The lowest order in this expansion [5] is the GW A.Calculations within the GW A are usually done in twosteps. First, a density functional theory (DFT) [6] calcu-lation is performed and the DFT orbitals and eigenvaluesare used to construct a first guess G , for the Green func-tion and a first guess W , for the screened interaction.In a second step, the self-energy Σ = − G W is con-structed and the Dyson equation is solved for the Green ∗ See J. Chem. Phys. for the published version. function. In principle, this new Green function shouldbe used to calculate a new self-energy and this processshould be iterated to self-consistency [5]. However, oneusually stops after the first iteration. The correspondingapproximation for the Green function is known as the G W approximation and has become one of the mostaccurate methods for the calculation of spectral proper-ties and band gaps of solids [3, 4]. One reason for notgoing beyond the first iteration of the G W method isthe large computational cost involved. There are fur-ther indications that a full self-consistent solution wouldworsen the spectral properties as a consequence of a can-cellation between dressing of Green functions and vertexcorrections [7]. This was investigated for the electrongas [8] and the Hubbard model [9]. However, this prob-lem has not been investigated in detail for real systemsmainly due to the computational cost involved.The G W approximation has, however, two unsatisfac-tory aspects. The first aspect is related to the satisfactionof conservation laws. Baym [10] has shown that the self-energy expressions that can be obtained as a functionalderivative of a functional Φ[ G ] of the Green function,i.e. Σ = δ Φ /δG , have the important property that theylead to conserving many-body approximations. Theseapproximations obey basic conservation laws, like theones for particle number, momentum, angular momen-tum and energy. The GW A is one of these conservingschemes [11, 12, 13]. However, the Φ-derivable approx-imations are only conserving when the Dyson equationfor the Green function is solved fully self-consistently.A lack of full self-consistency will generally result in aviolation of the conservation laws. For this reason theuse of conserving approximations, such as GW , is crucialin obtaining a correct description of transport phenom-ena within a nonequilibrium Green function approach[14, 15, 16, 17, 18]. Since it is one of our research goalsto study quantum transport, it will be necessary to con-sider the fully self-consistent GW (SC- GW ) approxima-tion [8, 19, 20, 21, 22, 23, 24].A second unsatisfactory aspect of nonself-consistentschemes, such as G W , is that the values of the ob-servables depend on the way they are calculated. Forinstance, the total energy can be calculated in differentways from the Green function and the self-energy: usingthe Galitskii-Migdal formula [25], a coupling constant in-tegration [22], a Luttinger-Ward expression [13, 26, 27,28] or various other expressions. For nonself-consistentcalculations all these expressions lead to different resultsand therefore to ambiguity in the value of the energy.It was, however, demonstrated in the work of Baym [10]that self-consistent Φ-derivable approximations are notonly conserving but also have the property that all thevarious ways in which the observables are calculated pro-vide the same result. This is another motivation for con-sidering fully self-consistent many-body schemes.We can therefore conclude that self-consistency is im-portant to obtain conserving and unambiguous results.However, the large computational cost of self-consistentschemes makes them unattractive for the calculation ofthe properties of large and extended systems. In orderto lower the computational effort it is possible to usepartial self-consistency which may result in a less severeviolation of conservation laws. One can, for instance,keep the screened interaction fixed during iteration ofthe Dyson equation. This leads to a scheme that canbe shown to still conserve the particle number and thathas been tested on the electron gas [23, 29]. Anotherapproach in which the self-consistency is constrained isthe so-called quasi-particle self-consistent GW (QS GW )method [30, 31, 32, 33]. In this approach a frequencyindependent self-energy of GW -form is constructed andused to solve a quasi-particle equation from which theGreen function and the screened interaction are con-structed iteratively. Due to the Hermitian nature of theself-energy the method leads to an orthonormal set ofquasi-particle states and thereby restricts the form ofthe Green function and the screened interaction. Thismethod has been succesful in improving the G W bandgaps and band widths for a large range of solids [32].One could further consider similar other approximationswithin a quasi-particle framework [34]. Such approxi-mations have been shown to improve the band struc-ture when local density approximation is a poor start-ing point. These methods are, however, not Φ-derivableand are in general not conserving. Extending methodsbased on quasi-particle equations to the time-dependentcase is not as straightforward as for the SC- GW , GW and G W methods, which are instead based on an equa-tion of motion for the Green function. For the same rea-son the computational schemes used in this paper (whichaims at an extension to the time-dependent case) wouldneed to be modified in order to do QS GW calculations.We therefore did not consider the QS GW method in this work. However, we propose another partially self-consistent scheme which is computationally cheaper thanthe GW method. In this approximation the correlationpart of the self-energy is fixed during the iteration cyclewhile only the Hartree and exchange parts are updatedself-consistently. In this paper we investigate this ap-proximation and other GW schemes at different levelsof self-consistency and test them on atoms and diatomicmolecules. We also present in more detail the computa-tional method behind the self-consistent GW calculationsthat we described briefly in an earlier Letter [35]. The pa-per is divided as follows: In Sec. II we briefly present thegeneral formalism and in Sec. III we describe in detail the GW approximation at different levels of self-consistency.We then present in Sec. IV the details of our compu-tational procedure. Finally, in Sec. V, we will discussthe results obtained with the GW A at different levels ofself-consistency for atoms and some diatomic molecules.These systems are well-suited to test the GW at differ-ent levels of self-consistency, but we are ultimately in-terested in applications in quantum transport theory formolecules attached to macroscopic leads. In such appli-cations the long range screening effects, as incorporatedin the GW A, are important. The investigations in thispaper are a first step in this direction and aim to getfurther insight into various aspects of the GW A that arerelevant in quantum transport theory.
II. GENERAL FORMALISM
We study finite many-particle systems using the Mat-subara formalism [1, 36] which can easily be extended toa nonequilibrium version of the theory [37, 38, 39]. Weconsider a many-body system in thermal equilibrium ata temperature T and chemical potential µ , and with theHamiltonian (in second quantization [1])ˆ H = Z d x ˆ ψ † ( x ) h ( r ) ˆ ψ ( x ) ++ 12 Z Z d x d x ˆ ψ † ( x ) ˆ ψ † ( x ) v ( r , r ) ˆ ψ ( x ) ˆ ψ ( x ) . (1)Here x = ( r , σ ) denotes the space- and spin coordinates.The two-body interaction v is taken to be of Coulombicform v ( r , r ) = 1 / | r − r | . We use atomic units ~ = m = e = 1 throughout this paper. The single particlepart of the Hamiltonian h ( r ) has the explicit form h ( r ) = − ∇ + w ( r ) − µ, (2)where w ( r ) is the external potential and where we ab-sorbed the chemical potential µ into h . The equilibriumexpectation value of an operator ˆ O in the grand canonicalensemble is then given by h ˆ O i = Tr { ˆ ρ ˆ O } , (3)where ˆ ρ = e − β ˆ H / Tr e − β ˆ H is the statistical operator, β = 1 /k B T the inverse temperature and k B is the Boltz-mann constant. The trace is taken over all states in Fockspace [1]. The Green function is then defined as G ( x τ , x ′ τ ) = − θ ( τ − τ ) h ˆ ψ H ( x τ ) ˆ ψ † H ( x ′ τ ) i + θ ( τ − τ ) h ˆ ψ † H ( x ′ τ ) ˆ ψ H ( x τ ) i , (4)where we define the Heisenberg form of the operators inthis equation to be ˆ O H = e τ ˆ H ˆ Oe − τ ˆ H . Since the Hamilto-nian is time-translation invariant, the equilibrium Greenfunction only depends on the difference between the timecoordinates: G ( x τ , x ′ τ ) = G ( x , x ′ ; τ − τ ). The Greenfunction satisfies the equation of motion h − ∂ τ − h ( r ) i G ( x , x ′ ; τ ) == δ ( τ ) δ ( x − x ′ ) + Z β dτ Z d x Σ[ G ]( x , x ; τ − τ ) G ( x , x ′ ; τ ) , (5)where the self-energy Σ[ G ]( x , x ′ ; τ ) incorporates themany-body interactions of the system. The self-energycan be approximated with the usual diagrammatic meth-ods [1, 2]. Since Σ[ G ] is a functional of the Green functionEq.(5) must be solved self-consistently. The self-energyis usually split into a Hartree part and an exchange-correlation part, according toΣ[ G ]( x , x ; τ ) = δ ( τ ) δ ( x − x ) v H ( r )+Σ xc [ G ]( x , x ; τ ) , (6)where the Hartree potential is defined as the potentialdue to the electron charge by v H ( r ) = Z d x ′ n ( x ′ ) v ( r , r ′ ) , (7)where we introduced the electron density n ( x ) = lim η → G ( x , x ; − η ) . (8)The main task is now to find an approximation for thisexchange-correlation part Σ xc of the self-energy and tosolve Eq.(5). We convert Eq.(5) to integral form [40] G ( x , x ; τ ) = G ( x , x ; τ )+ Z β dτ dτ Z d x d x G ( x , x ; τ − τ ) × (Σ[ G ]( x , x ; τ − τ ) − δ ( τ − τ )Σ ( x , x )) × G ( x , x ; τ ) . (9)Here we introduced a static reference self-energy Σ anda reference Green function G which is defined by theequation h − ∂ τ − h ( r ) i G ( x , x ′ ; τ ) == δ ( τ ) δ ( x − x ′ ) + Z d x Σ ( x , x ) G ( x , x ′ ; τ ) . (10) Φ GW = − (cid:1) − (cid:2) − (cid:3) + . . . Σ GW = (cid:4) + (cid:5) + (cid:6) + . . . FIG. 1: The GW self-energy Σ is the functional derivativeof a functional Φ[ G ]. In practice we solve first Eq.(10) for G and then wesolve Eq.(9) for G . It is clear from Eq.(5) that a fullyself-consistent solution of Eq.(9) does not depend on thereference Green function G . In this work we choosefor Σ a Hartree-Fock (HF) or a density functional self-energy. In the first case Σ = v H [ G ]+Σ x [ G ], consistingof Hartree and exchange parts, whereas in the second caseΣ = δ ( x − x ′ ) v Hxc [ G ]( x ), where v Hxc ( x ) is the sum ofthe Hartree and the exchange-correlation potential [6].From the Green function several observables can becalculated. To calculate the total energy E = T + V ne + U + U xc we use the fact that the exchange-correlationpart U xc of the interaction energy is given by [1, 2] U xc = 12 Z β dτ Z d x Z d x Σ xc ( x , x ; − τ ) G ( x , x ; τ ) . (11)The kinetic energy T , the nuclear-electron attrac-tion energy V ne , and the Hartree energy U =1 / R d r d r ′ n ( r ) v ( r , r ′ ) n ( r ′ ) can all be calculated directlyfrom the Green function. To calculate the ionization po-tentials from the Green function we used the extendedKoopmans theorem [41, 42, 43, 44, 45], a short deriva-tion of which is given in Appendix B. III. THE GW APPROXIMATION ATDIFFERENT LEVELS OF SELF-CONSISTENCYA. Fully self-consistent GW
Within the GW A the exchange-correlation part of theself-energy has the explicit form [5, 46, 47]Σ xc ( x , x ; τ ) = − G ( x , x ; τ ) W ( x , x ; τ ) , (12)in which W is a dynamically screened interaction corre-sponding to an infinite summation of bubble diagrams(see Fig. 1). From this figure we see that this self-energyis given as a functional derivative of a functional Φ[ G ]with respect to G and hence represents a conserving ap-proximation [10]. From the diagrammatic structure wesee that the screened potential W satisfies the equation W ( x , x ; τ ) = v ( r , r ) δ ( τ ) ++ Z d x d x Z β dτ ′ v ( r , r ) P ( x , x ; τ − τ ′ ) W ( x , x ; τ ′ ) , (13)where v is the bare Coulomb interaction and P is theirreducible polarization P ( x , x ; τ ) = G ( x , x ; τ ) G ( x , x ; − τ ) . (14)The problem is now completely defined. Equations (13)and Eq.(14) need to be solved self-consistently togetherwith Eqs.(12), (6) and (9). B. The G W and GW approximations The G W approximation, as mentioned before, is ob-tained from a single iteration of the Dyson equationEq.(9), starting from a refence Green function G . Forthis approximation the self-energy is given as Σ xc [ G ] = − G W where W is calculated by inserting G intoEq.(14) and solving Eq.(13) with this irreducible polar-ization. The Dyson equation (9) is then solved with thisself-energy to obtain an improved Green function G fromwhich spectral properties are calculated. In principle oneshould insert this Green function into the self-energy andsolve the Dyson equation again for a new Green function.This procedure should be continued until self-consistencyis achieved, but this is rarely done in practice for the rea-sons mentioned in the introduction.We further consider a partially self-consistent scheme inwhich we write the self-energy as Σ xc [ G, G ] = − GW ,where the Green function G is determined fully self-consistently by repeated solution of the Dyson equa-tion and where W is calculated from G in the sameway as for the G W approximation. This reduces thecomputational cost considerably as it avoids the self-consistent calculation of the screened interaction W . Thecorresponding approximation is known as the GW ap-proximation [29, 48]. This approximation was shownto be number conserving by Holm and von Barth [49]for the case of homogeneous systems. More preciselythey derived that the GW approximation satisfies theHugenholtz-van Hove theorem [50] for the homogeneouselectron gas. However, one can readily derive the num-ber conserving property for the inhomogeneous and time-dependent case. This requires nonequilibrium Greenfunctions in the proof, but this extension is straightfor-ward [51]. If we regard W as a given potential (albeitnonlocal in space and time), it is clear that Σ = δ Φ /δG for Φ[ G, W ] = − / GGW , where the trace denotesintegration over space-time variables. Since this Φ is in-variant under gauge transformations (the phases cancelat each vertex of Φ), we can follow the proof of Baym [10]and derive that GW is particle conserving. However, fortime-dependent and inhomogeneous systems W is notinvariant under spatial and time-translations, unlike thebare interaction v that usually appears in the functionalΦ[ G ]. Therefore the GW approximation will not be mo-mentum or energy conserving. C. The GW fc approximation The most time-consuming part of the GW calculationis the evaluation of the correlation part of the self-energywhich is nonlocal in time. We therefore propose anotherpartial self-consistent scheme in which we only evaluatethe time-local Hartree and exchange parts of the self-energy in a self-consistent manner. We therefore splitthe self-energy as followsΣ[ G, G ] = Σ HF [ G ] + Σ c [ G ] . (15)The first term in this equation represents the Hartree-Fock part of the self-energyΣ HF [ G ] = v H [ G ] + Σ x [ G ] , (16)which consists of a Hartree part and an exchange partΣ x [ G ] = − Gv . The last term in Eq.(15) represents thecorrelation part of the self-energy and has the explictform Σ c [ G ] = − G ( W − v ) , (17)where W is calculated from G in the same way as forthe G W approximation. The approximation for theself-energy of Eq.(15) will be denoted as the GW fc ap-proximation (where fc stands for fixed correlation). Thisapproximation is not conserving but, as we will see later,nevertheless produces observables in very close agree-ment with those obtained from a fully SC- GW calcu-lation. IV. COMPUTATIONAL METHODA. Numerical solution of the Dyson equation
In the following, we will describe the computationalmethods that we employed for calculating the Greenfunction and the screened interaction W . We considerthe case of spin-unpolarized systems where the Greenfunction has the form G ( x , x ′ ; τ ) = δ σσ ′ G ( r , r ′ ; τ ) . (18)The calculations are carried out using a set of basis func-tions such that the spin-independent part of the Greenfunction is expressed as G ( r , r ′ ; τ ) = X ij G ij ( τ ) φ i ( r ) φ ∗ j ( r ′ ) . (19)The basis functions φ i are represented as linear combi-nations of Slater functions ψ i ( r ) = r n i − e − λ i r Y m i l i (Ω)which are centered on the different nuclei and are charac-terized by quantum numbers ( n i , l i , m i ) and an exponent λ i . In these expressions and Y m i l i (Ω) are the usual spheri-cal harmonics. The molecular orbitals φ i and eigenvalues ǫ i are obtained from a Hartree-Fock or DFT Kohn-Shamcalculation in this basis. The particle number N is de-termined by the chemical potential. Since we considerclosed shell systems we have N/ ǫ i (some of which may be degener-ate). We therefore choose µ such that e i = ǫ i − µ < i ≤ N/ e i > i > N/
2. In the zero-temperaturelimit (we used β = 100) the observables are insensitiveto the value of µ , provided ǫ N/ < µ < ǫ N/ . Thereference Green function G corresponding to the Hamil-tonian h +Σ (either HF or DFT) is diagonal in the basis { φ i } i.e. in matrix form we have G ij, ( τ ) = δ ij G i, ( τ ),where G i, ( τ ) = θ ( τ )( n ( e i ) − e − e i τ + θ ( − τ ) n ( e i ) e − e i τ , (20)and n ( e j ) = ( e βe j + 1) − is the Fermi-Dirac distribution.The Dyson equation of Eq.(9) in basis representation hasthe form G ( τ ) = G ( τ )+ Z β dτ ′ Z β dτ ′′ G ( τ − τ ′ )Σ c [ G, G ]( τ ′ − τ ′′ ) G ( τ ′′ ) , (21)where we denoteΣ c [ G, G ]( τ ) = Σ[ G ]( τ ) − δ ( τ )Σ [ G ] , (22)and where all quantities are matrices. Since in the limit τ → − , G yields the density matrix, it is convenientto solve the Dyson equation for negative τ –values. Wetherefore rewrite Eq. (21) as G ij ( τ ) = δ ij G i, ( τ )+ X k Z − β dτ Z − β dτ G i, ( τ − τ )Σ cik ( τ − τ ) G kj ( τ ) , (23)with τ ∈ [ − β,
0] where we changed variables τ = τ ′ − β , τ = τ ′′ − β , and used G ( τ ) = − G ( τ + β ) with the samerelation for G [40]. We now discretize Eq. (23) using atrapezoidal rule on a time grid ( τ (0) = 0 , τ (1) . . . , τ ( m ) = − β ). Since the Green functions behave exponentiallynear the endpoints of the imaginary time interval [ − β, δ ij G i, ( τ ( p ) ) = X k,q (cid:20) δ ik δ pq − ∆ τ ( q ) Z ik ( τ ( p ) , τ ( q ) ) (cid:21) G kj ( τ ( q ) ) , (24)where we defined Z ik as Z ik ( τ ( p ) , τ ( q ) ) = Z − β dτ G i, ( τ ( p ) − τ )Σ cik ( τ − τ ( q ) ) . (25)The time steps are positive, where ∆ τ ( q ) = τ ( q − − τ ( q +1) except at the endpoints where ∆ τ (0) = τ (0) − τ (1) and ∆ τ ( m ) = τ ( m − − τ ( m ) . For a fixed j , Eq. (24) representsa set of linear equations of the form X Q A Q ,Q · x ( j ) Q = b ( j ) Q , (26)where A Q ,Q = A ( ip )( kq ) = δ ik δ pq − ∆ τ ( q ) Z ik ( τ ( p ) , τ ( q ) )and the vectors x ( j ) Q , b ( j ) Q are defined to be x ( j ) Q = x ( j ) kq = G kj ( τ ( q ) ) b ( j ) Q = b ( j ) ip = δ ij G i, ( τ ( p ) ) . The self-energy Σ c of Eq.(22) has the formΣ cij ( τ ) = Σ c,ij [ G ]( τ ) + δ ( τ ) (cid:2) Σ HFij [ G (0 − )] − Σ ij (cid:3) , (27)where Σ HF is the Hartree-Fock part of the self-energydefined in Eq.(16) and Σ c [ G ] the remaining correlationpart. The convolution integral (25) can therefore be sim-plified to Z ik ( τ ( p ) , τ ( q ) ) = G i, ( τ ( p ) − τ ( q ) ) (cid:2) Σ HFik [ G (0 − )] − Σ ik (cid:3) + Z − β dτ G i, ( τ ( p ) − τ )Σ c,ik ( τ − τ ( q ) ) . (28)When we specify the explicit form of Σ c , the solution ofthe Dyson equation is reduced to a calculation of Eq.(28)together with the linear system of equations (26). Whatremains to be discussed is the calculation of the self-energy itself. This is discussed in the next section. B. Numerical calculation of the screened potential:The product basis technique
To calculate the self-energy we need to solve the equa-tion for the screened interaction. The screened interac-tion has a singular time-local part representing the bareinteraction v . It is therefore convenient to subtract v from W and to treat its contribution to the self-energyexplicitly (this is simply the exchange part of the self-energy). From the remaining time nonlocal part of W ,given by f W ( r , r ; τ ) = W ( r , r ; τ ) − δ ( τ ) v ( r , r ), wecan calculate the correlation part of the self-energyΣ c ( r , r ; τ ) = − G ( r , r ; τ ) f W ( r , r ; τ ) . (29)After this quantity has been calculated it can then simplybe added to the Hartree-Fock part of the self-energy toobtain the full self-energy Σ[ G ]. The time-nonlocal part f W of the screened interaction satisfies the equation f W ( r , r ; τ ) = Z d r d r v ( r , r ) P ( r , r ; τ ) v ( r , r ) ++ Z β dτ ′ Z d r d r v ( r , r ) P ( r , r ; τ − τ ′ ) f W ( r , r ; τ ′ )(30)where P ( r , r ; τ ) = 2 G ( r , r ; τ ) G ( r , r ; − τ ) . (31)The factor of 2 in this expression results from spin-integrations in the equation of W using the form of theGreen function of Eq.(18). We now insert into Eq.(31)the basis set expansion for the Green function of Eq.(19),to obtain P ( r , r ; τ ) = X ijkl P ijkl ( τ ) φ i ( r ) φ ∗ j ( r ) φ k ( r ) φ ∗ l ( r )(32)where P ijkl = 2 G ij ( τ ) G kl ( − τ ). By defining the two-electron integrals f W pqrs ( τ ) = Z d r d r φ ∗ p ( r ) φ ∗ q ( r ) f W ( r , r ; τ ) φ r ( r ) φ s ( r ) v pqrs = Z d r d r φ ∗ p ( r ) φ ∗ q ( r ) v ( r , r ) φ r ( r ) φ s ( r )we transform Eq.(30) into the equation f W pqrs ( τ ) = X ijkl v plis P ijkl ( τ ) v jqrk + X ijkl Z β dτ ′ v plis P ijkl ( τ − τ ′ ) f W jqrk ( τ ′ ) . (33)If we use the multi-indices Q = ( ps ), Q = ( rq ), Q =( il ) and Q = ( jk ), then we can write this equation in amore convenient form as f W Q Q ( τ ) = X Q Q v Q Q P Q Q ( τ ) v Q Q + X Q Q Z β dτ ′ v Q Q P Q Q ( τ − τ ′ ) f W Q Q ( τ ′ ) , (34)where we defined v il,kj = v ijkl and similarly for f W and P il,jk = P ijkl . We have now obtained an equation whichwe can solve with the same algorithm we used for theDyson equation.Note that in this case we effectively use a product ba-sis f q ( r ) = φ i ( r ) φ ∗ j ( r ), where q = ( ij ) is a multi-index.This product basis is nonorthogonal and its size is in gen-eral much larger than we need in practice due to lineardependencies. We thus follow a technique developed byAryasetiawan and Gunnarsson [52], which allows to re-duce significantly the size of the product basis { f q ( r ) } and the computational cost.The overlap matrix S for the set of orbitals f q ( r ) S qq ′ = h f q | f q ′ i , (35)is diagonalized by a unitary matrix U X q q U † qq h f q | f q i U q q ′ = σ q δ qq ′ , (36) where the eigenvalues σ q are positive since S is a positivedefinite matrix. We now define a new set of orthonormalorbitals g q as g q ( r ) = 1 √ σ q X q ′ U q ′ q f q ′ ( r ) , (37)with h g q | g q ′ i = δ qq ′ . Our strategy is use the orbitals g q as a new basis and discard the functions that correspondto σ q < ǫ (we used ǫ = 10 − ). This leads to a much re-duced basis as compared to the set of all functions f q . Asdescribed in Ref. 52, this corresponds to discarding func-tions that are nearly linearly dependent and contributelittle in the expansion. The quantities Σ, f W and P willbe represented in this new basis using f q ( r ) = X q ′ g q ′ ( r ) √ σ q ′ U † q ′ q . (38)For the irreducible polarization we then find fromEq. (32) that P ( r , r ; τ ) = X qq ′ P qq ′ ( τ ) f q ( r ) f ∗ q ′ ( r ) == X q q X qq ′ U † q q P qq ′ ( τ ) U q ′ q √ σ q σ q g q ( r ) g ∗ q ( r ) , (39)where P qq ′ ( τ ) = 2 G ij ( τ ) G kl ( − τ ) . (40)With q = ( il ) and q ′ = ( jk ) we have P ( r , r ; τ ) = X q q e P q q g q ( r ) g ∗ q ( r ) , (41)where e P q q = (cid:2) √ σU † P ( τ ) U √ σ (cid:3) q q , (42)and √ σ is the diagonal matrix ( √ σ ) pq = δ pq √ σ q . Tocalculate the screened potential we now insert Eq.(41)into Eq.(30) and readily obtain the matrix product f W qq ′ ( τ ) = h v e P ( τ ) v i qq ′ + h v e P ( τ − τ ′ ) f W ( τ ′ ) i qq ′ , (43)where we defined the matrices f W qq ′ = Z d r d r g ∗ q ( r ) f W ( r , r ; τ ) g q ′ ( r ) (44)and v qq ′ = Z d r d r g ∗ q ( r ) v ( r , r ) g q ′ ( r ) . (45)It is important to note that in Eq.(41) and Eq.(43) thesummation only runs over the indices q for which σ q > ǫ .We see from Eq.(42) that terms with σ q < ǫ contributelittle to the total sum. This leads to a considerable reduc-tion of the number of matrix elements for v , P and f W .Finally the correlation part of the self-energy of Eq.(29)is given byΣ c,ij ( τ ) = Z d r Z d r φ ∗ i ( r )Σ c ( r , r ; τ ) φ j ( r )= − X kl G kl ( τ ) X pq f W pq ( τ ) Z d r φ ∗ i ( r ) φ k ( r ) g p ( r ) × Z d r φ j ( r ) φ ∗ l ( r ) g ∗ q ( r )= − X kl G kl ( τ ) Z ik,jl , (46)where Z ik,jl = X pq √ σ p U ik,p f W pq ( τ ) U † q,jl √ σ q . (47)We can summarize our procedure as follows: in the firststep the overlap matrix S qq ′ of Eq.(35) is obtained and di-agonalized. Further, using and Eq.(37) and (45) the two-electron integrals in the new basis v pq are constructed for p and q such that σ p , σ q > ǫ . Subsequently, for the samevalues of p and q the matrix e P pq ( τ ) is constructed fromEq.(42) and f W pq ( τ ) is solved from Eq.(43). In the laststep, the matrix (47) is obtained and the self-energy iscalculated from Eq.(46) and further used in the solutionof the Dyson equation. V. RESULTS
The various GW schemes described in section III areapplied to a set of atoms and diatomic molecules usingthe computational method of section IV. Details on thebasis sets are provided in Ref. [53]. In general we foundthat, in single processor calculations, the computationalcost of the GW fc method is comparable to that of the G W method, and roughly twice as fast as the GW method. The the fully self-consistent GW calculationswere the most time-consuming. Particle number conservation . We start by investigat-ing the number conservation property of the different GW schemes. In Fig. 2 we display the particle num-ber obtained from the trace of the Green function forthe case of the hydrogen molecule H for different sep-arations of the nuclei. We display results for the caseof SC- GW , GW , GW fc and G W , in which the refer-ence Green function G is obtained from a Hartree-Fockcalculation. We see that the SC- GW and GW schemesyield an integer particle number of N = 2 for all internu-clear separations. This is a consequence of the numberconserving property of both approximations. This canbe seen as follows. If we would adiabatically switch-onthe two-particle interactions from zero to full couplingstrength within a conserving scheme then the particle P a r ti c l e nu m b e r GWGW GW fc G W FIG. 2: Particle number for H at different interatomic dis-tances within the SC- GW , GW , GW fc and G W approxi-mations. number would be conserved during the switching. Thisis because the conserving property is independent of thestrenght of the interaction and follows from the structureof the Φ-functional only. Therefore the particle numberof the final correlated state will be the same as the parti-cle number of the initially noninteracting system. Henceconserving schemes always yield integer particle numberfor finite systems at zero temperature. For the case ofthe hydrogen molecule this is N = 2 for all bond dis-tances. For the case of G W we see that the particlenumber conservation is violated as the particle numberdeviates from N = 2 for all bond distances, the largestdeviations occuring for the larger bond distances. Forthe larger separations left-right correlation [54] in the hy-drogen molecule, not incorporated in the Hartree-Fockpart of the self-energy, become increasingly important.This puts more demands on the quality of the correla-tion part of the self-energy and consequently nonconser-vation of the particle number becomes more apparentat longer bond distances. Although the violation seemssmall (about 0 .
01 electron at R = 4 .
5) it should be em-phasized that a change in particle number of 0 .
05 cangive large changes in the spectral features and conduc-tive properties for molecules attached to leads. A clearexample of this is presented in the work of Thygesen [16].For the GW fc (See sec. III C) we also observe a violationof the number conservation law with increasing error forlarger internuclear separations. The error with respectto G W is however reduced by a factor of 3 at R = 5 . Ground state energies . For the various GW schemes ofsection III we calculated the total energies of some atomsand diatomic molecules from Eq.(11). The referenceGreen function G for the nonself-consistent schemes wasobtained from a Hartree-Fock calculation. In Table Iwe show the results. From comparison with benchmarkconfiguration interaction (CI) results we see that the to-tal energies of atoms and molecules calculated within allschemes are not very accurate. However, as we will seelater, energy differences are much better produced. Wecan nevertheless make a number of useful observationsfrom the total energies. We first note that all approx-imations produce a total energy that is lower than thebenchmark CI result, with the G W generally produc-ing the lowest and thereby the worst values. Both the GW and the GW fc methods yield total energies in ex-cellent agreement with SC- GW results, where for mostsystems the difference is 10 − Hartree or less. This meansthat both the GW and the GW fc methods can be usedto make an accurate prediction for the SC- GW energyat a much lower computational cost than the fully self-consistent calculation. Binding curve . The calculation of binding curves is agood test for the quality of total energy calculations. InFig. 3 we display the binding curve of the H moleculefor the various GW schemes together with benchmarkCI results. The reference Green function G was takenfrom a Hartree-Fock calculation. We further checked thatusing a G obtained from an LDA calculation only influ-ences the results slightly. For the values of the energiesaround the bond minimum we see the same trend thatwe observed before: all GW schemes lead to a total en-ergy that is lower than the benchmark CI results with G W being the lowest. The total energies of the par-tially self-consistent schemes GW and GW fc are veryclose to the fully self-consistent GW results for all bonddistances. Although all GW schemes considerably im-prove the bonding curve obtained from an uncorrelatedHartree-Fock calculation it is clear that all these schemesdeviate considerably from the CI results in the infiniteatomic separation limit. To cure this feature one eitherhas to do a spin-polarized calculation or go beyond the GW approximation and include vertex diagrams in thediagrammatic expansion for the self-energy. The shapeof the binding curve around the bond minimum is wellreproduced by the SC- GW , GW and GW fc schemes,implying that these methods may be used to obtain ac-curate vibrational frequencies. Since the shape of thebonding curve is only determined by total energy differ-ences, this already indicates that these approximationsmay perform better in obtaining the energy differencesthan in obtaining total energies. Two-electron removal energies . To test the perfor-mance of the various GW schemes in obtaining energydifferences, we investigated the two-electron removal en-ergies of the beryllium and magnesium atom. Since theseatoms and their doubly ionized counterparts are closedshell they were suitable test systems. Moreover, theberyllium atom is a well-known case for which electroncorrelations play an important role due to strong mix-ing of the 2 s and 2 p states in a configuration expansion.In table II, we display the two-electron removal ener-gies for various GW schemes as well as for the Hartree-Fock approximation. The reference Green function G is again obtained from a Hartree-Fock calculation. Theself-consistent and partially self-consistent GW schemes E n e r gy ( a . u . ) CIGWGW GW fc G W FIG. 3: The total energy of the H molecule, as a functionof the interatomic distance, calculated from the GW approx-imation at various levels of self-consistency and CI [56]. yield results within 0 . G W approxi-mation does not improve at all on the HF approximationand gives considerably worse results than the other GW schemes. We further see that both the GW and the GW fc approximations give removal energies that are inexcellent agreement with the fully self-consistent GW re-sults. Ionization Potentials . In Table III we show the ioniza-tion potentials obtained with the various GW methodsfor a number of atoms and diatomic molecules. Theseionization potentials were obtained using the extendedKoopmans theorem, as explained in Appendix A. For G W the results shown in the first column were obtainedby using a reference Green function G from a local den-sity functional (LDA) calculation using the parametriza-tion of the exchange-correlation functional due to Voskoet al. [59]. In all other cases we used a reference Greenfunction from a Hartree-Fock calculation. We see thatthe ionization potentials of fully self-consistent GW agreewell with the experimental values, the main exceptionsbeing the H molecule and the Be atom, which show adeviation of respectively 0 . . GW and GW fc yield re-sults that are very close to the fully self-consistent re-sults. The G W approximation based on the LDA ref-erence Green function performs a bit worse than the self-consistent GW scheme. For He and LiH there is an errorof about 1 eV and for Ne and H an error of about 0 . G W calculation based on a HF refer-ence G instead improves the results for several systemsbut worsens the agreement for H which is 1 . G within the G W method is clearly unsatisfactory. Thepartially self-consistent approximations suffer much lessfrom this problem. For those schemes we found thatchanging the reference Green function from a HF one to TABLE I: Total energies (in Hartrees) calculated from the GW approximation at various levels of self-consistency comparedto CI values.System E G W [ G HF ] E GW [ G HF ] E GW fc [ G HF ] E GW SC CIHe -2.9354 -2.9271 -2.9277 -2.9278 -2.9037 Be -14.7405 -14.6882 -14.7032 -14.7024 -14.6674 Be -13.6929 -13.6886 -13.6887 -13.6885 -13.6556 Ne -129.0885 -129.0517 -129.0506 -129.0499 -128.9376 Mg -200.2924 -200.1759 -200.1775 -200.1762 -200.053 Mg -199.3785 -199.3451 -199.3454 -199.3457 -199.2204 H -1.1985 -1.1889 -1.1891 -1.1887 -1.133 LiH -8.1113 -8.0999 -8.0997 -8.0995 -8.040 From Ref. [55]. From Ref. [56]. From Ref. [57].TABLE II: Two-electron removal energies E N − − E N (in eV)calculated from the Hartree-Fock and from the GW approx-imation at various levels of self-consistency, compared to theexperimental values.System HF ∆ E G W ∆ E GW ∆ E GW fc ∆ E GWSC
Expt. Mg - Mg From Ref. [58].TABLE III: Ionization potentials (eV) calculated from theextended Koopmans theorem from various GW approaches.Sys. G ( LDA )0 W G ( HF )0 W GW GW fc SC-GW Expt. He 23.65 24.75 24.59 24.56 24.56 24.59Be 8.88 9.19 8.82 8.81 8.66 9.32Ne 21.06 21.91 21.90 21.82 21.77 21.56Mg 7.52 7.69 7.43 7.38 7.28 7.65H From Ref. [58] an LDA one, only slightly changes the results.
VI. SUMMARY AND CONCLUSIONS
We investigated the performance of the GW at differ-ent levels of self-consistency for the case of atoms anddiatomic molecules. Our main motivation for studyingfully self-consistent Φ-derivable schemes was that theyprovide unambiguous results for different observables andthe fact that they satisfy important conservation lawsthat are important in future nonequilibrium applicationsof the theory [18]. We adressed the question to what ex-tent partially self-consistent schemes can reproduce theresults of a fully self-consistent GW calculation. Wefound that both the GW method, as well as the GW fc scheme proposed by us, yield results in close agreementwith fully self-consistent GW calculations. We furtherchecked the number conservation properties of the var-ious schemes. The fully self-consistent GW scheme be-ing Φ-derivable does satisfy all conservation laws, butalso the partially self-consistent GW approximation wasshown to be number conserving. The nonself-consistent G W and the partially self-consistent GW fc approxima-tions both violate the number conservation laws but, dueto the partial self-consistency in GW fc , the errors aremuch reduced in this scheme. A major advantage of thelatter scheme is, however, that it produces results thatare close to the fully self-consistent GW results at a muchlower computational cost. It will therefore be very valu-able to test this method on solid state systems for whichself-consistent GW calculations are difficult to performdue to the large computational effort. In this way it willbe possible to get further insight into the performanceof self-consistent GW for a large class of extended sys-tems. Work on application of the fully self-consistent GW method to transport phenomena is in progress [18]. APPENDIX A: IONIZATION POTENTIALSFROM THE EXTENDED KOOPMANSTHEOREM
Here we give a brief description on the way we extractthe ionization energies from the Green function using theextended Koopmans theorem [41, 42, 43, 44, 45]. Asinput, this method only needs the Green function and itstime derivative at τ = 0 − on the imaginary time axis.We define an N − | Φ N − [ u i ] > = Z d x u i ( x ) ˆ ψ ( x ) | Ψ N >, (A1)where u i ( x ) is determined by requiring the functional E N − [ u i ] = h Φ N − [ u i ] | ˆ H | Φ N − [ u i ] ih Φ N − [ u i ] | Φ N − [ u i ] i , (A2)which describes the energy of the N − u i . This0amounts to minimizing the energy of the N − u i . We find Z d x h Ψ N | ˆ ψ † ( x ′ ) (cid:2) ˆ ψ ( x ) , ˆ H (cid:3) | Ψ N i u i ( x ) =( E N − E N − i ) Z d x h Ψ N | ˆ ψ † ( x ′ ) ˆ ψ ( x ) | Ψ N i u i ( x ) , (A3)where the last term contains the density matrix. Thisquantity is easily obtained from the Green function as ρ ( x , x ′ ) = h ψ N | ˆ ψ † H ( x ′ τ ) ˆ ψ H ( x τ ) | Ψ N i = lim η → G ( x , x ′ , − η )(A4) i.e. ρ ( x , x ′ ) = ˜ G ( x , x ′ ; 0 − ) or ρ ij = G ij (0 − ) in molecu-lar orbital basis [40]. Also the expectation value underthe integral on the righthand side of Eq.(A3), is easilyobtained from the Green function − ∂ τ G ( x , x ′ ; τ ) | τ =0 − = h Ψ N | ˆ ψ † ( x ′ ) (cid:2) ˆ ψ ( x ) , ˆ H (cid:3) | Ψ N i = ∆( x , x ′ ) . (A5)In this derivation we used a zero-temperature formula-tion but making a connection to the finite temperatureformalism is straightforward. When we take into accountthat, in the finite temperature formalism, we included thechemical potential in the one-body part of the Hamilto-nian (see Eq.(2), then from (A3) and (A5) we obtain theeigenvalue equation Z d x ∆( x , x ′ ) u i ( x ) == ( E N − E N − i − µ ) Z d x ρ ( x , x ′ ) u i ( x ) , (A6)where ρ and ∆ are calculated according to Eq.(A4, A5).A similar equation for the electron affinities can simi-larly be derived starting from an N + 1-state. Since bothmatrices ρ and ∆ are easily evaluated from the Greenfunction, Eq.(A6) provides an easy way to extract re-moval energies from knowledge of the Green function onthe imaginary time axis. For completeness we mention that the extended Koop-mans method also provides a simple way to extract quasi-particle or Dyson orbitals [45] and to construct the Greenfunction on the real frequency axis. The Dyson orbitalsare given by f i ( x ) = h Φ N − i | ˆ ψ ( x ) | Ψ N i == Z d x ′ u ∗ i ( x ′ ) h Ψ N | ˆ ψ † ( x ′ ) ˆ ψ ( x ) | Ψ N i == Z d x ′ ρ ( x , x ′ ) u ∗ i ( x ′ ) . (A7)In terms of these orbitals and the extended Koopmanseigenvalues the hole-part of the Green function is thengiven on the real frequency axis as G ( x , x ′ ; ω ) = X n f n ( x ) f ∗ n ( x ′ ) ω − ( E N − E N − n + µ ) + iη . Similar derivations can be carried out for the affinitiesand the corresponding Dyson orbitals from which theparticle-part of the Green function can be constructedon the real axis.
APPENDIX B: THE UNIFORM POWER MESH
The uniform power mesh (UPM) [20] is a one-dimensional grid on an interval [0 , β ] which becomes moredense at the endpoints. Therefore, it is well-suited to de-scribe the Green function on the imaginary time axis,since it behaves exponentially around τ = 0 and τ = ± β [20, 40]. The UPM is defined by two integers u and p andthe length of the interval β . The procedure to constructit is simple: we consider the 2( p −
1) intervals [0 , β j ] and[ β − β j , β ] for j = 1 , . . . , p − β j = β/ j , and divideeach of these intervals in 2 u subintervals of equal lenght.The endpoints of all these intervals define our grid whichhas 2 pu + 1 grid points. [1] A. L. Fetter and J. D. Walecka, Quantum Theory ofMany-Particle Systems (McGraw-Hill, New York, 1971).[2] E. K. U. Gross, E. Runge, and O. Heinonen,
Many-Particle Theory (Verlag Adam Hilger, Bristol, 1991).[3] F. Aryasetiawan and O. Gunnarsson, Rep. Prog. Phys , 237 (1998).[4] W. Aulbur, L. J¨onsson, and J. Wilkins, Solid StatePhysics , 1 (2000).[5] L. Hedin, Phys. Rev. , A796 (1965).[6] R. M. Dreizler and E. K. U. Gross, Density FunctionalTheory, An Approach to the Quantum Many-Body Prob-lem (Springer-Verlag, Berlin, 1990).[7] G. D. Mahan, Comm. Cond. Mat. Phys. , 333 (1994).[8] B. Holm and U. von Barth, Phys. Rev. B , 2108 (1998). [9] T. J. Pollehn, A. Schindlmayr, and R. W. Godby, J.Phys.: Condens. Matter , 1273 (1998).[10] G. Baym, Phys. Rev. , 1391 (1962).[11] G. Baym and L. P. Kadanoff, Phys. Rev. , 287 (1961).[12] M. Bonitz, K. Balzer, and R. van Leeuwen, Phys. Rev.B , 045341 (2007).[13] C-.O.Almbladh, U. Barth, and R. Leeuwen, Int. J. Mod.Phys. B , 535 (1999).[14] K. S. Thygesen and A. Rubio, Phys. Rev. B , 115333(2008).[15] K. S. Thygesen and A. Rubio, J. Chem. Phys. ,091101 (2007).[16] K. S. Thygesen, Phys. Rev. Lett. , 166804 (2008).[17] X. Wang, C. D. Spataru, M. S. Hybertsen, and A. J. Millis, Phys. Rev. B , 045119 (2008).[18] P. My¨ohanen, A. Stan, G. Stefanucci, and R. vanLeeuwen, Europhys. Lett. , 67001 (2008).[19] A. Schindlmayr and R. W. Godby, Phys. Rev. Lett ,1702 (1998).[20] W. Ku and A. G. Eguiluz, Phys. Rev. Lett. , 126401(2002).[21] K. Delaney, P. Garc´ıa-Gonz´alez, A. Rubio, P. Rinke, andR. W. Godby, Phys. Rev. Lett. , 249701 (2004).[22] B. Holm, Phys. Rev. Lett. , 788 (1999).[23] P. Garc´ıa-Gonz´alez and R. W. Godby, Phys. Rev. B ,075112 (2001).[24] B. Holm and U. von Barth, Physica Scripta T109 , 135(2004).[25] V. M. Galitskii and A. B. Migdal, Zh. Eksp. Teor. Fiz. , 139 (1958), [ Sov. Phys. JETP , , 96 (1958)].[26] N. E. Dahlen, R. van Leeuwen, and U. von Barth, Phys.Rev. A , 012511 (2006).[27] R. van Leeuwen, N. E. Dahlen, and A. Stan, Phys. Rev.B , 195105 (2006).[28] M. P. Agnihotri, W. Apel, and W. Weller, Phys. Stat.Sol. B , 421 (2008).[29] U. von Barth and B. Holm, Phys. Rev. B , 8411 (1996).[30] T. Kotani and M. van Schilfgaarde, Solid St. Comm. ,461 (2002).[31] S. V. Faleev, M. van Schilfgaarde, and T. Kotani, Phys.Rev. Lett. , 126406 (2004).[32] M. van Schilfgaarde, T. Kotani, and S. Faleev, Phys. Rev.Lett. , 226402 (2006).[33] S. V. Faleev, M. van Schilfgaarde, T. Kotani, F. L´eonard,and M. P.Desjarlais, Phys.Rev. B , 033101 (2006).[34] F. Bruneval, N. Vast, and L. Reining, Phys.Rev. B ,045102 (2006).[35] A. Stan, N. E. Dahlen, and R. van Leeuwen, Europhys.Lett. , 298 (2006).[36] T. Matsubara, Progr. Theoret. Phys , 351 (1955).[37] L. V. Keldysh, Zh. Eksp. Teor. Fiz. , 1515 (1964), [ Sov.Phys. JETP , , 1018 (1965)].[38] P. Danielewicz, Ann. Phys. (N. Y.) , 239 (1984).[39] M. Wagner, Phys. Rev. B , 6104 (1996).[40] N. E. Dahlen and R. van Leeuwen, J. Chem. Phys. , 164102 (2005).[41] J. Katriel and E. R. Davidson, Proc. Natl. Acad. Sci.,USA , 4403 (1980).[42] O. W. Day, D. W. Smith, and R. C. Morrison, J. Chem.Phys. , 115 (1975).[43] D. W. Smith and O. W. Day, J. Chem. Phys. , 113(1975).[44] D. Sundholm and J. Olsen, J. Chem. Phys. , 3999(1993).[45] R. C. Morrison and P. W. Ayers, J. Chem. Phys. ,6556 (1995).[46] L. Hedin, A. Johansson, B. I. Lundqvist, S. Lundqvist,and V. Samathiyakanit, Arkiv f¨or Fysik , 97 (1968).[47] H. N. Rojas, R. W. Godby, and R. J. Needs, Phys. Rev.Lett. , 1827 (1995).[48] C. Fortmann, J. Phys. A: Math. Theor. , 445501(2008).[49] B. Holm, Ph.D. thesis, Lund University, Sweden (1997).[50] N. H. Hugenholtz and L. van Hove, Physica , 363(1958).[51] N. E. Dahlen, A. Stan, and R. van Leeuwen, J. Phys.Conf. Ser. , 324 (2006).[52] F. Aryasetiawan and O. Gunnarsson, Phys. Rev. B ,16214 (1994).[53] See EPAPS document for basis sets.[54] E. J. Baerends and O. V. Gritsenko, J. Phys. Chem. A , 5383 (1997).[55] S. J. Chakravorty, S. R. Gwaltney, E. R. Davidson, F. A.Parpia, and C. F. Fischer, Phys. Rev. A , 3649 (1993).[56] R. van Leeuwen, Ph.D. thesis, Vrije Universiteit, Ams-terdam (1994).[57] X. Li and J. Paldus, J. Chem. Phys. , 2470 (2003).[58] S.G. Lias, R.D. Levin, and S.A. Kafafi, Ion Energet-ics Data in NIST Chemistry WebBook, NIST Stan-dard Reference Database Number 69, Eds. P.J. Lin-strom and W.G. Mallard, March 2003, National Instituteof Standards and Technology, Gaithersburg MD, 20899(http://webbook.nist.gov).[59] S. H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys.58