Koopmans' condition in self-interaction corrected density functional theory
Peter Klüpfel, Mai Phuong Dinh, Paul-Gerhard Reinhard, Eric Suraud
KKoopmans’ condition in self-interaction corrected density functional theory
P. Kl¨upfel,
1, 2, 3
P. M. Dinh,
2, 4
P.-G. Reinhard, and E. Suraud
2, 4 CNRS, Laboratoire Collisions Agr´egats R´eactivit´e (IRSAMC), F-31062 Toulouse, France CNRS, Laboratoire de Physique Th´eorique (IRSAMC), F-31062 Toulouse, France Icelandic Center of Computational Science, University of Iceland, 107 Reykjav´ık, Iceland Universit´e de Toulouse, UPS, Laboratoire de Physique Th´eorique (IRSAMC), F-31062 Toulouse, France Institut f¨ur Theoretische Physik, Universit¨at Erlangen, D-91058 Erlangen, Germany (Dated: August 10, 2018)We investigate from a practitioner’s point of view the computation of the ionization potential(IP) within density functional theory (DFT). DFT with (semi-)local energy-density functionals isplagued by a self-interaction error which hampers the computation of IP from the single-particleenergy of the highest occupied molecular orbital (HOMO). The problem may be cured by a selfinteraction correction (SIC) for which there exist various approximate treatments. We compare theperformance of the SIC proposed by Perdew and Zunger with the very simple average-density SIC(ADSIC) for a large variety of atoms and molecules up to larger systems as carbon rings and chains.Both approaches to SIC provide a large improvement to the quality of the IP if calculated from theHOMO level. The surprising result is that the simple ADSIC performs even better than the originalPerdew-Zunger SIC (PZSIC) in the majority of the studied cases.
I. INTRODUCTION
Density Functional Theory (DFT) [1–3] has becomea standard theoretical tool for the investigation of elec-tronic properties in many physical and chemical systems.It provides fairly reliable results with moderate computa-tional effort. Practical implementations of DFT employsimple and robust approximations for the exchange andcorrelation functional. The simplest one is the Local Den-sity Approximation (LDA), which has been proven veryuseful in calculations of electronic structure and dynam-ics. Typical applications are reaching from first principlecalculations on the electronic ground state and molecu-lar geometries [4], over dynamic studies of near equilib-rium situations (e.g. optical response, direct one-photonprocesses) to highly non-linear dynamical scenarios [5–7].However, LDA is plagued by a self-interaction error [8]:The Kohn-Sham (KS) mean field is computed from thetotal density which includes all occupied single-particlestates, including the state on which the LDA field ac-tually acts. Locality of the energy functional leads toa wrong asymptotic KS field. This is still a great hin-drance in many applications, for instance a possibly largeunderestimation of IP and the absence of Rydberg or ex-citonic series in the static KS spectrum [9, 10], the polar-izability in chain molecules [11, 12] or the spectral andfundamental gap in solids [13, 14]. Another challengingapplication is the description of molecules or clusters de-posited on surfaces [15, 16]. In dynamic scenarios, theself-interaction error also dramatically affects ionizationdynamics, especially close to thresholds, e.g., in a time-dependent DFT model of electron emission [7, 17–19].In practice, the wrong asymptotics of the KS fieldstems from an incomplete cancellation of the self-interaction error between the Hartree potential and theapproximate exchange and correlation field. Such a spu-rious self-interaction is avoided by a complicated non-locality in exact KS-DFT [20–22]. For the total energy, the requirement for non-locality can be incorporated intogradients of density leading to the generalized gradientapproximation (GGA) [23–25]. This approximation in-deed served to lift DFT to a quantitative level in molec-ular physics and chemistry problems, but is insufficientto restore proper asymptotics of the mean field.Although there are approaches to improve the asymp-totic KS potential [26] those are often too demanding forpractical calculations, in particular in the time-domain.The aforementioned examples show that there is still aneed for robust and practical approaches to improve onthe asymptotic KS potential, such as self-interaction cor-rection (SIC).The original proposal for a SIC [8, 27] by Perdew andZunger (PZ) has been developed at various levels of re-finements and proved to be useful over the years, in par-ticular for structure calculations in atomic, molecular,cluster and solid state physics, see e.g. [28–35]. Thisoriginal PZSIC scheme, however, leads to an orbital-dependent mean field which causes several formal andtechnical difficulties [29, 32, 33, 35, 36]. There are at-tempts to circumvent the orbital dependence by treatingSIC with optimized effective potentials (OEP)[37], fora review see [38]. However the resulting formalism is,again, quite involved and usually treated approximately[39]. A very robust and simple SIC is the average densitySIC (ADSIC) which was proposed already very early [40],taken up in [41], and has been used since in many appli-cations to cluster structure and dynamics in all regimes.ADSIC takes the total density divided by the electronnumber as a reference for the single electron density ineach state. Non-locality is incorporated in the schemeby the global density integral providing the total parti-cle number. However the ADSIC functional, unlike thePZSIC one, is a functional of the total density and thus,the subsequent mean field is not orbital dependent any-more. Having the correct total charge, ADSIC providesthe proper asymptotics for the mean field. It is argued a r X i v : . [ phy s i c s . a t m - c l u s ] A ug that the approximation by one and the same (average)single-particle density for each state is only applicablein simple metals where all electronic states cover thesame region of space, that is, in the case of a cluster,the whole cluster itself [41, 42]. Later studies revealedthat ADSIC is also an efficient correction scheme for non-metallic systems with delocalized electrons, such as or-ganic molecules [43, 44].The aim of this paper is to investigate the performanceof ADSIC in direct comparison to PZSIC for a large va-riety of atoms and molecules in their ground state. Thesample covers systems of different binding types, and notonly metallic ones. We will compare ADSIC with a mereDFT treatment using (semi-)local functionals and withPZSIC, also occasionally with Hartree-Fock. The com-parison focuses on the proper description of the IP. Westart from atoms as elementary building blocks of anymolecule, proceed to a large variety of molecules fromsimple dimers to more complex organic structures, andfinally discuss carbon rings and chains with a systematicvariation of sizes. Such a strategy allows us to cover var-ious binding types but also various geometries and evendimensionality. II. IONIZATION POTENTIALA. Definitions
The key quantity of this survey is the ionization energy I , commonly called ionization potential (IP). The IP ofa N electron system is given by the energy difference I ≡ I ∆ = E ( N − − E ( N ) . (1)The energies E ( p ) correspond to ground-state configura-tions of a p particle system in a given external potential,typically the Coulomb potential created by the nuclearcharges. Both energies, E ( N ) as well as E ( N − n is related to the IP by n ( r ) | r |→∞ ∼ exp (cid:104) − √ I | r | (cid:105) . (2)In an exact KS-DFT, the asymptotic decay of the totaldensity is defined merely by the highest occupied KS or-bital, i.e., the HOMO [4]. In combination with the properasymptotic behavior of the KS potential ( v s ( r ) → I ≡ I ε = − ε HOMO . (3)For an exact exchange-correlation functional, both def-initions of the ionization energy, i.e., Eqs. (1) and (3),coincide, i.e. they obey I ∆ = − ε HOMO . The identifi-cation of the negative HOMO energy with the IP wasreferred to as Koopmans’ theorem [45] long before thefundamental concepts where extended rigorously to theframework of DFT [46, 47]. For approximate energy func-tionals Koopmans’ condition does not necessarily hold [2]and the deviation from the ideal behavior can be used todefine the Non-Koopmans (NK) energy∆ E NK = I ε − I ∆ . (4)A value ∆ E NK = 0 signals that Koopmans’ theorem isfulfilled. In such a situation, the properties of the HOMOlevel are closely related to ionization and electron emis-sion. We know that LDA produces rather large violationsof Koopmans’ theorem and thus exhibits sizable ∆ E NK .SIC should reduce that, and the amount of reduction isone measure of the efficiency of the actual SIC scheme.It is also interesting to compare the performance ofa calculation with respect to data. Thus we considerin addition to the NK energy the bare error in the IPrelative to experiments or other reference data. As wehave two definitions of the IP, there are two bare errorsin an approximate theory :∆ I ε = I ε − I ref , ∆ I ∆ = I ∆ − I ref . (5)An experimental reference energy may be hampered byuncertainties, as ionic relaxation throughout the ioniza-tion process can lead to situations which lie between theidealized vertical IP (for very fast ionization) and hori-zontal one (for very slow ionization). Here only verticalIP in the ground-state geometry is considered. Reliableatomic coordinates for small molecules are given by theMP2 optimized structures as provided in the G2 dataset[48]. G2 theoretical and experimental energies may differby several tens to hundrets of meV [49]. This error canbe considered negligible on the scale of the expected er-rors, stemming from the approximate nature of the usedexchange-correlation functionals and the use of pseudo-potentials [50]. Experimental data for vertical IP there-fore appear as a safe choice of reference [51].In contrast to the errors in the IP (5), the NK energydoes not require any reference data that may be ham-pered by experimental uncertainties. It therefore pro-vides a rather rigorous criterion for the quality of energyfunctional approximations. B. Impact of a proper description of IP
The two definitions (1) and (3) for the IP are equallyjustified in an exact calculation. However, for (semi-)local functionals (as in LDA and GGA), it is usuallyfound that only the energy difference (1) can be usedto extract a good estimate for the IP. The estimate (3)from the single-particle spectrum requires the proper 1 /r -asymptotics of the exchange-correlation potential (forneutral systems). This is not provided in calculationsbased on typical semi-local functionals.While energy differences often allow reliable estimatesalready with semi-local functionals, they require two cal-culations, which is more involved than a straightforwardextraction from the HOMO level. This alone would notbe a priori a major hindrance. But there are manysituations where the extraction via energy difference isnot an option : In periodic calculations (e.g., on sur-faces), a rigorous calculation of the IP (called work func-tion in this case) from an energy difference is hard toachieve because one cannot easily model a single excesscharge in a periodic setup. The same situation appliesfor calculation of band-gaps in solids. In dynamical situ-ations, as described by time-dependent DFT, an accuratemodeling of the ionization process requires an accuratestatic single-particle spectrum. As the propagation of thesingle-particle states is driven by the time-dependent KSHamiltonian, the energy differences in its spectrum anda proper position of the IP is more essential than thetotal energy. Ionization properties are also mostly de-fined by the HOMO level which thus has to be correctlydescribed.A way to illustrate the self-interaction error is toconsider the energy E ( ν ) as a function of a fractionalparticle number ν . The ionization process goes along ν = N −→ N −
1. An exact functional produces a linearbehavior [20], as : E ( ν ) = (1 − x ) E ( N )+ xE ( N − , x = N − ν ∈ [0 , . (6)A similar linear behavior is also observed for ionizationfrom an anion to the neutral system ( N + 1 −→ N ).This exact E ( ν ) is shown as (red) solid line in Figure1. The remarkable feature is a discontinuous derivative,that is a kink, at ν = N . Semi-local functionals dealwith smooth functions (no kinks, no discontinuities) andproduce smooth trends as shown in the (green) dottedline. The definition (3) of the IP through the HOMOenergy coincides with the slope of the total energy forfractional particle numbers at ν (cid:37) N : I ε = − dEdν (cid:12)(cid:12)(cid:12)(cid:12) ν (cid:37) N (7)as indicated in the figure. The linear trend of the exactenergy naturally guarantees I ε = I ∆ , while the convexcurve from LDA necessarily implies I ε < I ∆ .The (blue) dashed line in Figure 1 finally shows re-sults from exact exchange in Hartree-Fock (HF). Full HF FIG. 1. (Color online) Illustration of the ground state en-ergy E as a function of fractional occupation number ν . TheIP from the HOMO level, I ε , corresponds to the left-handedderivative (slope = black line) of the energy at ν = N . TheIP from energy differences I ∆ , is associated with the exactlylinear behavior (red/solid line). The LDA result provides asmooth curve (green/dotted line). Hartree-Fock (blue/dashedline) also has a discontinuous derivative at ν = N as the exacttrend, but tends to overestimate the kink. is free from self-interaction error. Thus it qualitativelyyields the correct result, namely the kink at ν = N . Ithowever differs from the linear trend in between the in-teger particle numbers. This leads to an overestimationof the IP from the HOMO, I ε > I ∆ . We will address thisquestion in the last example of carbon chains, see Figure8. III. VARIOUS SCHEMES FOR A SIC
As SIC is rather an ad-hoc measure to cure the self-interaction problem, various recipes and approximationsare used, depending on the field of application. In thissection, we briefly summarize the PZSIC and ADSICwhich we will use later on in the extensive comparisonof results.
A. Perdew-Zunger SIC
As already mentioned in the introduction, a very pop-ular approach to the definition of a one-particle self-interaction error and a corresponding correction was pre-sented by Perdew and Zunger [8]. The self-interaction er-ror is given by accumulating the contributions from theindividual orbital densities n i ( r ) = | ϕ i ( r ) | for a set ofsingle-particle states ϕ N = ( ϕ , . . . , ϕ N ). It reads E SI [ ϕ N ] = N (cid:88) i =1 ( E H [ n i ] + E xc [ n i ]) , (8a)where E H is the Coulomb Hartree energy and E xc thedensity functional for exchange and correlations. Notethat this is not a functional of density alone. In fact, E SI [ ϕ N ] depends on the detailed orbitals. The PZSIC isdefined by subtracting the self-interaction error from theoriginal functional, i.e., E PZSIC [ ϕ N ] = E H [ n ] + E xc [ n ] − E SI [ ϕ N ] , (8b)where n = (cid:80) Ni =1 n i is the total electronic density.The mean-field equations are derived in straightfor-ward manner by variation of the SIC energy E PZSIC withrespect to the occupied single-particle orbitals ϕ i . Itturns out that mean-field Hamiltonian depends explicitlyon the particular single-particle state on which it acts.This emerges because the PZSIC energy functional is notinvariant under unitary transformations amongst the oc-cupied states. There are several ways to deal with sucha state-dependent Hamiltonian [28–31]. A particularlyefficient way is to use two different sets of single-particlestates which are connected by a unitary transformationamongst occupied states. That is actually the solutionscheme which we are using, for details see [52]. In thisapproach, the HOMO level is defined as usual, in thebasis-set which diagonalizes the Hamiltonian matrix. B. Average Density SIC
The average density SIC (ADSIC) starts from the SICenergy (8b) and simplifies it by assuming that indis-tinguishable electrons are represented by equal single-particle densities. In such an extreme simplification, oneexpresses them as the one-particle fraction of the totalspin-density n i ( r ) = n σ i ( r ) /N σ i where σ i is the spin ofstate i and N σ i the number of particles with spin σ i . Insuch a scheme, the standard PZSIC functional is repre-sented by the ADSIC functional : E ADSIC [ n ↑ , n ↓ ] = E H [ n ] + E xc [ n ] − (cid:88) σ ∈{↑ , ↓} N σ ( E H [ n σ ] + E xc [ n σ ]) (9)where n = n ↑ + n ↓ . This is a spin-density functionaland can be treated in the same manner as any LDA orGGA scheme. This makes it extremely simple and ef-ficient to use in atomic and molecular systems. How-ever, the ADSIC functional contains a cumbersome non-locality as it explicitly depends on the particle number N σ = (cid:82) d r n σ ( r ). This inhibits an application in peri-odic systems, where N σ is infinite. IV. RESULTSA. Numerical scheme and pseudo-potentials
The calculations use a representation of the single-particle wave functions on a coordinate-space grid with a spacing of 0.2 ˚A. Densities and fields where representedon a refined grid of 0.1 ˚A to account for the higher Fouriercomponents in products of single-particle states. Thecore electrons are handled within the frozen-core approx-imation by a real-space implementation of the projectoraugmented wave (PAW) method [53] using a develop-ment version of GPAW [54]. The projectors and partialwaves of the PAW method are taken as provided withinthe GPAW repositories for bare LDA exchange and cor-relation, i.e. without accounting for a SIC.This corresponds to use pseudo-potentials developedfor LDA applications in the context of PZSIC or ADSICwithout readjustment of the pseudo-potential parame-ters. This minor inconsistency is acceptable in variousapplications of SIC [32, 55, 56]. Here such an improve-ment is avoided in favor of using a unique set of pseudo-potentials for all energy functionals.For the following survey, we show results from LDAusing the PW92 parameterization [57]. For most of theexamples below, we have also performed GGA calcula-tions with the PW91 functional [24]. Even if the GGAslightly improves the overall quality of the IP, in partic-ular if calculated from energy differences, with very fewexceptions, the effect of the gradient dependence is lessthan 0.5 eV. Thus it neither affects the overall magnitudeof errors or change the general trends that are discussedin the following sections. We therefore focus on the LDApart in this survey.
B. Atoms
The first step is to investigate the performance of bothSIC approaches for atoms. The latter ones are the basicbuilding blocks of molecules and solids. Thus they mustbe correctly described before we can proceed to morecomplex scenarios. The electronic structure of atoms in-corporates single-electron states with similar shape butdifferent spatial extensions. Thus atoms are a criticaltest case for SIC which is known to strongly depend onthe level of localization.Figure 2 shows the IP as such for neutral atoms fromhydrogen ( Z = 1) to argon ( Z = 18). All methods yieldvery similar IP if it is evaluated as I ∆ , i.e. as the en-ergy difference (1). Results differ more for I ε computedfrom the HOMO according to Eq. (3). Here, the bare(semi-)local energy functionals underestimate the ioniza-tion energy of by 30-40%. The defect is well known andcan be traced back to the wrong asymptotic behavior ofthe exchange-correlation potential for | r | → ∞ [8]. Ob-viously, both SIC approaches cure the problem and yieldexcellent agreement with experimental data. ADSIC isslightly superior in case of open shell atoms, while PZSICslightly overestimates the IP. Accounting for GGA (notshown here) has an insignificant effect for both I ∆ as wellas I ε . The results do not sufficiently differ to justify aseparate plot.Figure 3 shows the same data of figure 2 but in terms HH eL i B e B CN O F N e N a M g A l S i PS C l A r I ε ( e V ) exp.ADSICPZSICLDA 0 5 10 15 20 25 I ∆ ( e V ) FIG. 2. (Color online) Ionization potentials, I ∆ from Eq. (1),and I ε from Eq. (3), for neutral atoms from hydrogen to argonand different approaches to self-interaction correction : Av-erage density SIC (squares), Perdew-Zunger SIC (diamonds),and the uncorrected local-density approximation (open cir-cles). Experimental data are displayed as closed circles [51]. of errors with respect to experimental data and of theNK energy. This reveals some differences between PZSICand ADSIC where, somewhat surprisingly, the techni-cally much simpler ADSIC visibly yields smaller NK en-ergies and errors ∆ I ε .At this point, it is worth recalling that the ADSICscheme can be derived as an approximation to PZSIC,assuming a most delocalized representation of the single-particle densities. In ADSIC, orbitals are extending overthe whole atom in stark contrast to the localized orbitalsthat are commonly found in PZSIC calculations [33]. Itthus appears that a significant higher level of delocaliza-tion is actually desirable, which confirms previous con-cerns that PZSIC orbitals are too localized. C. Simple Molecules
As a next step, we consider simple molecules, as manydimers, and a few more complex ones. The selection hasbeen adapted from [25]. It covers systems which do nothave the problem of spatial symmetry breaking by anunrestricted mean-field calculation. Reference data wastaken from [51].Figure 4 shows the IP for a chosen set of molecules. Atfirst glance, the results resemble those for atoms in Fig-ure 2. The I ∆ provides reasonable results for all methodswhile I ε shows dramatic differences between LDA and theSIC models. However, taking a closer look, we also seethat results for ADSIC and PZSIC show larger differ- -8-6-4-2 0 2 HH eL i B e B CN O F N e N a M g A l S i PS C l A r ∆ E N K ( e V ) -6-4-2 0 2 4 6 ∆ I ∆ ( e V ) ADSICPZSICLDA-8-6-4-2 0 2 ∆ I ε ( e V ) FIG. 3. (Color online) Errors in calculated IP compared withexperimental values for the series of atoms depicted in Figure2. Upper and middle panels : error from I ε and I ∆ respec-tively, according to Eq. (5). Lower panel : non-Koopman’senergy defined in Eq. (4). H L i H CH NH O H H O H F L i L i F B e C H C H HCNC O N N O O F P C l I ε ( e V ) I ∆ ( e V ) exp.ADSICPZSICLDA FIG. 4. (Color online) Same as in Figure 2, but for a set ofsimple molecular systems. ences than in the case of atoms. Somewhat surprisingly,ADSIC comes again much closer to experimental datathan PZSIC.Figure 5 shows the data from the previous figure interms of energy differences, the errors ∆ I as comparedto reference data and the NK energy. The results con- -8-6-4-2 0 2 4 H L i H CH NH O H H O H F L i L i F B e C H C H HCNC O N N O O F P C l ∆ E N K ( e V ) -6-4-2 0 2 4 6 ∆ I ∆ ( e V ) ADSICPZSICLDA-8-6-4-2 0 2 4 ∆ I ε ( e V ) FIG. 5. (Color online) Same as in Figure 3, but for a set ofsimple molecular systems. firm the impressions indicated in the comparison of theIP as such : The ∆ I ε shows significant differences be-tween PZSIC and ADSIC since the latter one generallyperforms better.One may argue that comparison with reference datais also influenced by other details of the calculations orthe choice of the reference data. The NK energy (lowestpanel) is free of these uncertainties. ADSIC clearly de-livers the smallest NK energies. This was seen alreadyfor atoms. But here in the case of molecules the effect iseven more pronounced as PZSIC shows larger deviations. D. Systematic sets of molecules
In this section, we look at a systematic variation ofmolecules around basic carbohydrates. The first family(CH x ) represents a variation of the number of C-H bonds. The second family changes the character (single, double,triple bonds) in C H n . The third family is similar tothe first one but replacing carbon by a heavier element(silicon) with the same number of valence electrons, whilethe fourth series replaces the carbon atom by nitrogen(which has a different valence). In the final series one ofthe single bonded hydrogen atoms in CH is substitutedby a different group. The last element of the third andfourth series are not strictly within the systematics.We have seen in the previous systems that energy dif-ferences are showing more details than the energies assuch. We thus proceed here immediately to energy dif-ferences which are compiled in Figure 6. The results are -8-6-4-2 0 2 4 CHCH CH CH C HC H C H C H C H C H S i H S i H S i H S i H NHNH NH N H CH - S HCH - O HCH - C l CH - CH CH - H ∆ E N K ( e V ) -6-4-2 0 2 4 6 ∆ I ∆ ( e V ) ADSICPZSICLDA-8-6-4-2 0 2 4 ∆ I ε ( e V ) FIG. 6. (Color online) Same as in Figure 3, but for the familiesof molecules with systematically varied properties. very similar to the previous case of simple molecules, seeFigure 5. Some of the deviations are however larger thanin the previous case. This indicates that these complexmolecules are more critical test cases. Even in this moredemanding scenario, we find again that ADSIC performssuperior with respect to the deviation from reference dataand even more so for the NK energies.Thus we find that the ADSIC which assumes orbitaldensities that are delocalized over the whole moleculeyields a systematic improvement over PZSIC. This issomehow surprising in view of the deficits of ADSIC, inparticular as its inability to describe dissociation and thelack of size consistence are directly attributed to a toohigh level of delocalization.For infinite matter, ADSIC is not applicable due to theexplicit dependence on the particle number. Already forlarger systems, the explicit dependence on the total par-ticle number quickly renders the SIC contribution to theenergy functional an inefficient approach to cure prob-lems of the LDA. The observation that delocalization onthe length-scale of small molecules is in fact favorablefor the quality of the NK energy and IP calls for moresystematic investigations.
E. Carbon rings and chains
The self-interaction error on the IP for the CoulombHartree term is typically of order of e /R where R isthe radius of the system. The error for the exchange-correlation potential can be estimated within ADSIC as v xc [ n/N ]. Both shrink with increasing system size. Inorder to explore the evolution of the self-interaction er-rors with increasing size, we consider carbon rings andchains. For the latter ones, we only consider odd num-bers of atoms because only these have stable electronicconfigurations for spin saturated ground states. The car-bon atoms have more or less constant bond length. Thismeans that increasing the number of carbon atoms in-duces a (linear) growth of the geometrical extension, ei-ther of the chain or the ring.The upper two panels of Figure 7 show the IP forcarbon rings as a function of the number of atoms.Comparing ADSIC and LDA, we see again the equallygood performance for I ∆ , and the large self-interactionerror in I ε for LDA while ADSIC behaves very well. Thereference data, here calculated LDA values I ∆ , show apronounced step structure due to the successive filling ofthe electronic shells. Large I indicate particularly sta-ble electronic structures, i.e., shell closures. The suddenreductions shows that a new, and less bound, electronicshell has to be opened to place the given number of elec-trons. LDA and ADSIC reproduce the shell effect. Onfirst glance, the PZSIC results are quite surprising sincethey deviate even qualitatively from the other results,as they show less pronounced shell effects, at least withincreasing chain length.It shall be noted, that missing points also indicate thatreliable minimization of the PZSIC energy becomes chal-lenging for anionic configurations, where various localminima exist. The local minima correspond to differ-ent, almost energetically equivalent, configurations withdifferent levels of delocalization of the excess electron inthe spin-majority channel. The effect is worse for mid-shell systems but less problematic for closed shell config-urations. No such problem exists for ADSIC due to theabsence of the orbital-dependence.The lowest panel of Figure 7 shows the NK energies. -4-3-2-1 0 1 2 3 4 5 10 15 20 25 30 ∆ E N K ( e V ) number of atoms 0 2 4 6 8 10 I ∆ ( e V ) ADSICPZSICLDA 0 2 4 6 8 10 I ε ( e V ) I ∆ (LDA) FIG. 7. (Color online) Non-Koopmans energies (bottom), andionization potentials I ∆ (from energy differences, middle) and I ε (from the HOMO, top), computed in various schemes forcarbon rings of various size (5 ≤ N atoms ≤ I ∆ from LDA (see middle panel) is superimposed tothe I ε calculated in LDA, ADSIC and PZSIC. The ∆ E NK from LDA starts large but shrinks with in-creasing size as one could have expected. The ADSICresult is small throughout, but has a slight tendency toincrease with size, and of course, never becoming largerthan the error from LDA. However, the ∆ E NK fromPZSIC is generally large and even grows with system size.This finding is rather cumbersome, as it confirms that thedifference between the behavior of LDA and ADSIC onthe one hand, and PZSIC on the other hand, actuallystems from misconceptions in the PZSIC partially com-pensated in the approximate ADSIC.The significant and positive NK energy indicates thatstrong correlation effects, which are underestimated insemi-local exchange and correlation, are overestimatedby the PZSIC. The screening of such strong correlationeffects has to be reintroduced in the self-interaction cor-rected approach, e.g., by the assumption of more delo-calized states, as in case of the ADSIC.The convergence of ADSIC and LDA yet illustrates thecollapse of ADSIC as a working SIC scheme for extendedsystems, where N (cid:29)
1, contrary to the case of small N where the NK energy is still improved significantly. Thealmost constant but finite NK energy indicates that, al-though ADSIC is not capable of a complete curing of thenon-linear dependence of the LDA energy functional forfractional occupation, it at least provides a scheme thatyields similar magnitudes of errors for compact systemsand extended ones, whenever LDA is by itself considereda reasonable approximation there.Figure 8 shows IP and ∆ E NK for carbon chains. In -4-2 0 2 4 2 4 6 8 10 12 ∆ E N K ( e V ) number of atoms 0 2 4 6 8 10 12 14 I ∆ ( e V ) ADSICPZSICHFLDA 0 2 4 6 8 10 12 14 I ε ( e V ) I ∆ (LDA) FIG. 8. (Color online) Same as in Figure 7 but for linear car-bon chains (3 ≤ N atoms ≤ this case, we added also results from a pure Hartree-Fock(HF) calculation. PZSIC looks more agreeable here thanfor rings. But note that we have considered rather shortchains. There are again large differences between PZSICand ADSIC. This time, however, they are distributed al-most symmetrically around zero error (see lowest panel).No clear preference can be deduced in this example.The largest errors appear here for exact exchange inHF. Starting out perfect for the smallest chain C , TheNK energy jumps already for C and continues to grow further. This sounds, at first glance, very surprising asHF is free of any self-interaction error. However, removalof one electron causes polarization effects on the meanfield of the remaining electrons. These may be small incompact molecules. But polarizability grows huge par-ticularly in chains. The missing correlations from polar-ization effects are the source of the increasing NK errorwith HF. This can also be seen from the IP as such. TheHF result deviates much from the others for I ∆ (middlepanel), while it nicely stays in between for I ε . The miss-ing polarization effects explain the mismatch for HF. Thefact that the DFT based methods perform better indi-cate that some polarization effects are properly modeledin DFT, although it is also known that DFT underesti-mated the polarizability in some chain molecules [58]. F. Discussion
To summarize the results presented in the above fig-ures, we have computed average errors for each group ofsystem considered, atoms, simple molecules, and familiesof systematically varied molecules. Thereby, we distin-guish between mean error, mean absolute error and theerror fluctuations defined asME(∆ O ) = 1 N samp (cid:88) i ∆ O i , (10a)MAE(∆ O ) = 1 N samp (cid:88) i | ∆ O i | , (10b) σ (∆ O ) = 1 N samp (cid:88) i | ∆ O i − ME(∆ O ) | , (10c)where O is one of the considered observables, that is I ∆ , I ε , or E NK . The index i runs over the N samp samples ina given group, and ∆ O i = O i − O (ref) i stands for the ob-servables deviation from the reference data O (ref) i . Theresulting averages for each group are listed in Table I.Computation of IP as I ∆ , i.e. from energy differences, isalways a safe procedure yielding reliable results alreadywith LDA. Computation as I ε via the HOMO is possi-ble with good accuracy in both SIC models. The greatsurprise is that the very simplistic ADSIC approach per-forms very well for the I ε , typically even better thanPZSIC. The same conclusion is deduced from the non-Koopmans energy ∆ E NK . This was already seen fromthe above figures and is corroborated in Table I on aquantitative level.The excellent performance of ADSIC from compactsystems both in terms of accuracy and the small vio-lation of Koopmans’ condition is remarkable. Still, oneshould keep in mind the known deficiencies of the ap-proach. Most notably is the violation of size consistencywhich becomes apparent in the dissociation of a molecule.Consider a dimer with total electron number N whichdissociates into one part containing N electrons and an-other one with N electrons. The ADSIC for the com-pound involves, of course, the total electron number N . I ∆ I (cid:15) E NK ME MAE σ ME MAE σ ME MAE σ atomsLDA 0.2 0.3 0.3 -5.0 5.0 1.5 -5.1 5.1 1.6PZSIC 0.3 0.3 0.3 0.7 0.8 0.7 0.4 0.5 0.5ADSIC 0.4 0.5 0.4 0.2 0.4 0.4 -0.3 0.4 0.3small moleculesLDA 0.3 0.4 0.4 -4.6 4.6 0.7 -4.9 4.9 1.0PZSIC 0.5 0.6 0.5 1.4 1.5 0.9 0.9 1.0 0.7ADSIC 0.6 0.7 0.5 0.4 0.6 0.6 -0.2 0.4 0.3systematic mol.LDA 0.3 0.4 0.3 -4.1 4.1 0.5 -4.3 4.3 0.5PZSIC 0.5 0.6 0.4 1.4 1.4 0.7 0.9 0.9 0.4ADSIC 0.5 0.6 0.4 0.1 0.5 0.5 -0.5 0.5 0.1TABLE I. Mean error (ME), mean absolute error (MAE) anderror fluctuations σ as defined in eqs. (10) for IP as wellas NK energy for the data sets shown in figures 3, 5 and 6.Redundant data in 6 is only considered once in the averages.. Since we follow the dissociation path continuously, wenecessarily have to keep using N in the correction. Af-ter all, we end up with two isolated atoms which wouldbe treated by one common correction still regulated bythe total N . This is, of course, wrong as we know thateach single atom has to be separately corrected with itsown N i . The case is even worse in violent dynamics lead-ing to multi-fragmentation. The problem could alreadyhave been spotted from the fact that the dependence on N = (cid:82) d r n ( r ) implies a non-locality which becomes in-creasingly itching if n ( r ) ceases to be compact, but israther distributed over several regions of space.Fully accomplished dissociation and multi-fragmentation are, of course, extreme limits. Thedefects of ADSIC in this respect tend to show up earlier,for example, in the Born-Oppenheimer energies alongthe dissociation path. Thus one should not use ADSICfor computing large-amplitude molecular vibrationswithout careful checking its range of validity for thegiven application. Problems may also show up inmolecules which combine very different length scales as,e.g., in NaH O where the Na atom adds a rather diluteelectron distribution to the otherwise compact H O.In spite of the encouraging results presented above,one should check the NK energy ∆ E NK for each newapplication again.These known shortcomings should not hinder us to ap-preciate the good performance attained by ADSIC instructural and low energy dynamical situations. As il-lustrated all along the present work, ADSIC provides aremarkably robustness in terms of Koopmans’ violation.This implies, in particular, that it can be safely usedin the perturbative dynamical regime where only a tinyfraction of an electron is emitted. V. SUMMARY AND CONCLUSIONS
We have compared the performance of two differ-ent approaches to self-interaction correction regardingcalculated ionization potentials and violation of Koop-mans’ theorem. We have focused the discussions on twoSIC procedures: the original Perdew-Zunger approach(PZSIC) and the average density version thereof (AD-SIC). A wide range of electronic systems has been con-sidered ranging from atoms and simple molecules up tosystematics of moderate size molecules, in particular car-bon systems. The overall survey is thus quite general, sothat the conclusions attained have a safe ground, beyondany specific effect.We find in all examples considered here that ADSICprovides more reliable estimates of IP and a smaller vio-lation of Koopmans’ theorem. This is a welcome result inview of the remarkable simplicity (and correlatively lowcomputational price) of ADSIC.We have also explored a known collapse of the ADSICapproach for extended systems. It was shown that PZSICalso fails to cure flaws of LDA in this regime. An opti-mistic interpretation of the data obtained on the exam-ple of carbon chains allows to conclude that an efficientorbital-density dependent SIC should provide weak local-ization of the single-electron states over several atoms.Such a weak localization is in line with the excellent per-formance of bare ADSIC in case of the smaller moleculesstudied here.Whereas the results of this survey question the qual-ity of PZSIC as a benchmark approach to a SIC, theysimultaneously encourage the educated use of the muchsimpler ADSIC approach. However, as also noted, AD-SIC certainly does not provide the ultimate SIC schemeas it fails by construction, for example in the modeling ofdissociation processes or strong ionization. The limits ofADSIC with respect to dissociation or molecular struc-tural rearrangement need to be explored further. Still itremains a viable and robust option for many dynamicalsituations, especially in the case of perturbative ioniza-tion, where the ionization potential, precisely the nega-tive HOMO level, plays a central role. This implies thatADSIC remains the favorably self-interaction correctionin the calculation of reliable photo-electron spectra andangular distributions of emitted electrons, which repre-sent an ever-growing issue in the dynamics of irradiatedclusters and molecules.Future work should also aim at investigating to whichextent the level of localization can be controlled withinthe PZSIC scheme by modifying the functional form, e.g.,within the framework of GGA-SIC or by implying alter-native localization criteria during the optimization of in-ternal degrees of freedom, i.e., the unitary transformationamongst the single-particle states.0
ACKNOWLEDGMENTS
The authors acknowledge support from Institut Uni-versitaire de France. One of us (PK) also thanks theLaboratoire de Physique Theorique de Toulouse for its hospitality, the Centre National de la Recherche Scien-tifique for financial support, and P. Wopperer and S.Kl¨upfel for fruitful discussions. Allocations of compu-tational resources at the Regional Compute Center Er-langen (RRZE), Calcul en Midi-Pyr´en´ees (CALMIP) andunder the Nordic High Performance Computing (NHPC)project are gratefully acknowledged. [1] R. G. Parr and W. Yang,
Density-Functional Theory ofAtoms and Molecules (Oxford University Press, Oxford,1989).[2] R. M. Dreizler and E. K. U. Gross,
Density FunctionalTheory: An Approach to the Quantum Many-Body Prob-lem (Springer-Verlag, Berlin, 1990).[3] W. Kohn, Rev. Mod. Phys. , 1253 (1999).[4] R. G. Parr and W. Yang, Density-Functional Theoryof Atoms and Molecules (International Series of Mono-graphs on Chemistry) (Oxford University Press, USA,1994).[5] P.-G. Reinhard and E. Suraud,
Introduction to ClusterDynamics (Wiley, New York, 2004).[6] M. A. L. Marques, C. A. Ullrich, and F. Nogueira (edts.),
Time-dependent density functional theory , Lecture Notesin Physics, Vol. 706 (Springer, Berlin, 2006) p. 391.[7] T. Fennel, K.-H. Meiwes-Broer, J. Tiggesb¨aumker, P. M.Dinh, P.-G. Reinhard, and E. Suraud, Rev. Mod. Phys. , 1793 (2010).[8] J. P. Perdew and A. Zunger, Phys. Rev. B , 5048(1981).[9] H. B. Shore, J. H. Rose, and E. Zaremba, Phys. Rev. B , 2858 (1977).[10] K. Schwarz, Chem. Phys. Lett. , 605 (1978).[11] S. J. A. van Gisbergen, P. R. T. Schipper, O. V. Grit-senko, E. J. Baerends, J. G. Snijders, B. Champagne,and B. Kirtman, Phys. Rev. Lett. , 694 (1999).[12] S. K¨ummel, L. Kronik, and J. P. Perdew, Phys. Rev.Lett. , 213002 (2004).[13] M. S. Hybertsen and S. G. Louie, Phys. Rev. B , 5390(1986).[14] R. M. Nieminen, Current Opinion in Solid State and Ma-terials Science , 493 (1999).[15] M. B¨ar, L. V. Moskaleva, M. Winkler, P.-G. Reinhard,N. R¨osch, and E. Suraud, Eur. Phys. J. D , 507 (2007).[16] P. M. Dinh, P.-G. Reinhard, and E. Suraud, Phys. Rep. , 43 (2009).[17] A. Pohl, P.-G. Reinhard, and E. Suraud, Phys. Rev.Lett. , 5090 (2000).[18] A. Pohl, P.-G. Reinhard, and E. Suraud, Phys. Rev. A , 023202 (2004).[19] U. De Giovannini, D. Varsano, M. A. L. Marques, H. Ap-pel, E. K. U. Gross, and A. Rubio, Phys. Rev. A ,062515 (2012).[20] J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz,Phys. Rev. Lett. , 1691 (1982).[21] J. P. Perdew and M. Levy, Phys. Rev. Lett. , 1884(1983).[22] L. J. Sham and M. Schl¨uter, Phys. Rev. Lett. , 1888(1983).[23] A. D. Becke, Phys. Rev. A , 3098 (1988).[24] J. P. Perdew, “Unified theory of exchange and correlation beyond the local density approximation,” in ElectronicStructure of Solids 91 , edited by P. Ziesche and H. Eschrig(Akademie Verlag, Berlin, 1991) p. 11.[25] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.Lett. , 3865 (1996).[26] M. E. Casida and D. R. Salahub, J. Chem. Phys. ,8918 (2000).[27] J. P. Perdew, Chem. Phys. Lett. , 127 (1979).[28] M. R. Pederson, R. A. Heaton, and C. C. Lin, J. Chem.Phys. , 1972 (1984).[29] S. Goedecker and C. J. Umrigar, Phys. Rev. A , 1765(1997).[30] V. Polo, E. Kraka, and D. Cremer, Mol. Phys. , 1771(2002).[31] O. A. Vydrov and G. E. Scuseria, J. Chem. Phys. ,8187 (2004).[32] D. Hofmann, S. Kl¨upfel, P. Kl¨upfel, and S. K¨ummel,Phys. Rev. A , 062514 (2012).[33] S. Kl¨upfel, P. Kl¨upfel, and H. J´onsson, Phys. Rev. A ,050501 (2011).[34] S. Kl¨upfel, P. Kl¨upfel, and H. J´onsson, J. Chem. Phys. , 124102 (2012).[35] A. Svane, Phys. Rev. B , 4275 (1996).[36] J. Messud, P. M. Dinh, P.-G. Reinhard, and E. Suraud,Phys. Rev. Lett. , 096404 (2008).[37] R. T. Sharp and G. K. Horton, Phys. Rev. , 317 (1953).[38] S. K¨ummel and L. Kronik, Rev. Mod. Phys. , 3 (2008).[39] J. B. Krieger, Y. Li, and G. J. Iafrate, Phys. Rev. A ,101 (1992).[40] E. Fermi and E. Amaldi, Accad. Ital. Rome , 117 (1934).[41] C. Legrand, E. Suraud, and P.-G. Reinhard, J. Phys. B , 1115 (2002).[42] C. Legrand, E. Suraud, and P.-G. Reinhard, J. Phys. B , 1115 (2002).[43] I. Ciofini, H. Chermette, and C. Adamo, Chem. Phys.Lett. , 12 (2003).[44] I. Ciofini, C. Adamo, and H. Chermette, Chem. Phys. , 67 (2005).[45] T. Koopmans, Physica , 104 (1933).[46] P. Hohenberg and W. Kohn, Phys. Rev. , 864 (1964).[47] W. Kohn and L. J. Sham, Phys. Rev. , 1133 (1965).[48] L. A. Curtiss, K. Raghavachari, P. C. Redfern, and J. A.Pople, Journal of Chemical Physics , 1063 (1997).[49] L. A. Curtiss, K. Raghavachari, G. W. Trucks, and J. A.Pople, J.Chem.Phys. , 7221 (1991).[50] C. Hartwigsen, S. Goedecker, and J. Hutter, Phys. Rev.B , 3641 (1998).[51] Russell D. Johnson III, “NIST Computational ChemistryComparison and Benchmark Database NIST StandardReference Database Number 101 Release 14,” (2006).[52] J. Messud, P. M. Dinh, P.-G. Reinhard, and E. Suraud,Ann. Phys. (N.Y.) , 955 (2008). [53] P. E. Bl¨ochl, Phys. Rev. B , 17953 (1994).[54] J. J. Mortensen, L. B. Hansen, and K. W. Jacobsen,Phys. Rev. B , 035109 (2005).[55] C.-H. Park, A. Ferretti, I. Dabo, N. Poilvert, andN. Marzari, (2011), arXiv:1108.5726 [cond-mat.mtrl-sci].[56] A. Valdes, J. Brillet, M. Graetzel, H. Gudmundsdot-tir, H. A. Hansen, H. Jonsson, P. Kl¨upfel, G.-J. Kroes,F. Le Formal, I. C. Man, R. S. Martins, J. K. Norskov, J. Rossmeisl, K. Sivula, A. Vojvodic, and M. Zach, Phys.Chem. Chem. Phys. , 49 (2012).[57] J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson,M. R. Pederson, D. J. Singh, and C. Fiolhais, Phys. Rev.B , 6671 (1992).[58] S. J. A. van Gisbergen, F. Kootstra, P. R. T. Schipper,O. V. Gritsenko, J. G. Snijders, and E. J. Baerends,Phys. Rev. A57