Atomic-orbital-based approximate self-interaction correction scheme for molecules and solids

We present an atomic-orbital-based approximate scheme for self-interaction correction (SIC) to the local-density approximation (LDA) 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 nonvariational 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 atomiclike corrective scheme $\mathrm{LDA}+U$. We find that when no re-scaling is applied, i.e., when we consider the entire atomic correction, the Kohn-Sham highest occupied molecular orbital (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 band gaps 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 localdensity approximation ͑LDA͒ 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 nonvariational 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 atomiclike 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 highest occupied molecular orbital ͑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 band gaps 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

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, biomolecules, 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 ͑Ref. 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 regard 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 noninteracting system, v͑r͒ is the external potential, U is the Hartree electrostatic energy, and E xc is 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 parametrizations 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 it is 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 selfinteraction of occupied orbitals is canceled exactly by the nonlocal exchange. However, when using semilocal 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 selfinteraction 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 charge-transfer insulators. Moreover, the KS highest occupied molecular orbital ͑HOMO͒, the only KS eigenvalue that can be rigorously associated to a singleparticle 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 Semilocal 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 semilocal 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 ͑SICs͒ 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 semilocal 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 densityfunctional 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 orbital dependent. 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 single-particle 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 singleparticle equations, where the effective single-particle potential v eff,n ͑r͒ is de- These are solved in the standard KS way for atoms, with good results in terms of quasiparticle 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 nonorthogonality 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 size-inconsistent, 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 atomiclike 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 SICs.
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 xc SIC and canonical ͑Kohn-Sham͒ delocalized 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 ⑀ nm ,SIC , 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 Eq. ͑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 n SIC = u͓͑ n ͔ ; r͒ + v xc ,LDA ͓͑ n ↑ ,0͔ ; r͒. 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 Eq. ͑10͒ is that of using the so-called "unified Hamiltonian" method. 30 This is defined as ͑we drop the spin index ͒

͑13͒
where P n = ͉ n ͗͘ n ͉ is the projector over the occupied orbital n , and Q is the projector over the unoccupied ones Q =1 − ͚ n occup P n . The crucial point is that the diagonal elements of the matrix ⑀ nm ,SIC and their corresponding orbitals n are, respectively, eigenvalues and eigenvectors of H u , from which the whole ⑀ nm ,SIC can be constructed. Finally, and perhaps most important, at the minimum of the SIC functional, the canonical orbitals diagonalize the matrix ⑀ nm ,SIC , whose eigenvalues can now be interpreted as an analog 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 Eq. ͑10͒, where ⑀ nm ,SIC = ⑀ nm ,SIC ͑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-metal 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-metal monoxides, 42,44,45 rare-earth materials, 46 diluted magnetic semiconductors, 47 Fe 3 O 4 , 48 heavy element 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 atomiclike contribution, instead of being neglected. This atomic SIC ͑ASIC͒ scheme is certainly a drastic approximation, since it implicitly assumes that the transformation M minimizing the SIC functional leads to atomiclike 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 semicore 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 self-consistent SIC approach. Our practical implementation is also described.

A. ASIC potential
The starting point of our analysis is the SIC Schrödingerlike equation ͑10͒ for the canonical orbitals. Let us assume, as from Ref. 51, that the rotation M transforming localized orbitals ͑to be determined͒ into canonical orbitals ͓see Eq. ͑11͔͒ does not mix core and valence states. We also assume that core electrons are well localized into atomiclike wavefunctions and that they can be effectively described by a norm-conserving pseudopotential.
Let us now assume that M is known and so are the localized orbitals m . In this case the canonical orbitals diagonalize ⑀ nm ,SIC and Eq. ͑10͒ simply reduces to with ⌬v n SIC defined in Eq. ͑12͒. The Hamiltonian H 0 + ⌬v n SIC can be then rewritten in a convenient form as where v m SIC is the self-interaction potential for the localized orbital m , and P m is the projector over the same state,

͑16͒
Two main approximations are then taken in the ASIC approach. 1,25 First the localized states m are assumed to be atomiclike orbitals ⌽ m ͑ASIC orbitals͒ and the SIC potential is approximated as P m ⌽ is the projector of Eq. ͑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 atomiclike 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 xc SIC are Wannier-like functions. 60 The second approximation taken in the ASIC method is that of replacing the nonlocal projector P m ⌽ with its expectation value. The idea is that the SIC potential for the canonical orbitals ⌬v n SIC is formally a weighted average of the SIC potential for the localized orbitals v m SIC . For the exact SIC method the weighting factor is the nonlocal 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 Vogel and co-workers [25][26][27] this was simply set to 1. Here we consider the orbital occupation p m of the given ASIC orbital ⌽ m , i.e., we replace P m ⌽ with its expectation value 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 Ref. 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 atomiclike 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 widegap insulators.

B. Implementation
The final form of the SIC potential to subtract from the LSDA ͑local spin-density approximation͒ one ͓Eq. ͑19͔͒ is that of a linear combination of nonlocal 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 plane-wave and Gaussian orbital implementations are available, 1,25-27 while here we present our scheme based on the pseudoatomic orbital ͑PAO͒ ͑Ref. 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 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 norm-conserving 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 ṽ m SIC ͑r͒, which have a −2 / r asymptotic behavior, in a suitable nonlocal form. This is obtained with the standard Kleinman-Bylander 64 scheme and the final ASIC potential ͓Eq. ͑19͔͒ is written as where the ASIC projectors are given by ͑22͒ and the normalization factors are

͑23͒
The orbital functions ⌽ m Ј are atomiclike 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 pseudoatomic 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 ␦ SICm Ͻ0.1 mRy. Typically, when reasonable cutoff radii ͑6 -9 bohr͒ are used, we find that the atomic SIC-LSDA eigenvalues are reproduced to within 1 -5 mRy for the most extended shells and to within 0.1 mRy for more confined shells. Thus ␦ SICm are rather well converged already for cutoff radii defined by a PAO energies shifts 28 of around 20 mRy, 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 split-norm 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 subspace 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 xc LSDA 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 Eq. ͑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 Eq. ͑28͒ needs to be calculated. This is easily done by calculating both U and E xc LSDA 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 atomiclike 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 wherein we separate out the index for the atomic position I from the magnetic quantum number m, and introduce the off-diagonal populations p mn 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 mn I matrix. In this case the energy becomes simply . This reflects the fact that the SIC is defined only for occupied KS orbitals. An important consequence is that the opening of band gaps 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 role of the parameter ␣ in determining the band structure 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 quasiparticle 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 g KS , defined as the difference between the orbital energy of the HOMO and lowest occupied molecular orbital ͑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 This is valid even for the exact XC potential, and therefore in principle one has to give up KS band structures 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 g KS 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 g KS 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 Fig. 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 g KS increases as ␣ increases, as a result of the stronger SIC. The E g KS ͑␣͒ 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 ͑Ref. 41͒ demonstrate that the SIC band structure is rather insensitive of the localized orbitals once these have an atomiclike form.
For the other compounds the localized orbitals 's are not necessarily atomiclike 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 g exp to the calculated E g KS 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 g exp = E g KS . 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 materials 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 corresponds 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 bandwidth for the semiconductors investigated as calculated from LSDA ⌬E v LSDA and ASIC for both ␣ =1 ͑⌬E v ASIC 1 ͒ and ␣ = ␣ * ͑⌬E v ASIC ␣ * ͒. We also report the experimental values ⌬E v exp 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 bandwidth, it does just the opposite for KCl, SrO, and CaO.

B. Wide-gap semiconductors: ZnO and GaN
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 photoemission data disagree quite remarkably from LSDA calculations. In Fig. 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 position of the semicore Zn 3d states. The LDA binding energy for the Zn 3d states is ϳ5.5 eV while 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 bandwidth 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. First, 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. Second, 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. These results corroborate those of Refs. 1, 25, and 26 where in ZnO and GaN have previously been discussed in a pseudopotential based SIC context.

C. Transition-metal oxide: MnO
Transition-metal oxides like MnO and NiO are characterized by partially filled 3d orbitals and an associated local  Note that these are for the rhombohedral unit cell with four 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 antiferromagnetically coupled to each other. Also in this case, ASIC is a considerable improvement over LSDA. The size of the fundamental gap now resembles the experimental one and the VBT recovers some p character. We point the reader once again to Ref. 1 wherein transition-metal oxides have been discussed in much more detail.

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 ␣. How-  4. ͑Color online͒ Calculated band structure of antiferromagnetic MnO obtained from LSDA ͑left͒ and ASIC ͑right͒. In our calculation we obtain an LSDA band gap of ϳ0.65 eV whereas the ASIC band gap is much improved at ϳ3.5 eV. The VBT is aligned at 0 eV. ever, here we limit ourselves only to ␣ =1 ͑ASIC 1 ͒ and ␣ =1/2 ͑ASIC 1/2 ͒.
In Table III and Fig. 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, is 3.56 eV for LSDA, 1.69 eV for ASIC 1/2 , and only 0.58 eV for ASIC 1 ͑N runs over the molecules of 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 Fig. 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

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 Eq. ͑42͒ indicates, the EA is in principle accessible 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 noninteger 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 ⑀ N+1 HOMO ͒ of several singly negatively charged molecules with the experimental electron affinities. We also report the LUMO energies for the corresponding neutral species  ͑denoted as ⑀ N LUMO ͒. LSDA relaxed geometries for the neutral molecule are used for both the neutral and charged cases. We find that various ⑀ N+1 HOMO 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, 2.31, and 0.83 eV for LSDA, ASIC 1/2 , and ASIC 1 , respectively. Notice that ⑀ N+1 HOMO from LSDA is positive in most cases as the states are unbound.
In Fig. 6 we present our data together with ⑀ N+1 HOMO 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 SICs is less clear. As a general consideration, also for the case of vertical ionization energies selfconsistent SICs seem 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 ASICs. As already mentioned, even for the exact XC functional, the KS gap E g KS = ⑀ LUMO − ⑀ HOMO does not account for the actual quasiparticle gap E g = I N − A N . This in turn is the sum of E g KS 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  In Fig. 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 Eq. ͑41͒. The ASIC HOMO eigenvalue roughly agrees with the PZ-SIC eigenvalue at integer occupation but behaves linearly through noninteger 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 ASICs.
In view of the foregoing discussion, the actual size of the HOMO-LUMO gap in ASICs becomes significant with a di-rect 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 ⑀ ASIC LUMO ͑N͒ usually differs from ⑀ LSDA LUMO ͑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 atomicorbital 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 ASICs 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 unoccupied 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, it is 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 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 applications 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. FIG. 7. 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 state.