Atomic self-interaction correction for molecules and solids

We present an atomic orbital based approximate scheme for self-interaction correction (SIC) to the local density approximation of density functional theory. The method, based on the idea of Filippetti and Spaldin [Phys. Rev. B 67, 125109 (2003)], is implemented in a code using localized numerical atomic orbital basis sets and is now suitable for both molecules and extended solids. After deriving the fundamental equations as a non-variational approximation of the self-consistent SIC theory, we present results for a wide range of molecules and insulators. In particular, we investigate the effect of re-scaling the self-interaction correction and we establish a link with the existing atomic-like corrective scheme LDA+U. We find that when no re-scaling is applied, i.e. when we consider the entire atomic correction, the Kohn-Sham HOMO eigenvalue is a rather good approximation to the experimental ionization potential for molecules. Similarly the HOMO eigenvalues of negatively charged molecules reproduce closely the molecular affinities. In contrast a re-scaling of about 50% is necessary to reproduce insulator bandgaps in solids, which otherwise are largely overestimated. The method therefore represents a Kohn-Sham based single-particle theory and offers good prospects for applications where the actual position of the Kohn-Sham eigenvalues is important, such as quantum transport.

We present an atomic orbital based approximate scheme for self-interaction correction (SIC) to the local density approximation of density functional theory. The method, based on the idea of Filippetti and Spaldin [ [1] Phys. Rev. B 67, 125109 (2003)], is implemented in a code using localized numerical atomic orbital basis sets and is now suitable for both molecules and extended solids. After deriving the fundamental equations as a non-variational approximation of the selfconsistent SIC theory, we present results for a wide range of molecules and insulators. In particular, we investigate the effect of re-scaling the self-interaction correction and we establish a link with the existing atomic-like corrective scheme LDA+U . We find that when no re-scaling is applied, i.e. when we consider the entire atomic correction, the Kohn-Sham HOMO eigenvalue is a rather good approximation to the experimental ionization potential for molecules. Similarly the HOMO eigenvalues of negatively charged molecules reproduce closely the molecular affinities. In contrast a re-scaling of about 50% is necessary to reproduce insulator bandgaps in solids, which otherwise are largely overestimated. The method therefore represents a Kohn-Sham based single-particle theory and offers good prospects for applications where the actual position of the Kohn-Sham eigenvalues is important, such as quantum transport.

I. INTRODUCTION
Density functional theory (DFT), in both its static [2] and time-dependent [3] forms, has become by far the most successful and widely used among all the electronic structure methods. The most obvious reason for this success is that it provides accurate predictions of numerous properties of atoms, inorganic molecules, bio-molecules, nanostructures and solids, thus serving different scientific communities.
In addition DFT has a solid theoretical foundation. The Hohenberg-Kohn theorem [2] establishes the existence of a unique energy functional E[ρ] of the electron density ρ which alone is sufficient to determine the exact ground-state of a N -electron system. Although the energy functional itself is not known, several of its general properties can be demonstrated rigorously. These are crucial for constructing increasingly more predictive approximations to the functional and for addressing the failures of such approximations [4].
Finally, but no less important, the Kohn-Sham (KS) formulation of DFT [5] establishes a one to one mapping of the intrinsically many-body problem onto a fictitious single-particle system and offers a convenient way for minimizing E[ρ]. The degree of complexity of the Kohn-Sham (KS) problem depends on the approximation chosen for the density functional. In the case of the local density approximation (LDA) [5] the KS problem typically scales as N 3 , where the scaling is dominated by the diagonalization algorithm. However, clever choices with regards to basis sets and sophisticated numerical methods make order-N scaling a reality [6,7].
The energy functional E[ρ ↑ , ρ ↓ ] (ρ σ , σ =↑, ↓ is the spin density, ρ = σ ρ σ ) can be written as where T S is the kinetic energy of the non-interacting system, v(r) the external potential, U the Hartree electrostatic energy and E xc the exchange and correlation (XC) energy. This last term is unknown and must be approximated. The construction of an approximated functional follows two strategies: empirical and "constraint satisfaction".
Empirical XC functionals usually violate some of the constraints imposed by exact DFT, and rely on parameterizations obtained by fitting representative data. One includes in this category, XC functionals which borrow some functional dependence from other theories. This is for instance the case of the celebrated LDA+U scheme [8,9], where the Hubbard-U energy takes the place of the LDA energy for certain "strongly correlated" atomic orbitals (typically d and f shells). The method however depends on the knowledge of the Coulomb and exchange parameters U and J, which vary from material to material, and can also be different for the same ion in different chemical environments [10,11].
In contrast the construction based on "constraint satisfaction" proceeds by developing increasingly more sophisticated functionals, which nevertheless satisfy most of the properties of exact DFT [12]. It was argued that this method constructs a "Jacob's ladder" [13], where functionals are assigned to different rungs depending on the number of ingredients they include. Thus the LDA, which depends only on the spin-densities is on the first rung, the generalized gradient approximation (GGA) [14], which depends also on ∇ρ σ , is on the second rung, the so-called meta-GGA functionals [15], which in addition to ρ σ and ∇ρ σ depend on either the Laplacian ∇ 2 ρ σ or the orbital kinetic energy density, are on the third rung and so on. The higher its position on the ladder the more accurate a functional becomes, but at the price of increasing computational overheads. Therefore its worth investigating corrections to the functionals of the lower rungs, which preserve most of the fundamental properties of DFT and do not generate massive additional numerical overheads.
One of the fundamental problems intrinsic to the semilocal functionals of the first three rungs is the presence of self-interaction (SI) [16]. This is the spurious interaction of an electron in a given KS orbital with the Hartree and XC potential generated by itself. Such an interaction is not present in the Hartree-Fock method, where the Coulomb self-interaction of occupied orbitals is cancelled exactly by the non-local exchange. However when using semi-local functionals such a cancellation is not complete and the rigorous condition for KS-DFT for the orbital density ρ σ n = |ψ σ n | 2 of the fully occupied KS orbital ψ σ n is not satisfied. A direct consequence of the self-interaction in LDA/GGA is that the KS potential becomes too repulsive and exhibits an incorrect asymptotic behavior [16]. This "schizophrenic" (self-interacting) nature of semilocal KS potentials generates a number of failures in describing elementary properties of atoms, molecules and solids. Negatively charged ions (H − , O − , F − ) and molecules are predicted to be unstable by LDA [17] and transition metal oxides are described as small-gap Mott-Hubbard antiferromagnets (MnO, NiO) [18] or even as ferromagnetic metals (FeO, CoO) [18] instead of chargetransfer insulators. Moreover the KS highest occupied molecular orbital (HOMO), the only KS eigenvalue that can be rigorously associated to a single particle energy [19,20,21], is usually nowhere near the actual ionization potential [16].
Finally XC functionals affected by SI do not present a derivative discontinuity as a function of the occupation [19,20]. Semi-local functionals in fact continuously connect the orbital levels of systems of different integer occupation. This means for instance that when adding an extra electron to an open shell N -electron system the KS potential does not jump discontinuously by I N − A N where I N and A N are respectively the ionization potential and the electron affinity for the N -electron system. This serious drawback is responsible for the incorrect dissociation of heteronuclear molecules into charged ions [22] and for the metallic conductance of insulating molecules [23].
The problem of removing the SI from a semi-local potential was acknowledged a long time ago when Fermi and Amaldi proposed a first crude correction [24]. However the modern theory of self-interaction corrections (SIC) in DFT is due to the original work of Perdew and Zunger from almost three decades ago [16]. Their idea consists in removing directly the self-Hartree and self-XC energy of all the occupied KS orbitals from the LDA XC functional (the same argument is valid for other semi-local functionals), thus defining the SIC functional as (3) Although apparently simple, the SIC method is more involved than standard KS DFT. The theory is still a density functional one, i.e. it satisfies the Hohenberg-Kohn theorem, however it does not fit into the Kohn-Sham scheme, since the one-particle potential is orbitaldependent. This means that one cannot define a kinetic energy functional independently from the choice of E xc [16]. Two immediate consequences are that the ψ σ n are not orthogonal, and that the orbital-dependent potential can break the symmetry of the system. This last aspect is particularly important for solids since one has to give up the Bloch representation.
In this paper we explore an approximate method for SIC to the LDA, which has the benefit of preserving the local nature of the LDA potential, and therefore maintains all of the system's symmetries. We have followed in the footsteps of Filippetti and Spaldin [1], who extended the original idea of Vogel and co-workers [25,26,27] of considering only the atomic contributions to the SIC. We have implemented such a scheme into the localized atomic orbital code Siesta [28] and applied it to a vast range of molecules and solids. In particular we have investigated in detail how the scheme performs as a singleparticle theory, and how the SIC should be rescaled in different chemical environments.

II. REVIEW OF EXISTING METHODS
The direct subtraction proposed by Perdew and Zunger is the foundation of the modern SIC method. However the minimization of the SIC functional (3) is not trivial, in particular for extended systems. The main problem is that E xc itself depends on the KS orbitals. Thus it does not fit into the standard KS scheme and a more complicated minimization procedure is needed.
Following the minimization strategy proposed by Levy [29], which prescribes to minimize the functional first with respect to the KS orbitals ψ σ n and then with respect to the occupation numbers p σ n , Perdew and Zunger derived a set of single-particle equations where the effective single-particle potential v σ eff,n (r) is de- and These are solved in the standard KS way for atoms, with good results in terms of quasi-particle energies [16]. In this particular case the KS orbitals ψ σ n show only small deviations from orthogonality, which is re-imposed with a standard Schmidt orthogonalization.
The problem of the non-orthogonality of the KS orbitals can be easily solved by imposing the orthogonality condition as a constraint to the energy functional, thus obtaining the following single-particle equation Even in this case where orthogonality is imposed, two major problems remain: the orbitals minimizing the energy functional are not KS-type and in general do not satisfy the system's symmetries. If one insists in minimizing the energy functional in a KS fashion by constructing the orbitals according to the symmetries of the system, the theory will become sizeinconsistent, or in other words it will be dependent on the particular representation employed. Thus one might arrive at a paradox, where in the self-interaction of N hydrogen atoms arranged on a regular lattice of large lattice spacing (in such a way that there is no interaction between the atoms) vanishes, since the SIC of a Bloch state vanishes for N → ∞. However the SIC for an individual H atom, when calculated using atomic-like orbitals, accounts for essentially all the ground-state energy error [16]. Therefore a size-consistent theory of SIC DFT must look for a scheme where a unitary transformation of the occupied orbitals, which minimizes the SIC energy is performed. This idea is at the foundation of all modern implementations of SIC.
Significant progress towards the construction of a sizeconsistent SIC theory was made by Pederson, Heaton and Lin, who introduced two sets of orbitals: localized orbitals φ σ n minimizing E SIC xc and canonical (Kohn-Sham) de-localized orbitals ψ σ n [30,31,32]. The localized orbitals are used for defining the densities entering into the effective potential (5), while the canonical orbitals are used for extracting the Lagrangian multipliers ǫ σ,SIC nm , which are then associated to the KS eigenvalues. The two sets are related by unitary transformation ψ σ n = m M σ nm φ σ m , and one has two possible strategies for minimizing the total energy. The first consists in a direct minimization with respect to the localized orbitals φ σ n , i.e. in solving equation (8) when we replace ψ with φ and the orbital densities entering the definition of the one-particle potential (5) are simply ρ σ n = |φ σ n | 2 . In addition the following minimization condition must be satisfied where v SIC . An expression for the gradient of the SIC functional, which also constrains the orbitals to be orthogonal has been derived [33] and applied to atoms and molecules with a mixture of successes and bad failures [34,35,36].
The second strategy uses the canonical orbitals ψ and seeks the minimization of the SIC energy by varying both the orbitals ψ and the unitary transformation M. The corresponding set of equations is where H σ 0 is the standard LDA Hamiltonian (without SIC). Thus the SIC potential for the canonical orbitals appears as a weighted average of the SIC potential for the localized orbitals. The solutions of the set of equations (10) is somehow more appealing than that associated to the localized orbitals since the canonical orbitals can be constructed in a way that preserves the system's symmetries (for instance translational invariance).
A convenient way for solving the equation (10) is that of using the so called "unified Hamiltonian" method [30]. This is defined as (we drop the spin index σ) whereP n = |ψ σ n ψ σ n | is the projector over the occupied orbital ψ σ n , andQ is the projector over the unoccupied onesQ = 1 − occup nP n . The crucial point is that the diagonal elements of the matrix ǫ σ,SIC nm and their corresponding orbitals ψ σ n are respectively eigenvalues and eigenvectors of H u , from which the whole ǫ σ,SIC nm can be constructed. Finally, and perhaps most important, at the minimum of the SIC functional, the canonical orbitals diagonalize the matrix ǫ σ,SIC nm , whose eigenvalues can now be interpreted as an analogue of the Kohn-Sham eigenvalues [32].
It is also interesting to note that an alternative way for obtaining orbital energies is that of constructing an effective SI-free local potential using the Krieger-Li-Iafrate method [37]. This has been recently explored by several groups [38,39,40] When applied to extended systems the SIC method demands considerable additional computational overheads over standard LDA. Thus for a long time it has not encountered the favor of the general solid state community. In the case of solids the price to pay for not using canonical orbitals is enormous since the Bloch representation should be abandoned and in principle infinite cells should be considered. For this reason the second minimization scheme, in which the canonical orbitals are in a Bloch form, is more suitable. In this case for each k-vector one can derive an equation identical to equation (10), where ǫ σ,SIC nm = ǫ σ,SIC nm (k) and n is simply the band index [41]. The associated localized orbitals φ for instance can be constructed as Wannier functions and the minimization scheme proceeds in a similar way to that done for molecules.
The problem here is that in practice, the cell needed to describe the localized states φ may be considerably larger than the primitive unit cell. This is not the case for ionic insulators [41], where the localized orbitals are well approximated by atomic orbitals. Such a simplification is however not valid in general. For example supercells as large as 500 atoms have been considered for describing the localized d shells in transition metals oxides [42,43,44]. Despite these difficulties the SIC scheme has been applied to a vast range of solid state systems with systematic improvement over LDA. These include, in addition to transition metals monoxides [42,44,45], rare-earth materials [46], diluted magnetic semiconductors [47], Fe 3 O 4 [48], heavy elements compounds [49], just to name a few.
In order to make the SIC method more scalable several approximations have been proposed. One possibility is that of incorporating part of the SIC into the definition of the pseudopotentials [50]. The idea consists in subtracting the atomic SI from the free atom, and then transferring the resulting electronic structure to the definition of a standard norm-conserving pseudopotential. This approximation is sustained by the fact that the transformation M, which relates canonical and localized orbitals does hardly mix core and valence states [32]. Thus the SIC contribution to the total energy can be separated into the contributions from the core and the valence and in first approximation the latter can be neglected [51]. The benefit of this method is that translational invariance is regained and the complicated procedure of minimizing M is replaced by a pseudopotential calculation.
A further improvement over the pseudopotential method was recently presented by Vogel and co-workers [25,26,27] and then extended by Filippetti and Spaldin [1]. The method still assumes separability between the core and the valence contributions to the SIC, but now the SIC for the valence electrons is approximated by an atomic-like contribution, instead of being neglected. This atomic SIC (ASIC) scheme is certainly a drastic approx-imation, since it implicitly assumes that the transformation M minimizing the SIC functional leads to atomic like orbitals.
In the work of Vogel this additional contribution is not evaluated self-consistently for the solid, while the implementation of Filippetti assumes a linear dependance of the SIC over the orbital occupation. In spite of the approximations involved, the method has been applied successfully to a range of solids including II-VI semiconductors and nitrites [1,25,26], transition metal monoxides [1,52], silver halides [27], noble metal oxides [53], ferroelectric materials [1,54,55], high-k materials [56] and diluted magnetic semiconductors [57,58]. Interestingly most of the systems addressed by the ASIC method are characterized by semi-core d orbitals, for which an atomic correction looks appropriate, and a similar argument is probably valid for ionic compounds as recently demonstrated for the case of SiC [59].
Here we further investigate the self-consistent ASIC method [1] by examining both finite and extended systems, and by critically considering whether a scaling factor, additional to the orbital occupation, is needed for reproducing the correct single particle spectrum.

III. FORMALISM AND IMPLEMENTATION
In this section we derive the fundamental equations of the ASIC method, while looking closely at the main approximations involved in comparison to the fully selfconsistent SIC approach. Our practical implementation is also described.

A. The ASIC potential
The starting point of our analysis is the SIC Schrödinger-like equation (10) for the canonical orbitals. Let us assume, as from reference [51], that the rotation M transforming localized orbitals (to be determined) into canonical orbitals (see equation (11)) does not mix core and valence states. We also assume that core electrons are well localized into atomic-like wave-functions and that they can be effectively described by a normconserving pseudopotential.
Let us now assume that M is known and so are the localized orbitals φ σ m . In this case the canonical orbitals diagonalize ǫ σ,SIC nm and the equation (10) simply reduces to with ∆v SIC n defined in equation (12). The Hamiltonian H σ 0 + ∆v SIC n can be then re-written in a convenient form as where v SIC m is the self-interaction potential for the localized orbital φ σ m , andP φ m is the projector over the same statê Two main approximations are then taken in the ASIC approach [1,25]. First the localized states φ σ m are assumed to be atomic-like orbitals Φ σ m (ASIC orbitals) and the SIC potential is approximated as withṽ σSIC (16) obtained by replacing the φ's with the ASIC orbitals Φ, and α is a scaling factor. Note that the orbitals Φ m are not explicitly spin-dependent and one simply has Φ σ m =Φ m p σ m with p σ m the orbital occupation (p σ m = 0, 1). The factor α is an empirical factor, which accounts for the particular choice of ASIC orbitals. This first approximation is expected to be accurate for systems retaining an atomic-like charge density as in the case of small molecules. It is also formally exact in the one-electron limit (for α = 1). In the case of extended solids the situation is less transparent, since in general the functions minimizing E SIC xc are Wannier-like functions [60].
The second approximation taken in the ASIC method is that of replacing the non-local projectorP Φ m with its expectation value. The idea is that the SIC potential for the canonical orbitals ∆v SIC n is formally a weighted average of the SIC potential for the localized orbitals v SIC m . For the exact SIC method the weighting factor is the non-local projector M nm φ σ m ψ σ n . This means that the SIC potential for a given canonical orbital ψ n is maximized in those regions where the overlap with some of the localized orbitals φ m is maximum. In the ASIC method such nonlocal projector is replaced more conveniently by a scalar. In the original proposal by Voegl et al. [25,26,27] this was simply set to one. Here we consider the orbital occupation p σ m of the given ASIC orbital Φ m , i.e. we replacê P Φ m with its expectation valuê where f σ n is the occupation number of the Kohn-Sham orbital ψ σ n . The final form of ASIC potential is then Let us now comment on the empirical scaling factor α. In reference [1] α was set to 1/2 in order to capture eigenvalue relaxation. This choice however violates the one-electron limit of the SIC potential, which is correctly reproduced for α = 1. We can then interpret α as a measure of the deviation of the ASIC potential from the exact SIC potential. Ultimately α reflects the deviation of the actual ASIC projectors |Φ Φ| from the localized orbitals defining the SI corrected ground state. One then expects α to be close to 1 for systems with an atomic-like charge density, and to vanish for metals, whose valence charge density resembles that of a uniform electron gas [61]. Intermediate values are then expected for situations different from these two extremes, and we will show that a values around 1/2 describe well a vast class of mid-to wide-gap insulators.

B. Implementation
The final form of the SIC potential to subtract from the LSDA (local spin density approximation) one (equation (19)) is that of a linear combination of non-local pseudopotential-like terms. These are uniquely defined by the choice of exchange and correlation potential used and by ASIC orbitals Φ m . The practical way of constructing such potentials, i.e. the way of importing the atomic SIC to the solid state, depends on the specific numerical implementation used for the DFT algorithm. At present only plane-wave implementations are available [1,25,26,27], while here we present our new scheme based on the pseudo atomic orbital (PAO) [65] code Siesta [28].
We start by solving the atomic all-electron SIC-LSDA equation for all the species involved in the solid state calculation. Here we apply the original Perdew-Zunger (PZ-SIC) formalism [16] and we neglect the small nonorthogonality between the Kohn-Sham orbitals. Thus we obtain a set of SI corrected atomic orbitals Φ m , which exactly solve the atomic SIC-LSDA problem. The atomic orbitals Φ m describing the valence electrons are then used to define the ASIC potentialsṽ σSIC with ρ σ m = |Φ m | 2 . At the same time a standard LSDA calculation for the same atoms is used to construct the pseudopotentials describing the core electrons. These are standard normconserving scalar relativistic Troullier-Martins pseudopotentials [62] with nonlinear core corrections [63]. Thus we usually neglect the SIC over the core states, when constructing the pseudopotentials. This is justified by the fact that the eigenvalues for the SIC-LSDA-pseudoatom, i.e. for the free atom where the effects of core electrons are replaced by LSDA pseudopotentials but SIC is applied to the valence electrons, are in excellent agreement with those obtained by all-electron SIC-LSDA calculations [25].
The final step is that of recasting the ASIC potentials v SIC m (r), which have a −2/r asymptotic behaviour, in a suitable non-local form. This is obtained with the standard Kleinman-Bylander [64] scheme and the final ASIC potential (equation (19)) is written as where the ASIC projectors are given by and the normalization factors are The orbital functions Φ ′ m are atomic-like functions with a finite range, which ensure that the ASIC projectors γ m do not extend beyond that range. These are constructed in the same way as the Siesta basis set orbitals, i.e. as solutions of the pseudo-atomic problem with an additional confining potential at the cutoff radius r cutoff [65]. The choice of the appropriate cutoff radius for the SIC projectors should take into account the two following requirements. On the one hand it should be sufficiently large to capture most of the SIC corrections. A good criterion [1] is that the SIC-LSDA contribution to the orbital energy of the free atom is reproduced within some tolerance. On the other hand the cutoff should be reasonably short so as not to change the connectivity of the matrix elements of the PAO Hamiltonian. In other words we need to ensure that orbitals otherwise considered as disconnected in evaluating the various parts of the Hamiltonian matrix are not considered connected for the v σ ASIC matrix elements alone. As a practical rule we set the cutoff radius for a particular orbital of a given atom to be either equal to the longest among the cutoff radii of the PAO basis set for that particular atom (typically the first ζ of the lowest angular momentum), or, if shorter, the radius at which δε σ SIC m < 0.1mRy. Typically, when reasonable cutoff radii (6 to 9 Bohr) are used, we find that the atomic SIC-LSDA eigenvalues are reproduced to within 1 to 5 mRy for the most extended shells and to within 0.1 mRy for more confined shells. Thus δε σ SIC m are rather well converged already for cutoff radii defined by a PAO energies shifts [28] of around 20mRy, although usually smaller PAO energy shifts are necessary for highly converged total energy calculations.
Finally the matrix elements of the SIC potential are calculated as usual over the Siesta basis set. Additional basis functions χ m are constructed from the confined localized atomic orbitals described before using the splitnorm scheme [65,66,67]. The density matrix ρ σ is represented over such basis ρ σ µν and the orbital populations are calculated as where S mµ is a matrix element of the overlap matrix. Note that in principle the orbital population should be constructed for the atomic SIC orbital Φ m . However, we notice that p σ m is rather insensitive to the specific choice of orbital, once this has a reasonable radial range. For practical numerical reasons in the present implementation, we always use the orbital populations projected onto the basis set sub-space consisting of the most extended first-ζ orbitals of the atomic species involved. The matrix elements of the SIC potential are simply and the ASIC-KS equation takes the final form with v PP the pseudopotential.

C. Total Energy
The energy corresponding to the SIC-LSDA functional is given by [16] with E LSDA xc the LSDA exchange and correlation energy density. The orbital densities entering in the SI term are those associated to the local orbitals φ. As we have already mentioned, this functional needs to be minimized with respect to the φ's, which are an implicit function of the total spin density ρ σ . In the ASIC approximation these orbitals are not minimized, but taken as atomic functions. This means that in the present form the theory is not variational, in the sense that there is no functional related to the KS equation (27) by a variational principle. With this in mind we adopt the expression of equation (28) as a suitable energy, where the orbital densities are those given by the ASIC orbitals In our implementation the LSDA KS energy E LSDA is directly available as calculated in the Siesta code [28] and thus, only the second term of equation (28) needs to be calculated. This is easily done by calculating both U and E LSDA xc on an atomic radial grid for each atomic orbital in the system.

D. ASIC and LDA+U
We now compare our ASIC method with another atomic-like correction to LSDA, namely the LDA+U method [8,9]. In LDA+U one replaces the LSDA exchange and correlation energy associated to the "correlated" orbitals (d or f shells), with the Hubbard-U energy. Thus the functional becomes where the Hubbard energy E U and the double counting term E DC depend on the orbital populations p σ m of the correlated orbitals. Several forms for the LDA+U functional have been proposed to date. A particularly simple and transparent one [10,68], which is also rotationally invariant, redefines the U parameter as an effective parameter U eff = U − J and the functional takes the form where in we separate out the index for the atomic position I from the magnetic quantum number m, and introduce the off-diagonal populations Note that the LDA+U functional is SI free for those orbitals that are corrected.
Although a rotationally invariant form of the ASIC potential can be easily derived, we assume here for simplicity that the system under consideration is rotationally invariant, or alternatively that we have carried out a rotation, which diagonalizes the p I σ mn matrix. In this case the energy becomes simply with p I σ m = p I σ mm . It is then easy to compute the KS potential and the orbital energy These need to be compared with the ASIC potential (equation (19)) and orbital energy where the last term follows from ∂E ∂p Iσ m = C Iσ m and from equation (28). The main difference between the ASIC and LDA+U method is in the way in which unoccupied states are handled. In fact, while LDA+U corrects unoccupied states and pushes the orbital energies upwards by ∼ U eff /2, ASIC operates only on occupied states, that are shifted towards lower energies by C I σ m . This reflects the fact that the SIC is defined only for occupied KS orbitals. An important consequence is that the opening of bandgaps in the electronic structures, one of the main features of both the LDA+U and ASIC schemes, is then driven by two different mechanisms. On the one hand in LDA+U , gaps open up since occupied and unoccupied states are corrected in opposite directions leading to a gap of ∼ U eff . On the other hand ASIC is active only over occupied states and gaps open only if occupied and unoccupied bands have large differences in their projected atomic orbital content. Thus one should not expect any corrections for covalent materials where conduction and valence bands are simply bonding and antibonding states formed by the same atomic orbitals. This is for instance the case of Si and Ge. In contrast ASIC will be extremely effective for more ionic situations, where the orbital contents of conduction and valence bands are different.
Finally, by comparing the corrections to the orbital energy of a fully occupied state, one finds which establishes an empirical relation between the Hubbard energy and the ASIC correction. Since U is sensitive to the chemical environment due to screening [10], while all the other quantities are uniquely defined by an atomic calculation, we can re-interpret the parameter α as empirically describing the screening from the chemical environment within the ASIC scheme.

IV. RESULTS: EXTENDED SYSTEMS
The test calculations that we present in this work are for two classes of materials: extended and finite. First we investigate how our implementation performs in the solid state. In particular we discuss the rôle of the parameter α in determining the bandstructure of several semiconductors, considering both the KS band-gap and the position of bands associated with tightly bound electrons.

A. Estimate of α for semiconductors
The quasi-particle band gap E g in a semiconductor is defined as the difference between its ionization potential I and electron affinity A. These can be rigorously calculated from DFT as the HOMO energy respectively of the neutral and negatively charged systems. This actual gap cannot be directly compared with the KS band-gap E KS g , defined as the difference between the orbital energy of the HOMO and LUMO states of the neutral system. In fact, the presence of a derivative discontinuity in the DFT energy as a function of the electron occupation establishes the following rigorous relation [20,69] with This is valid even for the exact XC potential, and therefore in principle one has to give up KS bandstructures as a tool for evaluating semiconductor band-gaps. The size of ∆ xc is however not known for real extended systems and the question of whether most of the error in determining E g from E KS g is due to the approximation in the XC potential or due to the intrinsic ∆ xc is a matter of debate.
In general SI-free potentials bind more than LSDA and one expects larger gaps. Surprisingly, functionals based on exact exchange, provide KS gaps rather close to the experimental values [70,71]. The reason for such a good agreement is not fully understood, but it is believed that the exact KS gaps should be smaller than the actual ones.
With this in mind, we adopt a heuristic approach and we use the KS band-gap as a quality indicator for interpreting the parameter α and for providing its numerical value for different classes of solids. Here we investigate the dependence of E KS g over α and we determine the value for α yielding the experimental band-gap. Assuming that ∆ xc does not vary considerably across the materials investigated, this will allow us to relate α to the degree of localization in a semiconductor and to extract the value useful for ASIC to be an accurate single-particle theory.
In figure 1 we present the band-gap of four representative semiconductors as a function of α together with the value needed to reproduce the experimental band-gap. LSDA corresponds to α = 0 and while α = 1 accounts for the full ASIC. In general E KS g increases as α increases, as  a result of the stronger SIC. The E KS g (α) curve is almost linear with a slope, which appears to be material-specific.
For the most ionic compound, NaCl, the experimental gap is reproduced almost exactly by α = 1, i.e. by the full ASIC. This is somehow expected since the charge density of solid NaCl is rather close to a superposition of the Na + and Cl − ionic charge densities. In this case of strongly localized charge densities the ASIC approximation is rather accurate yielding results substantially identical to those obtained with full self-consistent PZ-SIC [41]. Indeed earlier calculations for LiCl [41] demonstrate that the SIC band-structure is rather insensitive of the localized orbitals φ once these have an atomic-like form.
For the other compounds the localized orbitals φ's are not necessarily atomic-like functions and deviations from α = 1 are expected. Interestingly we find that, for all the materials investigated, a value of around 1/2 reproduces the experimental band-gap rather accurately. As an illustration, in table I we compare the experimental band-gap E exp g to the calculated E KS g for ASIC (α = 1) and LDA, for several semiconductors ranging from ionic salts to wide-gap II-VI and III-V semiconductors. We also report the value of α = α * needed for E exp g = E KS g . Clearly for all the strongly ionic compounds (LiCl, NaCl and KCl) the full ASIC correction α = 1 reproduces quite accurately the experimental gap and agrees with previous self-consistent SIC calculations [81]. For all the other compounds a value of around 1/2 is always adequate, confirming the initial choice of Filippetti and Spaldin. For these materilas we do not find any particular regularity. In general α is large when the experimental gap is large, however there is no direct connection between α and the ionicity or covalency of a compound. In fact, the improvement of the band-gap is not simply due to a rigid shift of the valence band, but usually cor-responds to a general improvement of the whole quasiparticle spectrum. Examples for ZnO and GaN will be presented in the next section.
As a further proof of this point in table II we present the valence band-width for the semiconductors investigated as calculated from LSDA ∆E LSDA v and ASIC for both α = 1 (∆E ASIC1 v ) and α = α * (∆E ASIC α * v ). We also report the experimental values ∆E exp v whenever available, although a direct comparison with experiments is difficult, since these values are rather imprecise and sometimes not known. The general feature is that ASIC produces only minor corrections over LSDA, and that these corrections do not follow a generic trend. Thus, while for the nitrites ASIC always increases the band-width, it does just the opposite for KCl, SrO and CaO. Having established α = 1/2 as an appropriate value for II-VI and III-V semiconductors, we now look at the whole band-structure (not just the fundamental gap) for a few test cases. Here we consider ZnO and GaN for which photo-emission data disagree quite remarkably from LSDA calculations. In figure 2 we compare the band structure of wurtzite ZnO obtained respectively from LSDA and our ASIC.
In ZnO, the valence band top (VBT) is predominantly oxygen 2p in character and the conduction band minimum (CBM) is essentially zinc 4s. With a value of ∼ 0.5 for the scaling parameter α, the ASIC band gap closely matches the experimental gap of E g =3.43 eV, whereas the LDA band-gap is very small (∼ 0.85 eV). Some part of the LDA band gap error in ZnO can be traced to an underestimation of the semi-core Zn 3d states. The LDA binding energy for the Zn 3d states is ∼ 5.5eV while A LM Γ A HK Γ photoemission results place them at around ∼ 7.8 eV. ASIC however rectifies the problem and is in very good agreement with experiment. This results furthermore in the removal of the spurious Zn 3d -O 2p band mixing seen in LDA. An additional feature is that the band-width of the valence band increases considerably as an effect of the downshift of the d manifold. Its worth mentioning that the positions of the Zn 3d levels obtained from ASIC in the case of ZnS, ZnSe, and ZnTe also agree remarkably well with experiment.
The wide-gap III-V semiconductor GaN presents similar phenomenology to that of ZnO. Figure 3 compares the band structure for wurtzite GaN obtained from LSDA and ASIC. When compared to X-ray photoemission spectra [89], the LSDA band structure of GaN has several shortcomings. Firstly, the band-gap between N 2p bands (VBT) and Ga 4s bands (CBM) is underestimated at around 2.2 eV against the experimental value of 3.4 eV. Secondly, the 3d states of Ga are too shallow in LSDA, leading to a spurious 3d-2s hybridization. As a result the Ga 3d states overlap with and split the N 2s bands. ASIC rectifies the picture on both counts by improving the band gap and lowering the position of the Ga 3d bands with respect to the N 2s bands. Note that these are for the rhombohedral unitcell with 4 atoms per cell. The two Mn ions are antiferromagnetically aligned and the oxygen ions are nonmagnetic. This results in a layered ferromagnetic order of the (111) planes, which in turn are anti-ferromagnetic coupled to each other. Also in this case, ASIC is a considerable improvement over LSDA. The size of the fun-damental gap now resembles the experimental one and the VBT recovers some p character.

A. Ionization potentials
In view of the fact that the ASIC method gives improved eigenvalue spectra for several solid state systems, it is worth taking a cautious look at how it performs with molecules. This is particularly important for assessing whether the ASIC scheme can be adapted to work in DFT electron transport schemes based on the KS spectra [23,90]. In exact KS DFT only the highest occupied orbital eigenvalue (ǫ HOMO ) has a rigorous physical interpretation and corresponds to the negative of the first ionization potential [19,20]. More generally, for a N electron system, the following equations hold in exact KS-DFT where −I N and −A N are the ionization potential (IP) and the electron affinity (EA) respectively. Therefore we start our analysis by looking at these quantities as calculated by ASIC. Also in this case we investigate different values of α. However here we limit ourselves only to α = 1 (ASIC 1 ) and α = 1/2 (ASIC 1/2 ). In table III and figure 5 we compare the experimental negative IP for several molecules with the corresponding (ǫ HOMO ) obtained using LSDA and ASIC. It is clear that LSDA largely underestimates the removal energies in all the cases and that the values obtained from ASIC 1/2 are also consistently lower than the experimental value. However, as made evident by the figure the agreement between ASIC 1 and experiments is surprisingly good. In fact the mean deviation δ(X) (X = LSDA, ASIC 1/2 , ASIC 1 ) from experiment  table III). It is worth noting that we have excellent agreement over the whole range of molecules investigated going from N 2 to large fullerenes C 60 and C 70 .
For comparison in figure 5 we have also included results obtained with a full self-consistent PZ-SIC approach [35]. Surprisingly our atomic approximation seems to produce a better agreement with experiments than the self-consistent scheme, which generally overcorrects the energy levels. This is a rather general feature of the PZ-SIC scheme and it is generally acknowledged that some re-scaling procedure is needed [92,93].  [91], while the star symbol represents full PZ-SIC calculations from reference [35].

B. Electron affinities
In Hartree Fock theory where Koopmans' theorem holds [94], the lowest unoccupied molecular orbital (LUMO) energy (ǫ LUMO ), corresponds to the vertical EA of the N electron system, if one neglects electronic relaxation. No such interpretation exists for (ǫ LUMO ) in DFT and so the EA is not directly accessible from the ground state spectrum of the N electron system. However, as equation (42) indicates, the EA is in principle accessi-ble from the ground state spectrum of the N + 1 − f (0 < f < 1) electron system and asserts in particular that it must be relaxation free through non-integer occupation. Unfortunately, the LSDA/GGA approximate functionals usually perform rather poorly in this regard as the N + 1 electron state is unbound with a positive eigenvalue. So one resorts instead to extracting electron affinities from more accurate total energy differences [95], or by extrapolating them from LSDA calculations for the N electron system [96]. This failing of approximate functionals has been traced in most part to the SI error and so SIC schemes are expected to be more successful in describing the N + 1 electron spectrum.
In table IV we compare HOMO energies (denoted as ǫ HOMO N +1 ) of several singly negatively charged molecules with the experimental electron affinities. We also report the LUMO energies for the corresponding neutral species (denoted as ǫ LUMO N ). LSDA relaxed geometries for the neutral molecule are used for both the neutral and charged cases. We find that various ǫ HOMO N+1 obtained from ASIC 1 once again are in reasonably good agreement with corresponding experimental electron affinities while LSDA and ASIC 1/2 continue to be poor even in this regard. In this case δ(X) stands at 4.1 eV , 2.31 eV and 0.83 eV for LSDA, ASIC 1/2 and ASIC 1 respectively. Notice that ǫ HOMO N+1 from LSDA is positive in most cases as the states are unbound.
In figure 6 we present our data together with ǫ HOMO N+1 as calculated using the PZ-SIC [35]. Again ASIC 1 performs better than PZ-SIC, that also for the EA systematically overcorrects.

C. Vertical excitations
Having shown that ASIC offers a good description of both IP and EA for a broad range of molecules, we turn  our attention to the remaining vertical ionization potentials. As mentioned before, KS-DFT lacks of Koopmans theorem, and therefore the KS energies are not expected to be close to the negative of the vertical ionization potentials. However, at least for atoms, the introduction of SIC brings a remarkable cancellation between the negative relaxation energy and the positive non-Koopmans corrections [16]. For this reason the SIC KS eigenvalues are a good approximation to the relaxed excitation energies. As an example, in table V we present the orbital energies calculated with ASIC 1 and ASIC 1/2 for the N 2 molecule. These are compared with experimental data [97] and orbital energies obtained respectively with Hartree-Fock (HF), self-consistent SIC, and SIC where molecular orbitals are used instead of localized orbitals (D-SIC) [31].
Remarkably ASIC 1 seems to offer good agreement over the whole spectrum, improving considerably over LSDA and in some cases even over SIC and HF results. This improvement is not just quantitative, but also qualitative. For instance while rectifying the LSDA spectrum of the N 2 molecule, ASIC 1 preserves the correct order between 3σ g and 1π u orbitals, which are erroneously inverted by both SIC and HF. So why does ASIC perform better than the other methods with regards to removal energies? In LSDA, electron relaxation typically cancels only half of the non-Koopmans contributions, resulting in energies that are too shallow [16]. In contrast HF lacks energy relaxation and the orbital energies are too deep. The reason why ASIC 1 performs better than self-consistent SIC is less clear. As a general consideration, also for the case of vertical ionization energies self-consistent SIC seems to overcorrect the actual values. Thus the SIC potential appears too deep, and the averaging procedure behind the ASIC approximation is likely to make it more shallow.
As a further test we calculated the orbital energies for a few other molecules and compared them both with LSDA and experiment [98]. These are presented in table VI. Again the ASIC 1 results compare rather well with experiment, and we can conclude that the ASIC method offers a rather efficient and inexpensive theory for single particle vertical excitations.
D. HOMO-LUMO gap and discontinuity of the exchange and correlation potential We are now in a position to discuss the HOMO-LUMO gap in ASIC. As already mentioned, even for the ex-act XC functional, the KS gap E KS g = ǫ LUMO − ǫ HOMO does not account for the actual quasi-particle gap E g = I N − A N . This in turn is the sum of E KS g and the discontinuity of the exchange and correlation potential ∆ xc . Equivalently i.e. ∆ xc is the discontinuity in the eigenvalue of the LUMO state at N . Therefore, in order to extract the actual gap from the KS gap, provided that the spectrum is reasonably well described at integer electron numbers N , what remains is to model the derivative discontinuity at N and ensure that ǫ HOMO N ±f is relaxation free for (0 < f < 1). Local and semi-local (LSDA/GGA) XC functionals lack such a discontinuity, while self-interaction corrections are able to restore it, at least in part. For instance the PZ-SIC scheme is successful in this regard. In figure 7 we illustrate the ionization curve for the ethylene (C 2 H 4 ) molecule as the occupation of the HOMO state is varied from 0 to 1 in going from the ionized C 2 H + 4 to the netural C 2 H 4 configuration. It is seen that among the three schemes presented, only the PZ-SIC scheme approximately models the behaviour required by the equation (41). The ASIC HOMO eigenvalue roughly agrees with the PZ-SIC eigenvalue at integer occupation but behaves linearly through non-integer values. Thus we find that the derivative discontinuity for the molecule is smoothed out in ASIC, which still connects continuously different integer occupations. This is one of the limitations of the atomic representation employed in ASIC.
In view of the foregoing discussion, the actual size of the HOMO-LUMO gap in ASIC becomes significant with a direct bearing on the physics described. Ideally, we want ǫ LUMO (N ) (LUMO for the N -electron system) to be as close to ǫ HOMO (N + 1) so that the range of eigenvalue relaxation through fractional occupation numbers M ∈ (N, N + 1) is minimized. Looking at columns 6,7 and 8 in table IV however, we see that for almost all the molecules, this is hardly the case. The agreement between ǫ LUMO (N ) and -EA from experiment (≃ ǫ HOMO (N + 1)) for ASIC 1 is quite poor implying a considerable energy range spanning fractional particle number. We still expect this energy range to be smaller for ASIC 1 than LSDA. It is also apparent from the table IV that ǫ LUMO ASIC (N ) usually differs from ǫ LUMO LSDA (N ) and in fact by considerable magnitudes in some cases. Thus the ASIC 1 "correction" to the empty LUMO state does not vanish in contrast to the PZ-SIC scheme where, by definition, the empty eigenstates are SIC free.
Since the SIC operator v σ ASIC is constructed in an atomic orbital representation, the correction to any KS eigenstate ψ σ n either filled or empty is not necessarily zero unless ψ σ n only projects onto empty atomic orbitals. Also this correction to the LUMO with respect to the LSDA is negative in most cases, exceptions being NH 2 and CH 3 where it is desirably positive. Thus the fundamental HOMO-LUMO gap in ASIC is a combination of both the HOMO and LUMO corrections. Table  VII shows how this combination works out in ASIC 1/2 and ASIC 1 when compared to LSDA. The molecular test set is the same as that in table III. We see in almost all cases the ASIC gap is systematically larger than the LSDA one. This is expected because the correction to the HOMO is usually much stronger than that to the LUMO. In general, ASIC is expected to work well for systems where the occupied and un-occupied KS eigenstates of the extended system have markedly different atomic orbital signatures being derived predominantly from filled and empty atomic orbitals respectively. In such a case, the ASIC correction to the empty states would be nullified in being scaled by near-zero atomic orbital populations. In some cases, provided phase factors combine suitably, the correction to the empty states can even be positive with respect to the same in LSDA.

E. Final Remarks
Before we conclude, we discuss some general properties of the ASIC method which are relevant to any orbital dependent SIC implementation and also some possible pitfalls. As with other SIC schemes, ASIC is not invariant under unitary transformations of the orbitals used in constructing the SIC potential. Thus the ASIC correction is likely to change as the atomic orbitals used for projecting onto the KS eigenstates of the system are rotated or transformed otherwise. Unlike the Perdew-Zunger method however, there can be no variational principle over all possible unitary transformations of the atomic orbitals because in the general case they do not represent the Hamiltonian of the system under consideration. This also implies that if the scheme is used with a system that is already well described by LSDA, the "correction" additional to the LSDA result does not necessarily vanish. Simple metals and narrow gap systems are likely candidates for this scenario.
Furthermore, its pertinent to mention that ASIC becomes ineffective if not counterproductive for materials with homonuclear bonding, in which valence and conduction bands have the same atomic orbital character. In this situation the ASIC potential will shift the bands in an almost identical way, without producing any quantitative changes, such as the opening up of the KS gap. Note that this is a pitfall of the ASIC approximation, which distinguishes occupied from empty states only through their projected atomic orbital occupation, but not of the SIC in general. Typical cases are those of Si and Ge. The KS gap in Si goes from 0.48 eV in LSDA to only 0.57 eV for ASIC 1/2 , while Ge is a metal in both cases. In addition the LSDA calculated valence bandwidths of 12.2 eV for Si and 12.8 eV for Ge, in good agreement with experiments, are erroneously broadened to 14.3 eV and 14.8 eV respectively.

VI. CONCLUSIONS
In conclusion, we have implemented the ASIC scheme proposed by Filipetti and Spaldin within the pseudopotential and localized orbital framework of the Siesta code. We have then investigated a broad range of semiconductors and molecules, with the aim of providing a reasonable estimate for the scaling parameter α. We found that α = 1, which accounts for the full atomic SI, describes surprisingly well ionic semiconductors and molecules. In particular for molecules, both the IP and the EA can be obtained with good accuracy from the HOMO KS eigenvalues respectively for the neutral and singly charged molecule. This makes the ASIC scheme particularly suited for application such as quantum transport, where the position of the HOMO level determines most of the I-V curve.
In contrast III-V and II-VI semiconductors are better described by α = 1/2, which corrects the atomic SI for screening. This makes ASIC 1/2 an interesting effective band theory for semiconductors. The relation of the present scheme with the fully self-consistent SIC methods has been emphasized, and so has been that with LDA+U .