Benchmarking GW against exact diagonalization for semi-empirical models
aa r X i v : . [ c ond - m a t . m e s - h a ll ] D ec Benchmarking GW against exact diagonalization for semi-empirical models
K. Kaasbjerg and K. S. Thygesen Nano-Science Center, Niels Bohr Institute, University of CopenhagenUniversitetsparken 5, DK-2100 Copenhagen, Denmark Center for Atomic-scale Materials Design (CAMD),Department of Physics, Technical University of Denmark, DK - 2800 Kgs. Lyngby, Denmark (Dated: November 3, 2018)We calculate groundstate total energies and single-particle excitation energies of seven pi con-jugated molecules described with the semi-empirical Pariser-Parr-Pople (PPP) model using self-consistent many-body perturbation theory at the GW level and exact diagonalization. For the totalenergies GW captures around 65% of the groundstate correlation energy. The lowest lying exci-tations are overscreened by GW leading to an underestimation of electron affinities and ionizationpotentials by approximately 0.15 eV corresponding to 2.5%. One-shot G W calculations startingfrom Hartree-Fock reduce the screening and improve the low-lying excitation energies. The effect ofthe GW self-energy on the molecular excitation energies is shown to be similar to the inclusion offinal state relaxations in Hartree-Fock theory. We discuss the break down of the GW approximationin systems with short range interactions (Hubbard models) where correlation effects dominate overscreening/relaxation effects. Finally we illustrate the important role of the derivative discontinuityof the true exchange-correlation functional by computing the exact Kohn-Sham levels of benzene. PACS numbers: 31.15.bu,33.15.Ry,31.15.V-
I. INTRODUCTION
For more than two decades the many-body GWapproximation of Hedin has been the state of theart for calculating band structures of metals, semi-conductors, and insulators . With the entry ofnanoscience the use of the GW method has been ex-tended to low-dimensional systems such as molecules,carbon nanotubes, graphene and molecule-surface in-terfaces . In these systems the interplaybetween quantum confinement (in one or more di-mensions) and electronic correlation effects leads tonovel phenomena like the renormalization of molecu-lar electronic levels at surfaces by dynamical polariza-tion in the substrate . Very recently, the non-equilibrium version of the GW approximation has beenapplied to quantum transport and dynamics in molecularjunctions where dynamic correlationsseems to be particularly important.As the range of systems to which the GW approxima-tion is being applied continues to expand, critical inves-tigations of the performance of GW for other systemsthan the crystalline solids become important. Here wereport on benchmark GW calculations for π -conjugatedmolecules based on the semi-empirical Pariser-Parr-Pople(PPP) model . By comparing with exact results weobtain a direct and unbiased estimate of the quality ofthe GW approximation in molecular systems.Previous benchmark model studies of the GW approx-imation have all focused on Hubbard models with localinteractions with the conclusion that GW workswell for small interaction strengths but fails for largerinteractions strength. The use of GW in systems withlocal interactions is in fact unfortunate because the im-portance of electronic screening, which is the main effect described by GW, is weak in comparison to correlationeffects. In contrast to Hubbard models, the PPP descrip-tion includes long range interactions and its parametershave been fitted to yield realistic excitation energies ofconjugated molecules. It therefore provides a better andmore natural starting point for a study addressing the ac-curacy of GW for real molecules and nanostructures. Wemention that in a related work we have performed first-principles GW calculations for a series of 33 moleculesarriving at very similar conclusions regarding the perfor-mance of GW as those reported here. Ab-initio
GW calculations typically involve a numberof ”technical” approximations such as the plasmon poleapproximation, the neglect of off-diagonal matrix ele-ments in the GW self-energy, or analytic continuations.Moreover they are usually performed non-selfconsistentlyand are subject to basis set errors. In the present workthe GW calculations are carried out fully self-consistentlywithout any further approximations apart from the GWapproximation itself.In this work we calculate total energies and excitationspectra of the seven conjugated molecules listed in Tab. I.The excitation spectrum of a system can be obtainedfrom the spectral function A i ( ε ) = 2 π X n (cid:20) |h Ψ n ( N + 1) | c † i | Ψ ( N ) i| δ ( ε − ε n )+ |h Ψ n ( N − | c i | Ψ ( N ) i| δ ( ε − ε n ) (cid:21) , (1)which has peaks at the excitation energies ε n = E n ( N +1) − E ( N ) and ε n = E ( N ) − E n ( N −
1) correspond-ing to electronic addition and removal energies, respec-tively. Often in the GW literature, excitation energiesare often referred to as quasi-particle (QP) energies. Inthe expressions for the excitation energies E n ( N ) denotesthe energy of the n th excited N -electron state, | Ψ n ( N ) i ,with N referring to the neutral state of the system. Formolecules the first addition and the first removal energy,i.e. n = 0, corresponds to the electron affinity and theionization potential. In Hartree-Fock theory Koopman’stheorem states that the eigenvalues of the Hartree-FockHamiltonian equal the addition/removal energies calcu-lated without orbital relaxations in the charged states,i.e. ε HF n = h c † n Ψ HF0 ( N ) | H | c † n Ψ HF0 ( N ) i− E HF0 ( N ) for a vir-tual orbital n . In particular, the highest occupied molec-ular orbital (HOMO) and the lowest unoccupied orbital(LUMO) represent well defined approximations to theionization potential and electron affinities, respectively .This approximation neglects two important effects. Oneis the relaxation of the single-particle HF orbitals whenan electron is removed from or added to the molecule.The other is the correlation energy which by definition isomitted in HF theory. It is instructive to write the exactQP energies as the sum of the three contributions ε n = ε HF n + ∆ relax + ∆ corr , (2)The relaxation contribution is the correction that followsby calculating the QP energy from self-consistently deter-mined HF energies of the neutral and the charged states N ±
1. The last term ∆ corr is the remaining contribu-tion from the correlation energy. For the addition of anelectron, i.e. an unoccupied orbital, the relaxation andcorrelation contributions are given by∆ relax = E HF n ( N + 1) − E HF0 ( N ) − ε HF n (3)and∆ corr = [ E n ( N + 1) − E HF n ( N + 1)] − [ E ( N ) − E HF0 ( N )] . (4)In extended systems the potential due to a single de-localized electron/hole decreases with the size of the sys-tem. Hence, in such systems there will be no or littlerelaxation of the states due to the addition/removal ofan electron, and the majority of the correction to theQP energy will come from the correlation part ∆ corr .In molecules, nanostructures, molecules at surfaces, anddisordered systems with finite localization lengths, thisis not the case. Here, the introduction of an additionalelectron or hole will lead to a relaxation of the single-particle orbitals corresponding to a screening of the ad-ditional charge. As a consequence, the relaxation cor-rection ∆ relax to the QP energy cannot be neglected insuch systems. In fact, we find that ∆ relax is larger than∆ corr for all the molecules studied here, and that the GWexcitation energies correspond roughly to including only∆ relax in Eq. (2).The paper is outlined as follows. In Sec. II the PPPmodel Hamiltonian for conjugated molecules is intro-duced. In Secs. III A and III B we provide an overview ofthe theory and numerical implementation of the GW andexact calculations, and in Sec. III C we discuss the use of the von Neumann entropy as a measure of correlation.The results for total energies and spectral properties ofthe PPP model are presented in Secs. IV A and IV B, anda comparison is made to short ranged Hubbard modelsin Sec. IV C. In Sec.IV D we calculate the exact Kohn-Sham levels for the benzene molecule and compare to theexact QP levels. The conclusions are given in Sec. V. II. PARISER-PARR-POPLE HAMILTONIAN
The Pariser-Parr-Pople model is an effective π -electron description of conjugated molecules that in-cludes electron-electron interactions explicitly. The PPPHamiltonian is given by H = X i ε i ˆ n i − X h ij i σ t ij c † iσ c jσ + 12 X i = j V ij (ˆ n i − Z i )(ˆ n j − Z j ) + X i U i ˆ n i ↑ ˆ n i ↓ , (5)where c † i ( c i ) creates (annihilates) an electron in the p z orbital on atom i of the molecule, ˆ n i = ˆ n i ↑ + ˆ n i ↓ is thenumber operator, ˆ n iσ = c † iσ c iσ , Z i is the valence (i.e. thenumber of π electrons) of atom i , and h ij i denotes nearestneighbour hopping. The Ohno parametrization is usedfor the long range interactions V ij = 14 . q (28 . / ( U i + U j )) + R ij , (6)where R ij is the inter-atomic distance (in ˚A) and U i is theonsite Coulomb interaction (in eV). For large distancesthe Ohno parametrization recovers the 1 /r behavior ofthe Coulomb interaction while it for small distances rep-resents a screened interaction that interpolates to onsiteCoulomb interaction U i for R ij = 0. The onsite energy ε i , the hopping element t ij and the onsite Coulomb in-teraction U i are treated as fitting parameters. In thepresent work values for these parameters have been takenfrom the literature . Since existing parametershave been optimized to optical excitation spectra, an ex-act agreement with experimental values for the moleculargaps is not to be expected. III. METHODSA. GW approximation
Hedin’s equations provides a formally exact frame-work for the determination of the single-particle Greenfunction in a self-consistent manner. In the GW approx-imation, which follows by neglecting the so called vertexcorrections, the electronic self-energy Σ is given by theproduct of the Green function G and the screened inter-action W , and can be written symbolically asΣ = iGW, (7)where the Green function obeys the usual Dyson equation G = G + G Σ G . The screened interaction W is givenby the bare Coulomb interaction V and the polarizationin the random-phase approximation (RPA) P = − iGG through the Dyson-like equation W = V + V P W. (8)In fully self-consistent GW the set of coupled equationsfor Σ, G , P , and W are solved iteratively until theGreen function has converged. Due to the computa-tional requirement of a fully self-consistent GW scheme, ab-initio GW calculations are usually carried out non-selfconsistently. This approach, which is referred to asG W , starts from an approximate G , typically the non-interacting Kohn-Sham Green function, from which a sin-gle self-energy iteration is carried out to obtain the finalGreen function.
1. Numerical details
The GW calculations have been performed followingthe method described in detail in Ref. 37. Here we givea brief overview of the method for completeness.The retarded and advanced single-particle Green func-tions are given by G r/a ( ε ) = ( ε ± iη − H − V H − Σ GW ( ε )) − (9)where η is a small positive infinitesimal, H containsthe first two terms in Eq. (5), and V H is the Hartreepotential. We represent the Green functions and allother energy-dependent quantities on a uniform grid, − E m , − E m + dε, . . . , E m . The Fast Fourier Transformis used to switch between the energy and time rep-resentations. Since η determines the minimum widthof features in the Green function’s energy dependence,the energy grid spacing should obey dε ≪ η . All re-sults presented here have been converged with respectto η, dε, E m . Typical converged values are (in eV) η =0 . , dε = 0 . , E m = 50.The lesser/greater Green functions are given by G < ( ε ) = − f ( ε − µ )[ G r − G a ] (10) G > ( ε ) = (1 − f ( ε − µ ))[ G r − G a ] (11)where f ( ε − µ ) is the Fermi-Dirac function. The chemicalpotential µ is adjusted to yield the desired number ofelectrons in the system. The formulation in terms ofa fixed chemical potential rather than a fixed particlenumber is reminiscent of the fact that the method hasbeen developed for quantum transport. The one-bodydensity matrix is given by ρ ij = − i Z G
2. Total energy
The total energy can be split into kinetic (and ex-ternal), Hartree, and exchange-correlation energy E = E + E H + E xc . In terms of the Green function we have E + E H = Tr[ H ρ ] + 12 Tr[ V H ρ ] (20)For the exchange-correlation energy we have E xc = 12 i Z Tr[Σ r ( ε ) G < ( ε ) + Σ < ( ε ) G a ( ε )] dε, (21)where Σ is the exchange-correlation self-energy. In thiswork Σ is either the bare exchange, Σ x , yielding the HFapproximation, or the GW self-energy, Σ GW . The expres-sion (21) follows by expressing h ˆ V i in terms of the two-particle Green function, G , and then using the definingequation for the self-energy in terms of G . B. Exact diagonalization
The most direct way to the spectral properties of asystem is via the Lehmann representation of the Greenfunction in Eq. (1). However, since this requires the fullset of eigenstates and eigenvalues of the Hamiltonian,it is of limited practical use and other routes must betaken. The following section gives a brief overview ofthe Lanczos method for iterative diagonalization of largematrices.
1. Calculating the ground state - Lanczos algorithm
In exact diagonalization the given many-body Hamil-tonian is diagonalized directly in the Fock space which isspanned by many-particle states (Slater determinants).Since the dimensionality of the Fock space grows expo-nentially with the number of basis orbitals, symmetriesof the Hamiltonian can help to reduce the dimensionalityconsiderably. For the Pariser-Parr-Pople Hamiltonian inEq. (5) the number of up and down electrons, N ↑ and N ↓ , are good quantum numbers since their correspond-ing operators commute with the Hamiltonian. This im-plies that the exact diagonalization can be carried outin each of the ( N ↑ , N ↓ )-subblocks of the Fock space inde-pendently. The dimensionality of each ( N ↑ , N ↓ )-subblockis given by the number of ways N ↑ spin up electrons and N ↓ spin down electrons can be distributed over L basisorbitals, d ( N ↑ , N ↓ ) = L ! N ↑ !( L − N ↑ )! × L ! N ↓ !( L − N ↓ )! . (22)Very often the ground state is located in the half-filledsubblock, i.e. N ↑ = N ↓ = L/ L is the number ofbasis orbitals. For L = 16 the dimensionality of this sub-block is d = 165636900, implying that storing a vector indouble floating point precision requires ∼ K generatedby repeated applications of H on an arbitrary initial state | φ i , i.e. K = span {| φ i , H | φ i , H | φ i , · · · , H N | φ i} . (23)In the Krylov subspace the extreme eigenvalues of theHamiltonian converge fast with respect to the size N ofthe subspace, thus reducing the full diagonalization toa manageable diagonalization of a N × N matrix, with N ≪ d .In the Lanczos algorithm the Hamiltonian is pro-jected onto a specially constructed orthogonalised Krylovbasis in which the Hamiltonian has a tridiagonal repre-sentation. The basis vectors are generated recursivelyas | φ n +1 i = H | φ n i − a n | φ n i − b n | φ n − i , (24)where the coefficient are given by a n = h φ n | H | φ n ih φ n | φ n i and b n = h φ n | φ n ih φ n − | φ n − i (25) with initial conditions b = 0 and | φ − i = 0. At anypoint during the Lanczos iterations only three Lanczosvectors needs to be kept in memory, which makes the al-gorithm memory efficient. In the basis of the normalizedvectors (the basis vectors above are not normalized) theHamiltonian has the following tridiagonal representation H = a b · · · b a b ...0 b a . . . 0... . . . . . . b N · · · b N a N (26)which can be readily diagonalized with methods for tridi-agonal matrices. In practice the Lanczos iterations arecontinued until the desired eigenvalues have converged toa given tolerance. For the ground state energy E , typi-cal values for N range from a few to ∼
200 depending onthe system size.The ground state resulting from a diagonalization ofthe tridiagonal Hamiltonian in Eq. (26) is provided inthe Lanczos basis, i.e. | Ψ i = P n c n | φ n i . In order to beable to calculate the Green function, its representationin the original many-body basis is required. Since theLanczos vectors are not stored, the Lanczos iterationsmust be repeated (starting from the same initial vector)to obtain the expansion coefficients α i = P n c n h Φ i | φ n i in the original many-body basis {| Φ i i} di =1 .The most time consuming part of the Lanczos algo-rithm is the matrix-vector multiplication H | φ n i . An ef-ficient implementation of this part is hence crucial. Forthis purpose it is convenient to use the bit representationof an unsigned integer to code the basis states. Denotingthe integers with bit representations corresponding to thespin up and spin down occupations of a given basis statewith I ↑ and I ↓ , respectively, the integer representation ofthe basis state is I = I ↑ + 2 L I ↓ . With the binary rep-resentation of the basis states, the multiplication of theHamiltonian can be done efficiently using bitwise opera-tions.
2. Calculating the Green function
Having obtained the ground state, the Green functioncan now be calculated. From the Lehmann representa-tion it follows that it can be written as G rij ( ε ) = G eij ( ε ) + G hij ( ε ) (27)with the electron and hole Green functions defined by G eij ( ε ) = h Ψ N | c i ε − H + E N + iη c † j | Ψ N i (28)and G hij ( ε ) = h Ψ N | c † j ε + H − E N + iη c i | Ψ N i , (29)respectively. In the following we focus on the electronGreen function which is the matrix representation of theresolvent operator ( z − H ) − in the basis spanned bythe | i i = c † i | Ψ N i vectors. To obtain the i ’th diagonalelement, G eii ( ε ) = h i | ( z − H ) − | i i , (30)where z = ε + E N + iη , again the Lanczos algorithm isused to put H on a tridiagonal form, but this time theLanczos iterations are started from the normalized initialstate | φ i = | i i /b where b = h i | i i . Hence, in the gen-erated Krylov subspace the diagonal element in Eq. (30)corresponds to the matrix element b [( ε − H + E N + iη ) − ] of a tridiagonal matrix, which can be obtainedas the continued fraction G eii ( ε ) = b ε − a − b ε − a − b ε − a − · · · . (31)Again the Lanczos iterations are continued until the fre-quency dependent Green function element has converged. C. Von Neumann entropy
The following section demonstrates how a quantitativemeasure of the degree of correlations in a system can beobtained by considering the von Neumann entropy of thereduced single-particle density matrix ρ . The entropy isdefined by S [ ρ ] = − Tr[ ρ log ρ ] = − X n ρ n log ρ n , (32)where in the last equality ρ has been expressed in itsdiagonal representation, ρ = P n ρ n | n ih n | .In the basis of the atomic p z orbitals the matrix ele-ments of the reduced density matrix are given by (withthe spin index suppressed) ρ ij = h Ψ | c † j c i | Ψ i , (33)with the diagonal elements equal to the site occupations.In the diagonal representation ρ n thus represents the oc-cupation of the eigenstate | n i of the density matrix.We note that 0 ≤ S ≤ L log 2, where 2 L is the di-mension of the single-particle Hilbert space includingspin. The expression for S max follows because the num-ber of electrons equal L in all the systems, i.e. half filled“band”. When | Ψ i is a single Slater determinant (cor-responding to zero correlation) we have S = 0, and when | Ψ i has equal weight on a complete set of orthogonalSlater determinants (corresponding to maximal correla-tion) we have ρ n = 1 / n and thus S = L log 2.Thus the number 0 ≤ S/S max ≤ | Ψ i . IV. RESULTSA. Total energies
We first address the degree of correlation in the ex-act ground states by considering the von Neumann en-tropies of the corresponding density matrices. The cal-culated entropies are listed in Tab. I. Except for theHubbard description of benzene (see Sec. IV C) whichclearly presents strong correlations, the entropies of theground states are ∼
10% of their maximum value S max corresponding to weakly correlated systems. The finitevalues of the entropies reveal that none of the groundstates are single Slater determinants implying that theHartree-Fock ground state energies will be larger thanthe exact ones.We here follow the usual convention and define thecorrelation energy as the part of the total energy notincluded in Hartree-Fock, i.e. E corr = E exact − E HF . (34)Fig. 1 shows the exact correlation energies of the neu-tral molecules together with those obtained by evaluat-ing the total energy from Eqs. (20) and (21) with theself-consistently determined Green function and GW self-energy.For the series of molecules considered here the correla-tion energy constitute less than 0 .
5% of the total energies.Furthermore, as expected it decreases (in absolute size)with the number of atoms in the molecule. Clearly, theGW approximation performs reasonably well for all themolecules capturing on average 66% of the correlationenergy.
B. Spectral properties
For isolated systems such as molecules, true quasi-particles resembling single-particle excitations are char-
Formula
L S/S max E gap (eV)thiophene C H S 5 0.07 11.19pyridine C H N 6 0.11 10.61benzene C H H
12 0.10 9.24naphthalene C H
10 0.11 8.65anthracene C H
14 0.12 7.06OPV2 C H
14 0.10 8.30TABLE I: Chemical formula, number of p z orbitals ( L ) in-cluded in the PPP model and exact ground state entropies( S ) for the listed molecules. T h i o p h e n e P y r i d i n e B e n z e n e B i p h e n y l N a p h t h a l e n e A n t h r a c e n e O P V -2-1.8-1.6-1.4-1.2-1-0.8-0.6-0.4-0.20 E c o rr ( e V ) ExactGW
FIG. 1: (Color online) Exact and GW correlation energies ofthe neutral groundstate of the seven molecules. acterized by having a weight close to unity (for non-degenerate levels) in the spectral function, i.e. Z n = X i |h Ψ N +1 n | c † i | Ψ N i| ∼ . (35)This is equivalent to saying that there exists an orbital | ν i so that the excited state h Ψ N +1 n | can be written asthe single-particle excitation c † ν | Ψ N i . In Fig. 2 we showthe single-particle density of states (DOS), D ( ε ) = X i A i ( ε ) (36)for the OPV2 molecule on a logarithmic scale. The heightof the peaks reflects the value of Z n (modulo degenera-cies). The HF, and in particular, the GW approximationreproduce the lowest lying excitations quite well whilehigher excitations are poorly described. All the peaksin the HF spectrum have Z n = 1 while GW does shiftsome spectral weight from the main peaks to tiny satellitestructures (at higher energies than shown on the plot).However, the GW satellites do not correspond to featuresin the exact spectrum. This shows that excitations with Z n ≪
1, i.e. excitations which do not have single-particlecharacter, are not well described by GW whose main ef-fect is to improve the position of the HF single-particlepeaks.In the following we consider the lowest lying single-particle excitations of the molecules as obtained withHartree-Fock, G W and self-consistent GW. In theG W calculations the starting Green function G istaken to be the self-consistently determined Hartree-FockGreen function. Fig. 4 gives an overview of the calculatedexcitation energies relative to the exact ones. Energiescorresponding to electron removal and electron additionare located on the negative and positive half of the x -axis,respectively. From this plot clear trends in the calculatedexcitation energies emerge. ε (eV) Log D ( ε ) ( a r b . un i t s ) ExactHFGWOPV2
FIG. 2: (Color online) Single-particle DOS of the OPV2molecule. Note the logarithmic axis.
Within HF the occupied (unoccupied) levels are sys-tematically overestimated (underestimated), and the de-viation from the exact values worsens for the higher ly-ing excitations. A closer inspection of the figure revealsa few HF energies at ∼ ± ∼ ± not arise because HF gives a correct description of themany-body states and their energies. This was alreadyclear from the analysis above which showed that the exacteigenstates are not single Slater determinants and hencethe excitation energies in Eq. (2) have contributions fromboth ∆ relax and ∆ corr . The good agreement must there-fore be ascribed to cancellations between the relaxationand correlation contribution to the exact energies (this isdiscussed further in connection with Fig. 4).Both the G W and the GW give consistently betterenergies than HF – in particular for the higher lying ex-citations where the absolute errors are reduced to lessthan ∼ ∼ relax and ∆ corr to the excitation energies in Eq. (2),we plot in Fig. 4 the difference between the exact gapsand the gaps obtained from the (i) Hartree-Fock eigen-values, (ii) Hartree-Fock total energy differences with self-consistent relaxations in the N ± E gap = ε LUMO − ε HOMO can be expressed as E gap = ε HFLUMO − ε HFHOMO + ∆ gaprelax + ∆ gapcorr (37) -9 -8 -7 -6 -5 -4 -3 ε nexact (eV) -1-0.500.51 ∆ ε napp r o x ( e V ) HFG W [HF]GW Exact ref.
FIG. 3: (Color online) Energy of the 3 highest occupied and3 lowest unoccupied molecular orbitals relative to the exactvalues. While Hartree-Fock underestimates the occupied andoverestimates the unoccupied levels, self-consistent GW showsthe opposite trends but deviates on average less from the exactresult. where ∆ gaprelax and ∆ gapcorr are the gap equivalents of the cor-responding quantities in Eq. (2) and ε HFHOMO/LUMO arethe Hartree-Fock HOMO/LUMO eigenvalues. By defini-tion ∆ gaprelax is difference between the gaps obtained fromthe HF eigenvalues and relaxed HF total energy differ-ences. In Fig. 4 this is given by the vertical distancebetween the (blue) squares and circles. The correlationcontribution ∆ gapcorr can be read off as the difference be-tween the exact gap (dashed horizontal line) and the re-laxed HF total energy gap (blue squares). Inclusion ofrelaxation effects clearly reduces the HF gaps consider-ably implying that ∆ gaprelax <
0. This reduction is due tothe screening from the orbital relaxation which reducesthe Coulomb interaction with the added hole or electronand hence also the gap.In contrast to the HF (eigenvalue) gaps for which theagreement with the exact gap worsens as a function ofthe size of the molecules, the GW gaps follow more con-sistently the same trend and underestimates the exactgaps with 0 . − .
35 eV for all the molecules. The closeresemblance between GW and the relaxed HF result in-dicates that the effect of GW is mainly to account for thescreening effects included in HF via orbital relaxations,∆ relax . C. Long- versus short-range interactions
To demonstrate the shortcomings of the GW approx-imation for strongly correlated systems, we consider aHubbard model description of the benzene molecule. Itshould be noted that this Hubbard description of benzene T h i o p h e n e P y r i d i n e B e n z e n e B i p h e n y l N a p h t h a l e n e A n t h r a c e n e O P V -0.4-0.200.20.40.6 ∆ ε gap ( e V ) GWHF (eigenvalues)HF (with relaxations)
Exact ref.
FIG. 4: (Color online) The HOMO-LUMO gap relative to theexact values. In addition to the HF and GW single-particleenergies, the relaxed Hartree-Fock total energy differences, E HF0 ( N + 1) + E HF0 ( N − − E HF0 ( N ) are also shown. Theexcellent results of HF for the three smallest molecules is aresult of error cancellation between relaxation and correlationcontributions. is not intended as a realistic description of the benzenemolecule, rather it serves to illustrate the limitations ofthe GW approximation. The Hamiltonian is identical tothe PPP-Hamiltonian in Eq. (5), except that the longrange Coulomb interactions in the third term have beenomitted. The values for the hopping elements and theonsite Coulomb interaction are t = 2 .
539 and U = 10 . U/t -ratio of ∼ .
48 eV. This should be compared to 0 .
16 eV whichis the difference between the exact and the GW groundstate energy for the PPP description of benzene.The poor performance of both Hartree-Fock and GWfor the spectral properties of the Hubbard benzene is il-lustrated in Fig. 5 which shows the spectral function ascalculated with the two methods together with the exactone. Both Hartree-Fock and GW severely underestimatesthe position of the LUMO level and completely misses thedetails of the spectrum at higher energies.This clearly demonstrates that GW is of limited rele-vance when considering systems where correlation effects(∆ corr ) dominates over screening, or relaxation, effects(∆ relax ). ε (eV) Log D ( ε ) ( a r b . un i t s ) ExactHFGW Hubbard benzene
FIG. 5: (Color online) Single-particle DOS for the Hubbarddescription of the benzene molecule (only on-site interactionsfrom the PPP model are kept). Note the logarithmic axis.
D. Exact Kohn-Sham orbital energies
Within density functional theory (DFT) the eigenval-ues of the single-particle Kohn-Sham Hamiltonian, H KS ,are often interpreted as physical energies. In principlethe validity of this approximation depends on the sizeof the derivative discontinuity of the true and unknownexchange-correlation (xc-) functional . In practice theuse of semi-local xc-functionals represents an additionalapproximation. It is of general interest to investigateto what extent the disrepancy between KS energies andtrue QP energies result from the use of approximate func-tionals and to what extent this is a property of the ex-act functional. Below we compare the exact Kohn-Shamspectrum to the exact QP spectrum of the PPP benzenemolecule using the lattice version of DFT.The lattice version of DFT follows by extending thefundamental concepts of standard DFT, such as theHohenberg-Kohn theorem and the Kohn-Sham equa-tions, to model Hamiltonians as e.g. the PPP-Hamiltonian. In this reformulation of DFT the site oc-cupations n i replaces the continuous electron density n ( r )as the fundamental variable that determines the groundstate properties. The lattice version of the single-particleKohn-Sham Hamiltonian is given by the sum of the hop-ping terms (the kinetic energy) and a site dependentKohn-Sham potential V KS i which is constructed to yieldthe correct site occupations of the ground state, H KS = − X h ij i σ t ij c † iσ c jσ + X i V KS i ˆ n i . (38)For the present purpose the explicit form of the site po-tential V KS i is not important. The fact that the latticeversion of the Kohn-Sham potential is an onsite poten- tial, is equivalent to the restriction of the Kohn-Shampotential V xc ( r ) in the real-space formulation of DFT toa local potential.Due to the high symmetry of the benzene moleculeall sites in the PPP-Hamiltonian are equivalent implyingthat V KS i has the same value for all sites. Except for aconstant shift, the eigenvalues of the Kohn-Sham Hamil-tonian are therefore given by those of the hopping partof the Hamiltonian. The HOMO-LUMO gap calculatedfrom the Kohn-Sham eigenvalues is E KSgap = 5 .
08 eV whichis a severe underestimation of the true gap of 11 .
39 eV. Inline with previous studies we thus conclude that themain reason for the discrepancy between KS eigenvaluesobtained with approximate xc-functionals and the exactorbital energies is due to the derivative discontinuity ofthe exact xc-functional.
V. CONCLUSION
We have presented calculations for the total energy andcharged single-particle excitations in seven conjugatedmolecules described by the semi-empirical PPP modelwithin fully self-consistent GW and exact diagonaliza-tion. The results show that the GW approximation givesa consistently good description of both total energies andelectronic excitations with a slight tendency to overesti-mate (underestimate) the position of the latter for occu-pied (unoccupied) levels. We have found that the effectof the GW self-energy is similar to the inclusion of or-bital relaxations in the N ± VI. ACKNOWLEDGEMENTS
Financial support from the Danish Council for Pro-duction and Technology (FTP) under Grant 26-04-0181“Atomic scale modeling of emerging electronic devices”is acknowledged. KST acknowledges support from theDanish Center for Scientific Computing. The Center forAtomic-scale Materials Design (CAMD) is sponsored bythe Lundbeck Foundation. L. Hedin, Phys. Rev. , A796 (1965). M. S. Hybertsen and S. G. Louie, Phys. Rev. B , 5390(1986). W. G. Aulbur, L. Jonsson, and J. W. Wilkins, in
SolidState Physics , edited by H. Ehrenreich and F. Saepen (Aca-demic Press, New York, 2000), vol. 54, p. 1. F. Aryasetiawan and O. Gunnarsson, Rep. Prog. Phys. ,237 (1998). 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). 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). C. Freysoldt, P. Rinke, and M. Scheffler, Phys. Rev. Lett. , 056803 (2009). K. Kaasbjerg and K. Flensberg, Nano. Lett. , 3809(2008). K. S. Thygesen and A. Rubio, Phys. Rev. Lett. ,046802 (2009). C. D. Spataru, L. X. Benedict, and S. G. Louie, Phys. Rev.B , 205204 (2004). K. S. Thygesen and A. Rubio, J. Chem. Phys. , 091101(2007). P. Darancet, A. Ferretti, D. Mayou, and V. Olevano, Phys.Rev. B , 075102 (2007). K. S. Thygesen, Phys. Rev. Lett. , 166804 (2008). P. My"ohanen, A. Stan, G. Stefanucci, and R. vanLeeuwen, Europhys. Lett. , 67001 (2008). C. D. Spataru, M. S. Hybertsen, S. G. Louie, and A. J.Millis, Phys. Rev. B , 155110 (2009). M. P. von Friesen, C. Verdozzi, and C.-O. Almbladh,arXiv:0905.2061 (2009). N. E. Dahlen and R. van Leeuwen, Phys. Rev. Lett. ,153004 (2007). J. A. Pople, Trans. Faraday Soc. , 1375 (1953). R. Pariser and R. G. Parr, J. Chem. Phys. , 466 (1953). R. Pariser and R. G. Parr, J. Chem. Phys. , 767 (1953). A. Schindlmayr, T. J. Pollehn, and R. W. Godby, Phys.Rev. B , 12684 (1998). T. J. Pollehn, A. Schindlmayr, and R. W. Godby, J. Phys.:Condens. Matter , 1273 (1998). C. Verdozzi, R. W. Godby, and S. Holloway, Phys. Rev.Lett. , 2327 (1995). C. Rostgaard, K. W. Jacobsen, and K. S. Thygesen, Phys.Rev. B, submitted (2009). P. Fulde,
Electron Correlations in Molecules and Solids ,no. 100 in Springer Series in Solid-State Sciences (Springer,Berlin, 1995), 3rd ed. K. Ohno, Theor. Chim. Acc. , 219 (1964). I. R. Ducasse, T. E. Miller, and Z. G. Soos, J. Chem. Phys. , 4094 (1982). W. Barford and R. J. Bursill, Chem. Phys. Lett. , 535(1997). R. J. Bursill, C. Castleton, and W. Barford, Chem. Phys.Lett. , 305 (1998). W. Barford, R. J. Bursill, and M. Y. Lavrentiev, J. Phys.:Condens. Matter , 6429 (1998). M. Y. Lavrentiev, W. Barford, S. J. Martin, and H. Daly,Phys. Rev. B , 9987 (1999). K. S. Thygesen and A. Rubio, Phys. Rev. B , 115333(2008). P. Pulay, Chem. Phys. Lett. , 393 (1980). A. L. Fetter and J. D. Walecka,
Quantum theory of many-particle systems (McGraw-Hill, New York, 1971). E. Dagotto, Rev. Mod. Phys. , 763 (1994). J. P. Perdew and M. Levy, Phys. Rev. Lett. , 1884(1983). R. W. Godby, M. Schluter, and L. J. Sham, Phys. Rev.Lett. , 2415 (1986). K. Schonhammer, O. Gunnarsson, and R. M. Noack, Phys.Rev. B , 2504 (1995). P. Dayal, M. Troyer, and R. Villiger,
The It-erative Eigensolver Template Library (IETL) , . In the remaining of the paper HOMO and LUMO abbrevi-ations will be used to refer to the ionization potential andelectron affinity also outside Hartree-Fock theory. The present work has used the implementation from theIETL project44