Fully selfconsistent GW calculations for molecules
aa r X i v : . [ c ond - m a t . m t r l - s c i ] J a n Fully selfconsistent GW calculations for molecules
C. Rostgaard, K. W. Jacobsen, and K. S. Thygesen
Center for Atomic-scale Materials Design (CAMD),Department of Physics, Technical University of Denmark, DK - 2800 Kgs. Lyngby, Denmark (Dated: October 25, 2018)We calculate single-particle excitation energies for a series of 33 molecules using fully selfcon-sistent GW, one-shot G W , Hartree-Fock (HF), and hybrid density functional theory (DFT). Allcalculations are performed within the projector augmented wave (PAW) method using a basis setof Wannier functions augmented by numerical atomic orbitals. The GW self-energy is calculatedon the real frequency axis including its full frequency dependence and off-diagonal matrix elements.The mean absolute error of the ionization potential (IP) with respect to experiment is found to be4.4, 2.6, 0.8, 0.4, and 0.5 eV for DFT-PBE, DFT-PBE0, HF, G W [HF], and selfconsistent GW,respectively. This shows that although electronic screening is weak in molecular systems its inclu-sion at the GW level reduces the error in the IP by up to 50% relative to unscreened HF. In generalGW overscreens the HF energies leading to underestimation of the IPs. The best IPs are obtainedfrom one-shot G W calculations based on HF since this reduces the overscreening. Finally, we findthat the inclusion of core-valence exchange is important and can affect the excitation energies by asmuch as 1 eV. PACS numbers: 31.15.A-,33.15.Ry,31.15.V-
I. INTRODUCTION
Density functional theory (DFT) with the single-particle Kohn-Sham (KS) scheme is today the mostwidely used approach to the electronic structure problemof real materials in both solid state physics and quantumchemistry. While properties derived from total energies(or rather total energy differences) are accurately pre-dicted by DFT, it is well known that DFT suffers froma band gap problem implying that the single-particle KSeigenvalues cannot in general be interpreted as real quasi-particle (QP) excitation energies. In particular, semilocalexchange-correlation functionals severely underestimatethe fundamental gap of both insulators, semi-conductors,and molecules. The hybrid and screened hybrid functionals,which admix around 25% of the (screened) Fock ex-change with the local DFT exchange, generally improvethe description of band gaps in bulk semi-conductorsand insulators . However, the orbital energies obtainedfor finite systems using such functionals still underes-timate the fundamental gap, I p − E a , (the differencebetween ionization potential and electron affinity) byup to several electron volts. In fact, for molecules thepure Hartree-Fock (HF) eigenvalues are usually closerto the true electron addition/removal energies than arethe hybrid DFT eigenvalues. This is because HF isself-interaction free and because screening of the ex-change interaction is a relatively weak effect in molecularsystems. . On the other hand, in extended systemsthe effect of self-interaction is less important and thelong range Coulomb interaction becomes short rangeddue to dynamical screening. As a consequence HF breaksdown in extended systems leading to dramatically over-estimated band gaps and a qualitatively incorret descrip-tion of metals. The many-body GW approximation of Hedin has been widely and succesfully used to calculateQP band structures in metals, semi-conductors, andinsulators. The GW approximation can be viewedas HF with a dynamically screened Coulomb interac-tion. The fact that the screening is determined by thesystem itself instead of being fixed a priori as in thescreened hybrid schemes, suggests that the GW methodshould be applicable to a broad class of systems rang-ing from metals with strong screening to molecules withweak screening. With the entry of nanoscience the use ofGW has been extended to low-dimensional systems andnanostructures and more recently even nonequilib-rium phenomena like quantum transport . In viewof this trend it is important to establish the perfor-mance of the GW approximation for other systems thanthe crystalline solids. In this work we present first-principles benchmark GW calculations for a series ofsmall molecules. In a closely related study we comparedGW and Hartree-Fock to exact diagonalization results forsemi-empirical PPP models of conjugated molecules .The main conclusions from the two studies regarding thequalities of the GW approximation in molecular systems,are very consistent.Most GW calculations to date rely on one or severalapproximations of more technical character. These in-clude the plasmon pole approximation, the linearizedQP equation, neglect of off-diagonal matrix elementsin the GW self-energy, analytic continuations from theimaginary to the real frequency axis, neglect of corestates contributions to the self-energy, neglect of self-consistency. The range of validity of these approxima-tions has been explored for solid state systems by a num-ber of authors , however, much less is known abouttheir applicability to molecular systems . Our imple-mentation of the GW method avoids all of these tech-nical approximations allowing for a direct and unbiasedassessment of the GW approximation itself.Here we report on single-shot G W and fully self-consistent GW calculations of QP energies for a set of33 molecules. The calculated IPs are compared with ex-perimental values as well as single-particle eigenvaluesobtained from Hartree-Fock and DFT-PBE/PBE0 theo-ries. As additional benchmarks we compare to second-order M¨oller-Plesset (MP2) and DFT-PBE total energydifferences between the neutral and cation species. Spe-cial attention is paid to the effect of selfconsistency in theGW self-energy and the role of the initial Green func-tion, G , used in one-shot G W calculations. The useof PAW rather than pseudopotentials facilitate the in-clusion of core-valence exchange, which we find can con-tribute significantly to the HF and GW energies. Ourresults show that the GW approximation yields accuratesingle-particle excitation energies for small molecules im-proving both hybrid DFT and full Hartree-Fock results.The paper is organized as follows. In Sec. II we de-scribe the theoretical and numerical details behind theGW calculations, including the augmented Wannier func-tion basis set, the self-consistent solution of the Dysonequation, and the evaluation of valence-core exchangewithin PAW. In Sec. III we discuss and compare theresults of G W , GW, HF, PBE0, and PBE calculations.We analyze the role of dynamical screening, and discussthe effect of self-consistency in the GW self-energy. Weconclude in Sec. IV. II. METHODA. Augmented Wannier function basis
For the GW calculations we apply a basis set consist-ing of projected Wannier functions (PWF) augmentedby numerical atomic orbitals (NAO). The PWFs, φ i , areobtained by maximizing their projections onto a set oftarget NAOs, Φ Alm , subject to the condition that theyspan the set of occupied eigenstates, ψ n . Thus we maxi-mize the functionalΩ = X i X A,l,m |h φ i | Φ Alm i| (1)subject to the condition span { φ i } ⊇ span { ψ n } occ asdescribed in Ref. 44. The target NAOs are given byΦ Alm ( r ) = ζ Al ( r ) Y lm ( r ) where ζ Al is a modified Gaus-sian which vanish outside a specified cut-off radius, and Y lm are the spherical harmonics corresponding to the va-lence of atom A . The number of PWFs equals the num-ber of target NAOs. For example we obtain one PWFfor H ( l max = 0), and four PWFs for C ( l max = 1). ThePWFs mimick the target atomic orbitals but in additionthey allow for an exact representation of all the occupiedmolecular eigenstates. The latter are obtained from anaccurate real-space PAW-PBE calculation .The PWFs obtained in this way provide an exact rep-resentation of the occupied PBE eigenstates. However, this does not suffice for GW calculations because the po-larizability, P , and the screened interaction, W , do notlive in this subspace. Hence we augment the PWFs byadditional NAOs including so-called polarization func-tions which have l = l max + 1 and/or extra radial func-tions (zeta functions) for the valence atomic orbitals. Formore details on the definition of polarization- and higherzeta functions we refer to Ref. 46. To give an exam-ple, a double-zeta-polarized (DZP) basis consists of thePWFs augmented by one set of NAOs corresponding to l = 0 , ..., l max and one set of polarization orbitals. Notethat the notation, SZ, SZP, DZ, DZP, etc., is normallyused for pure NAO basis sets, but here we use it to de-note our augmented Wannier basis set. We find that theaugmented Wannier basis is significantly better for HFand GW calculations than the corresponding pure NAObasis.The GW and HF calculations presented in Sec. IIIwere performed using a DZP augmented Wannier basis.This gives a total of 5 basis functions per H, Li, and Na,and 13 basis functions for all other chemical elements con-sidered. In Sec. III C we discuss convergence of the GWcalculations with respect to the size of the augmentedWannier basis. B. GW calculations
The HF and GW calculations for isolated moleculesare performed using a Green function code developedfor quantum transport. In principle, this scheme is de-signed for a molecule connected to two electrodes withdifferent chemical potentials µ L and µ R . However, thecase of an isolated molecule can be treated as a specialcase by setting µ L = µ R = µ and modelling the cou-pling to electrodes by a constant imaginary self-energy,Σ L/R = iη . The chemical potential µ is chosen to lie inthe HOMO-LUMO gap of the molecule and the size of η , which provides an artifical broadening of the discretelevels, is reduced until the results have converged. In thislimit of small η the result of the GW calculation becomesindependent of the precise position of µ inside the gap.In Ref. 47 the GW-transport scheme was described forthe case of an orthogonal basis set and for a truncated,two-index Coulomb interaction. Below we generalize therelevant equations to the case of a non-orthogonal basisand a full four-index Coulomb interaction. Some cen-tral results of many-body perturbation theory in a non-orthogonal basis can be found in Ref. 48.The central object is the retarded Green function, G r , G r ( ε ) = [( ε + iη ) S − H KS + v xc − ∆ v H − Σ r xc [ G ]( ε )] − (2)In this equation all quantities are matrices in the aug-mented Wannier basis, e.g. H KS ,ij = h φ i | ˆ H KS | φ j i isthe KS Hamiltonian matrix and S ij = h φ i | φ j i is anoverlap matrix. The term ∆ v H represents the changein the Hartree potential relative to the DFT Hartreepotential already contained in H KS , see Appendix A.The local xc-potential, v xc , is subtracted to avoid dou-ble counting when adding the many-body self-energy,Σ xc [ G ]. As indicated, the latter depends on the Greenfunction and therefore Eq. (2) must in principle be solvedself-consistently in conjuction with the self-energy.In the present study Σ xc is either the bare exchangepotential or the GW self-energy. To be consistent withthe code used for the calculations, we present the equa-tions for the GW self-energy on the so-called Keldyshcontour. However, under the equilibrium conditions con-sidered here the Keldysh formalism is equivalent to theordinary time-ordered formalism.The GW self-energy is defined byΣ GW ij ( τ, τ ′ ) = i X kl G kl ( τ, τ ′ + ) W ik,jl ( τ, τ ′ ) , (3)where τ and τ ′ are times on the Keldysh contour, C .The dynamically screened Coulomb interaction obeys theDyson-like equation W ij,kl ( τ, τ ′ ) = V ij,kl δ C ( τ, τ ′ )+ X pqrs Z C d τ V ij,pq P pq,rs ( τ, τ ) W rs,kl ( τ , τ ′ ) , (4)and the polarization bubble is given by P ij,kl ( τ, τ ′ ) = − iG ik ( τ, τ ′ ) G lj ( τ ′ , τ ) . (5)In the limit of vanishing polarization, P = 0, W reducesto the bare Coulomb interaction V ij,kl = Z Z d r d r ′ | r − r ′ | φ i ( r ) φ ∗ j ( r ) φ ∗ k ( r ′ ) φ l ( r ′ ) (6)and the GW self-energy reduces to the exchange potentialof HF theory.From the above equations for the contour-orderedquantities, the corresponding real time components, i.e.the retarded, advanced, lesser, and greater components,can be obtained from standard conversion rules . Forcompleteness we give the expressions for the real timecomponents of the GW equations in Appendix A.The time/energy dependence of the dynamical quan-tities G , W , P , and Σ, is represented on a uniform grid.We switch between time and energy domains using theFast Fourier Transform in order to avoid time consumingconvolutions. A typical energy grid used for the GW cal-culations in this work ranges from -150 to 150 eV witha grid spacing of 0.02 eV. The code is parallelized overbasis functions and energy grid points. We use a Pulaymixing scheme for updating the Green function G r wheniterating Eq. (2) to self-consistency as described in Ref.47.We stress that no approximation apart from the fi-nite basis set is made in our implementation of the GWapproximation. In particular the frequency dependenceis treated exactly and analytic continuations from theimaginary axis are avoided since we work directly on thereal frequency/time axis. The price we pay for this is thelarge size of the energy grid. C. Spectral function
The single-particle excitation spectrum is contained inthe spectral function A ( ε ) = i ( G r ( ε ) − [ G r ( ε )] † ) . (7)For a molecule A ( ε ) shows peaks at the QP energies ε n = E n ( N + 1) − E ( N ) and ε n = E ( N ) − E n ( N − E n ( N ) denotes the energy of the n thexcited state of the system with N electrons and N refersto the neutral state.When the Green function is evaluated in a non-orthogonal basis, like the augmented Wannier basis usedhere, the projected density of states for orbital φ i be-comes D i ( ε ) = [ SA ( ε ) S ] ii / πS ii , (8)where matrix multiplication is implied. Correspond-ingly, the total density of states, or quasiparticle spec-trum, is given by D ( ε ) = Tr( A ( ε ) S ) / π. (9) D. Calculating Coulomb matrix elements
The calculation of all of the Coulomb matrix elements, V ij,kl , is prohibitively costly for larger basis sets. For-tunately the matrix is to a large degree dominated bynegligible elements. To systematically define the mostsignificant Coulomb elements, we use the product basistechnique of Aryasetiawan and Gunnarsson . In thisapproach, the pair orbital overlap matrix S ij,kl = h n ij | n kl i , (10)where n ij ( r ) = φ ∗ i ( r ) φ j ( r ) is used to screen for the sig-nificant elements of V .The eigenvectors of the overlap matrix Eq. (10) rep-resents a set of “optimized pair orbitals” and the eigen-values their norm. Optimized pair orbitals with insignif-icant norm must also yield a reduced contribution to theCoulomb matrix, and are omitted in the calculation of V . We limit the basis for V to optimized pair orbitalswith a norm larger than 10 − a − . This gives a signifi-cant reduction in the number of Coulomb elements thatneeds to be evaluated, and it reduces the matrix size of P ( ε ) and W ( ε ) correspondingly, see Appendix A.The evaluation of the double integral in Eq. (6) isefficiently performed in real space by solving a Poissonequation using multigrid techniques . E. Valence-core exchange
All inputs to the GW/HF calculations, i.e. the selfcon-sistent Kohn-Sham Hamiltonian, H KS , the xc potential v xc , the Coulomb matrix elements, V ij,kl , are calculatedusing the real-space PAW code GPAW .In GPAW, the core electrons (which are treated scalar-relativistically) are frozen into the orbitals of the freeatoms, and the Kohn-Sham equations are solved for thevalence states only. Unlike pseudo potential schemes,these valence states are subject to the full potential ofthe nuclei and core electrons. This is achieved by a parti-tioning scheme, where quantities are divided into pseudocomponents augmented by atomic corrections. The oper-ators obtained from GPAW are thus full-potential quan-tities, and the wave functions from which the Wannierbasis functions are constructed correspond to the all-electron valence states. Ref. 53 describes how the all-electron Coulomb elements can be determined within thePAW formalism.Since both core and all-electron valence states areavailable in the PAW method, we can evaluate the contri-bution to the valence exchange self-energy coming fromthe core electrons. As the density matrix is simply theidentity matrix in the subspace of atomic core states, thisvalence-core exchange readsΣ core x,ij = − core X k V ik,jk , (11)where i, j represent valence basis functions. We limitthe inclusion of valence-core interactions to the exchangepotential, neglecting it in the correlation. This is reason-able, because the polarization bubble, P , involving coreand valence states will be small due to the large energydifference and small spatial overlap of the valence andcore states. This procedure was used and validated forsolids in Ref. 42. We find that the elements of Σ core x,ij can be significant – on average 1.2 eV for the HOMO –and are larger (more negative) for the more bound or-bitals which have larger overlap with the core states. Ingeneral, the effect on the HOMO-LUMO gap is to en-large it, on average by 0.4 eV because the more boundHOMO level is pushed further down than the less boundLUMO state. In the case of solids, the role of valence-core interaction has been investigated by a number ofauthors . Here the effect on the QP band gapseems to be smaller than what we find for the molec-ular gaps. We note that most GW calculations rely onpseudopotential schemes where these valence-core inter-actions are not accessible. In such codes, the xc contri-bution from the core electrons are sometimes estimatedby Σ corexc ≈ v xc [ n ] − v xc [ n val ] where n val is the valenceelectron density, but as the local xc potential is a non-linear functional of the density, this procedure is not welljustified. Instead we subtract the xc potential of the fullelecton density n , and add explicitly the exact exchangecore contribution. III. RESULTS
In Fig. 1 we compare the calculated HOMO en-ergies with experimental ionization potentials for the33 molecules listed in Table I. The geometries of themolecules, which all belong to the G2 test set, are takenfrom Ref. 56. The different HOMO energies correspondto: DFT-PBE and DFT-PBE0 eigenvalues, Hartree-Fock eigenvalues, and fully selfconsistent GW. The GWenergies are obtained from the peaks in the correspond-ing density of states Eq. (9) extrapolated to η = 0 ( η gives an artificial broadening of the delta peaks). FIG. 1: (Color online) Calculated negative HOMO energyversus experimental ionization potential. Both PBE andPBE0 systematically understimates the ionization energy dueto self-interaction errors while HF overestimates it slightly.The dynamical screening from the GW correlation lowers theHF energies bringing them closer to the experimental values.Numerical values are listed in Table I.
We stress the different meaning of fully selfconsistentGW and the recently introduced method of quasiparti-cle selfconsistent GW : In fully selfconsistent GW theGreen function obtained from Dyson’s equation Eq. (2)with Σ xc [ G ] = Σ GW [ G ] is used to calculate the Σ GW ofthe next iteration. In QP-selfconsistent GW, Σ GW is al-ways evaluated using a non-interacting Green functionand the self-consistency is obtained when the differencebetween the non-interacting GF and the interacting GF,is minimal.Fig. 1 clearly shows that both the PBE and PBE0eigenvalues of the HOMO severely underestimates theionization potential. The average deviation from the ex-perimental values are 4.35 eV and 2.55 eV, respectively.The overestimation of the single-particle eigenvalues ofoccupied states is a well known problem of DFT and canbe ascribed to the insufficient cancellation of the self-interaction in the Hartree potential. Part of this self-interaction is removed in PBE0. However, the fact thatthe HF results are significantly closer to experiments in-dicates that the 25% Fock exchange included in the PBE0is not sufficient to cure the erroneous description of (oc-cupied) molecular orbitals. On the other hand PBE0gives good results for band gaps in semi-conductors andinsulators where in contrast full Hartree-Fock does notperform well.
We conclude that the amount of Fockexchange to be used in the hybrid functionals to achievegood quasiparticle energies is highly system dependent.A similar problem is encountered with self-interactioncorrected exchange-correlation functionals. FIG. 2: The deviation of the calculated HOMO energy fromthe experimental ionization potential in GW and HF, respec-tively. The vertical displacement of points from the line x = y gives the difference between the GW and HF energies and rep-resents the effect of screening. Notice that the GW correctionis always negative (corresponding to higher HOMO energy)and that it generally overcorrects the HF energies. Also no-tice that the GW correction is larger for molecules where HFpresents the largest overestimation of the ionization potential. As can be seen from Fig. 1, GW performs better thanHartree-Fock for the HOMO energy yielding a mean ab-solute error with respect to experiments of 0.5 eV com-pared to 0.81 eV with Hartree-Fock. As expected the dif-ference between HF and GW is not large on an absolutescale (around 1 eV on average, see Table II) illustratingthe fact that screening is weak in small molecules. On arelative scale selfconsistent GW improves the agreementwith experiments by almost 30% as compared to HF.To gain more insight into the influence of screening onthe orbital energies, we compare in Fig. 2 the deviationof the HF and GW energies from IP exp . The GW self-energy can be split into the bare exchange potential andan energy-dependent correlation partΣ GW ( r , r ′ ; ε ) = v x ( r , r ′ ) + Σ corr ( r , r ′ ; ε ) (12) Accordingly the quasiparticle energy can be written asthe bare HF energy and a correction due to the energy-dependent part of the GW self-energy (the dynamicalscreening term) ε QP n = ε HF n + ∆ GW n . (13)In Fig. 2 the line y = x corresponds to ∆ GW n = 0, andthe vertical displacement from the line thus representsthe effect of screening on the calculated HOMO energy.We first notice that the effect of screening is to shift theHOMO level upwards in energy, i.e. to reduce the ion-ization potential. This can be understood by recallingthat the Hartree-Fock eigenvalue represents the energycost of removing an electron from the HOMO when or-bital relaxations in the final state are neglected (Koop-mans’ theorem ). In Ref. 37 we showed, on the basisof GW and exact calculations for semi-empirical modelsof conjugated molecules, that ∆ GW n mainly describes theorbital relaxations in the final state and to a lesser extentaccounts for the correlation energy of the initial and finalstates. This explains the negative sign of ∆ GW n becausethe inclusion of orbital relaxation in the final state lowersthe energy cost of removing an electron. We note thatthis is different from the situation in extended, periodicsystems where orbital relaxations vanish and the main ef-fect of the GW self-energy is to account for correlationsin the initial and final states.In Table I we list the calculated HOMO energy for eachof the 33 molecules. In addition to selfconsistent GW wehave performed one-shot G W calculations based on theHF and PBE Green’s function, respectively. The bestagreement with experiment is obtained for G W [HF].This is because the relatively large Hartree-Fock HOMO-LUMO gap reduces the (over-)screening described by theresulting GW self-energy. There are not many GW cal-culations for molecules available in the litterature. BelowTable I we list the few we have found. As can be seenthey all compare quite well with our results given thedifferences in the implementation of the GW approxima-tion.For comparison we have included the HOMO energypredicted by second order Møller-Plesset theory (MP2)[taken from Ref. 59] with a Gaussian 6-311G ∗∗ basisset. These are generally very close to our calculated HFvalues, with a tendency to lower energies which worsensthe agreement with experiment slightly as compared toHF.We have also calculated the DFT-PBE total energy dif-ference between the neutral and cation species, E ( N ) − E ( N − TABLE I: Experimental ionization potential (first column) and HOMO energy calculated using different approximations forexchange and correlation. “X-eig” refers to a single-particle eigenvalue while “X-tot” refers to a total energy difference, E ( N ) − E ( N − W (PBE) energies have been obtained from the QP equation while the GW and G W energies areobtained from the DOS in Eq. (9). Last row shows the mean absolute error (MAE) with respect to experiments. All energiesare in eV.Molecule Expt. ( a ) PBE-eig PBE0-eig HF-eig GW G W (HF) G W (PBE)-QP MP2 ( a ) PBE-totLiH 7.90 4.34 5.81 8.14 8.0 ( b ) ( b ) H H ( c ) Cl 11.29 7.08 8.80 11.68 11.0 11.4 11.1 11.90 11.10CH OH 10.96 6.31 8.49 12.14 10.7 10.8 10.5 12.16 10.72CH SH 9.44 5.60 7.09 9.50 8.8 9.0 8.4 9.73 9.29Cl O CO 10.88 6.28 8.37 11.93 10.4 10.5 10.6 11.97 10.80HCN 13.61 9.05 10.67 13.19 12.7 13.2 12.4 13.33 13.67HF 16.12 9.61 12.47 17.74 16.0 15.6 15.7 17.35 16.27H O 12.62 7.24 9.59 13.88 12.3 12.1 11.9 ( d ) H ( e ) H ( a ) From Ref. 59. The MP2 calculations use a Gaussian 6-311G ∗∗ basis set. ( b ) To be compared with the GW value 7.85 and the G W (HF) value 8.19 reported in Ref. 24. ( c ) To be compared with the G W (LDA) value 14.3 reported in Ref. 21. ( d ) To be compared with the G W (LDA) value 11.94 reported in Ref. 22. ( e ) To be compared with the G W (LDA) values 12.7 and 12.66 reported in Refs. 21 and 22, respectively. In Table II we provide an overview of the compara-tive performance of the different methods. Shown is themean average deviation between the IPs calculated withthe different methods as well as the experimental values.Note that the numbers in the experiment row/column are the same as those listed in the last row of Table II.
TABLE II: Mean absolute deviation between the IPs of the 33 molecules calculated with the different methods and experiment.The mean absolute deviation with respect to experiment coincide with the last row in Table IMethod Expt. (a)
PBE-eig PBE0-eig HF-eig GW G W [HF] MP2 (a) PBE-totExpt. 0.00 4.35 2.55 0.81 0.5 0.4 0.82 0.24PBE 4.35 0.00 1.79 4.90 3.9 4.1 4.99 4.27PBE0 2.55 1.79 0.00 3.11 2.1 2.3 3.20 2.48HF 0.81 4.90 3.11 0.00 1.0 0.8 0.17 0.80GW 0.5 3.9 2.1 1.0 0.00 0.3 1.1 0.4G W [HF] 0.4 4.1 2.3 0.8 0.3 0.00 0.9 0.3MP2 0.82 4.99 3.20 0.17 1.1 0.9 0.00 0.84PBE-tot 0.24 4.27 2.48 0.80 0.4 0.3 0.84 0.00 ( a ) Data taken from Ref. 59. -50 -40 -30 -20 -10 0 ε (eV) -1001020 E ne r g y ( e V ) DOS (GW)DOS (HF) ε nQP ε nHF Re Σ ε − ε nHF -13 -12 -11 -10
FIG. 3: (Color online) Density of states for the NH moleculecalculated in HF and GW, respectively. Arrows mark the levelcorresponding to the HOMO in the two calculations. The in-tersection between the line y = ε − ε HF n and the real part of h ψ | Σ corr ( ε ) | ψ i (green curve) determines the posi-tion of the GW level. A. Linearized quasiparticle equation
In the conventional GW method the full Green func-tion of Eq. (2) is not calculated. Rather one obtains thequasiparticle energies from the quasiparticle equation ε QP n = ε n + Z n h ψ n | Σ GW ( ε n ) − v xc | ψ n i . (14)where ψ n and ε n are eigenstates and eigenvalues of anapproximate single-particle Hamiltonian (often the LDAHamiltonian), and Z n = h − ∂ h ψ DFT n | Σ GW ( ε ) | ψ n i ∂ε (cid:12)(cid:12)(cid:12) ε n i − . (15)Moreover the GW self-energy is evaluated non-selfconsistently from the single-particle Green function,i.e. Σ GW = iG W [ G ], with G ( z ) = ( z − H ) − .The quasiparticle equation (14) relies on the assump-tion that off-diagonal matrix elements, h ψ n | Σ GW ( ε n ) − v xc | ψ m i , can be neglected, and that the frequency de-pendence of Σ GW can be approximated by its first or-der Taylor expansion in a sufficiently large neighbor-hood of ε n . We have found that these two assumptionsare indeed fullfilled for the molecular systems studiedhere. More precisely, for the GW and G W (HF) self-energies, the QP energies obtained from Eq. (14) arealways very close to the peaks in the density of statesEq. (9). We emphasize that this result could well be re-lated to the rather large level spacing of small molecules,and may not hold for extended systems. An example ispresented in Fig. 3 which shows the full HF and GWdensity of states for NH together with the real part of h ψ | Σ corr ( ε ) | ψ i . As explained in the follow-ing section this is not quite the case for the G W (PBE)calculations. B. G -dependence As stated in the previous section the GW andG W (HF) energies can be obtained either from the fullspectral function or from the QP equation. In this case,returning to Table I, we see that G W (HF) yields sys-tematically larger IPs than GW. This is easy to under-stand since G HF describes a larger HOMO-LUMO gapthan G GW , and therefore produces less screening. Whenthe PBE rather than the HF Green function is used toevaluate the GW self-energy, we find that the spectralfunction obtained from Eq. (2) does not resemble asimple discrete spectrum. In fact the peaks are signif-icantly broadened by the imaginary part of Σ GW and itbecomes difficult to assign precise values to the QP ener-gies. Apart from the spectral broadening, the moleculargap is significantly reduced with respect to its value in theGW and G W (HF) calculations. Both of these effectsare due to the very small HOMO-LUMO gap describedby G PBE which leads to severe overscreening and QP life-time reductions. A similar effect was observed by Ku andEguiluz in their comparison of GW and G W (LDA) forSi and Ge crystals .The problems encountered when attempting to solvethe Dyson equation (2) using the G W (PBE) self-energyoccur due to the large mismatch between ε PBE n and ε QP n .On the other hand, in the QP equation, the GW self-energy is evaluated at ε n rather than ε QP n . As a conse-quence the unphysical broadening and overscreening isavoided and a well defined QP energy can be obtained(last column in Table I).To summarize, G can have a very large effect on theQP spectrum when the latter is obtained via the Dysonequation (2). In particular, the use of a G with a toonarrow energy gap (as e.g. the G PBE ) can lead to un-physical overscreening and spectral broadening. Whenthe QP levels are obtained from the QP equation, the G -dependence is less pronounced because Σ GW [ G ] isevaluated at ε n which is consistent with G .The self-consistent GW spectrum is independent of thechoice of G , but the number of iterations required toreach self-consistency is less when based on G HF . C. Basis set convergence
In Figs. 4 and 5 we show the energy of the three high-est occupied molecular orbitals of H O and CO obtainedfrom selfconsistent GW using various sizes of the aug-mented Wannier basis. Clearly, the polarization func-tions have relatively little influence on the QP energieswhile the first set of additional zeta functions reduce theQP energies by up to 0.5 eV. The differences betweenDZP and TZDP are less than 0.15 eV for all the levelswhich justifies the use of DZP basis.We have also compared the eigenvalues obtained fromselfconsistent HF calculations using the DZP augmentedWannier basis to accurate HF calculations performedwith the real-space code GPAW . Here we obtain aMAE of 0.09 eV for the energy of the HOMO level of the33 molecules. IV. CONCLUSIONS
As the range of systems to which the GW method isbeing applied continues to expand it becomes importantto establish its performance for other systems than thesolids. In this work we have discussed benchmark GWcalculations for molecular systems.The GW calculations were performed using a novelscheme based on the PAW method and a basis set consist-ing of Wannier functions augmented by numerical atomicorbitals. We found that a basis corresponding to double-zeta with polarization functions was sufficient to obtainGW energies converged to within 0.1 eV (compared totriple-zeta with double polarization functions). The GWself-energy was calculated on the real frequency axis in-cluding its full frequency dependence and off-diagonalelements. We thereby avoid all of the commonly usedapproximations, such as the the plasmon pole approx-imation, the linearized quasiparticle equation and ana-
FIG. 4: (Color online) Convergence of the three highest oc-cupied levels of H O obtained from GW calculations withdifferent sizes of the augmented Wannier function basis. SZdenotes the Wannier function basis, while e.g. DZDP denotesthe Wannier basis augmented by one extra radial function pervalence state and two sets of polarization functions.FIG. 5: (Color online) Same as Fig. 4 but for CO. lytical continuations from imaginary to real frequencies,and thus obtain a direct and unbiased assessment of theGW approximation itself. We found that the inclusionof valence-core exchange interactions, as facilitated bythe PAW method, is important and affect the HF/GWHOMO levels by -1.2 eV on average.The position of the HOMO for a series of 33 moleculeswas calculated using fully selfconsistent GW, single-shotG W , Hartree-Fock, DFT-PBE0, and DFT-PBE. BothPBE and PBE0 eigenvalues grossly overestimate theHOMO energy with a mean absolute error (MAE) withrespect to the experimental ionization potentials (IP) of4.4 and 2.5 eV, respectively. Hartree-Fock underesti-mates the HOMO energy but improves the agreementwith experiments yielding a MAE of 0.8 eV. GW andG W overcorrects the Hartree-Fock levels slightly lead-ing to a small overestimation of the HOMO energy witha MAE relative to experiments of 0.4-0.5 eV. This showsthat although screening is a weak effect in molecular sys-tems its inclusion at the GW level improves the electronremoval energies by 30-50% relative to the unscreenedHartree-Fock. The best IPs were obtained from one-shotG W calculations starting from the HF Green’s func-tion where the overscreening is least severe. Very similarconclusions were reached by comparing GW, G W andHF to exact diagonalization for conjugated molecules de-scribed by the semi-empirical PPP model. V. ACKNOWLEDGMENTS
We thank Mikkel Strange for useful discussions andassistance with the projected Wannier functions. We ac-knowledge support from the Danish Center for ScientificComputing and The Lundbeck Foundation’s Center forAtomic-scale Materials Design (CAMD).
Appendix A: The GW self-energy
Let U denote the rotation matrix that diagonalizes thepair orbital overlap S ij,kl = h n ij | n kl i , i.e. U † SU = σI .The columns of U are truncated to those which havecorresponding eigenvalues σ q < − a − . We then onlycalculate the reduced number of Coulomb elements V qq ′ = h n q | | r − r ′ | | n q ′ i , (A1)where n q ( r ) are the optimized pair orbitals n q ( r ) = X ij n ij ( r ) U ij,q / √ σ q , (A2)which are mutually orthonormal, i.e. h n q | n q ′ i = δ qq ′ .Determining the GW self-energy proceeds by calculat-ing first the full polarization matrix in the time domain P
Solid StatePhysics , edited by H. Ehrenreich and F. Saepen (AcademicPress, New York 2000), Vol. 54, 1. G. Onida, L. Reining, and A. Rubio, Rev. Mod. Phys. ,601 (2002) J. C. Grossman, M. Rohlfing, L. Mitas, S. G. Louie, andM. L. Cohen, Phys. Rev. Lett. , 472 (2001). P. H. Hahn, W. G. Schmidt, and F. Bechstedt Phys. Rev.B , 245425 (2005) T. Niehaus, M. Rohlfing, F. D. Sala, A. Di Carlo, and T.Frauenhaim, Phys. Rev. A , 022508 (2005) A. Stan, N. E. Dahlen, and R. van Leeuwen, Europhys.Lett. , 298 (2006) C. D. Spataru, S. Ismail-Beigi, X. L. Benedict, and S. G.Louie, Phys. Rev. Lett. , 077402 (2004) P. E. Trevisanutto, C. Giorgetti, L. Reining, M. Ladisa,and V. Olevano, Phys. Rev. Lett. , 226405 (2008) M. L. Tiago and J. R. Chelikowsky, Phys. Rev. B ,205334 (2006) J. B. Neaton, M. S. Hybertsen, and S. G. Louie Phys. Rev.Lett. , 216405 (2006) K. S. Thygesen and A. Rubio, Phys. Rev. Lett. ,046802 (2009) C. Freysoldt and P. Rinke and M. Scheffler, Phys. Rev.Lett. , 056803 (2009) J. M. Garcia-Lastra, C. Rostgaard, A. Rubio, and K. S.Thygesen Phys. Rev. B , 245427 (2009) K. S. Thygesen and A. Rubio, J. Chem. Phys. , 091101(2007). K. S. Thygesen, Phys. Rev. Lett. , 166804 (2008) P. My¨ohanen and A. Stan and G. Stefanucci and R. vanLeeuwen, Euro. Phys. Lett. , 67001 (2008) C. D. Spataru and M. S. Hybertsen and S. G. Louie andA. J. Millis, Phys. Rev. B , 155110 (2009) M. P. von Friesen and C. Verdozzi and C.-O. Almbladh,Phys. Rev. Lett. , 176404 (2009) K. Kaasbjerg and K. S. Thygesen (submitted) B. Holm and U. von Barth Phys. Rev. B , 2108 (1998) M. Usuda, N. Hamada, T. Kotani, and M. van Schifl-gaarde, Phys. Rev. B , 125101 (2002) W. Ku and A. G. Eguiluz, Phys. Rev. Lett. , 126401(2002) M. van Schiflgaarde, T. Kotani, and S. V. Faleev, Phys.Rev. B , 245125 (2006) M. Shishkin and G. Kresse, Phys. Rev. B , 035101 (2006) P. Rinke, A. Qteish, J. Neugebauer, C. Freysoldt, and M.Scheffler New J. Phys.
126 (2005) X. Qian, J. Li, L. Qi, C.-Z. Wang, T.-L. Chan, Y.-X. Yao,K.-M. Ho, and S. Yip, Phys. Rev. B , 245112 (2008) J. J. Mortensen, L. B. Hansen, K. W. Jacobsen, Phys. Rev.B , 035109 (2005) A. H. Larsen, M. Vanin, J. J. Mortensen, K. S. Thygesen,and K. W. Jacobsen Phys. Rev. B , 195112 (2009) K. S. Thygesen and A. Rubio, Phys. Rev. B , 115333(2008). K. S. Thygesen, Phys. Rev. B , 035309 (2006) H. Haug and A. -P. Jauho,
Quantum Kinetics in Transportand Optics of Semiconductors , Springer (1998) R.van Leeuwen, N.E.Dahlen, G.Stefanucci, C.-O.Almbladh, and U.von Barth, Lectures Notes inPhysics vol. 706 (Springer) A. Stan, N. E. Dahlen, and R. van Leeuwen, J Chem. Phys. , 114105 (2009) F. Aryasetiawan and O. Gunnarsson, Phys. Rev. B ,16214 (1994) M. Walter, H. H¨akkinen, L. Lehtovaara, M. Puska, J.Enkovaara, C. Rostgaard, and J. J. Mortensen, J Chem.Phys. , 244101 (2009) P. E. Bl¨ochl, Phys. Rev. B , 17953 (1994) E. Engel, Phys. Rev. B , 161205(R) (2009) L. A. Curtiss, K. Raghavachari, P. Redfern, and J. A.Pople, J. Chem. Phys. , 1063 (1997). J. P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. , 3865 (1996) M. van Schilfgaarde, T. Kotani, and S. Faleev, Phys. Rev.Lett. , 226402 (2006)59