{'fraction_non_alphanumeric': 0.03186414806334669, 'fraction_numerical': 0.010268685712303343, 'mean_word_length': 4.664930726147195, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A self-interaction corrected pseudopotential scheme for magnetic and strongly-correlated systems

Local-spin-density functional calculations may be affected by severe errors when applied to the study of magnetic and strongly-correlated materials. Some of these faults can be traced back to the presence of the spurious self-interaction in the density functional. Since the application of a fully self-consistent self-interaction correction is highly demanding even for moderately large systems, we pursue a strategy of approximating the self-interaction corrected potential with a non-local, pseudopotential-like projector, first generated within the isolated atom and then updated during the self-consistent cycle in the crystal. This scheme, whose implementation is totally uncomplicated and particularly suited for the pseudopotental formalism, dramatically improves the LSDA results for a variety of compounds with a minimal increase of computing cost.


I. INTRODUCTION
Among the innumerable successes of the local-density (LDA) and local-spin-density (LSDA) approximations to density functional theory (DFT) there are also well known, systematic failures 1,2 that compromise the accuracy of predictions for a range of properties, especially in magnetic and strongly-correlated 3 compounds.
Typical examples of LSDA failures are the series of transition-metal monoxides, 4 improperly described in LSDA as either small-gap Mott-Hubbard antiferromagnetic insulators (MnO, NiO) 5,6 or even ferromagnetic and non-magnetic metals (FeO, CoO, CuO), 5,7 whereas according to experiments these materials are charge-transfer antiferromagnetic wide-gap insulators. A similar situation occurs for the high-T c compounds La 2 CuO 4 and YBa 2 Cu 3 O 6 which are described in LSDA as nonmagnetic metals instead of as antiferromagnetic insulators, 8 and for the perovskite manganites (e.g. La x Ca 1−x MnO 3 ), 9,10 for which the LSDA fails to predict the correct magnetic and orbital orderings. In general, LSDA favors metallic and ferromagnetic ground states over the observed antiferromagnetic insulating ground states. This is particularly harmful in the case of hexagonal YMnO 3 , which is antiferromagnetic and ferroelectric, but is described as a metal within LSDA, 11 thus preventing the possibility of calculating any ferroelectric properties at all.
These failures can be, at least in part, attributed to the presence of the self-interaction (SI) in the LSDA energy functional, i.e. the interaction of an electron charge with the Coulomb and exchange-correlation potential generated by the same electron. The SI vanishes in the thermodynamic limit for delocalized states, but is present in systems characterized by spatially localized electron charges like 2p, 3d and 4f electrons. As a consequence of the SI, the binding energies, the on-site Coulomb energies (i.e. the Hubbard U parameter) and the exchange-splitting of the d and f states are underestimated, whereas the hybridizations of cation d and anion p states and the corresponding band widths are overemphasized.
Since, for materials with partially filled d states, the on-site Coulomb energy should be much larger than the charge-transfer energy between p and d electrons, it is clear that the suppression of U and the overestimation of the p-d hybridization may change dramatically the character of the band structure, with a tendency to favor metallic and ferromagnetic ground states over insulating antiferromagnetic states. Furthermore, within the LSDA Kohn-Sham (KS) description based on a local, orbitalindependent potential, orbital and charge orderings cannot be properly accounted for. The presence of SI and the possible strategies for eliminating the SI in density functional theories are longstanding issues which go back to the Thomas-Fermi and the Slater (Xα) approaches. 12 A large amount of literature has been produced [13][14][15][16][17] within the context of the LSDA. Particularly fundamental are the works by Perdew and Zunger 15,17 which proposed a form of selfinteraction correction (SIC) within LSDA (SIC-LSDA) and successfully applied it to the calculation of atomic properties.
More recently, Svane and co-workers presented a long series of works [18][19][20][21][22][23][24][25][26][27][28] which attempted the highly challenging application of the fully self-consistent SIC-LSDA to extended systems with encouraging results. However the full SIC-LSDA requires a large computing effort when applied to extended systems even for materials with small unit cells, and becomes prohibitive for larger systems.
A useful alternative to this approach has been suggested by Vogel and co-workers. [29][30][31] They approximated the SI part of the KS potential with a non-local, atomiclike contribution included within the pseudopotential construction. This scheme, applied to non-magnetic II-VI semiconductors and III-V nitrides, is capable of remarkable improvements upon the LSDA results while keeping the computational cost comparable to that of an ordinary pseudopotential calculation.
Inspired by these results, in this paper we develop an approach which is in the same spirit as that of Ref. 29, but can be applied to more general cases, and in particular to magnetic and highly-correlated systems where there is a coexistence of strongly localized and hybridized electron charges. The main innovation of our formalism is the introduction of orbital occupation numbers in the electron KS Hamiltonian: The pseudopotential-like SIC (pseudo-SIC) is rescaled by the occupation numbers calculated self-consistently within the crystal environment, so that the SI coming from localized, hybridized, or completely itinerant states may be discriminated accordingly. Also, these occupation numbers discriminate between valence and conduction bands so that only the former are corrected, in accord with the idea that conduction states are itinerant, thus SI-free. Furthermore, our scheme is generally applicable to insulators as well as metals, without the necessity of knowing the character of the material a priori. This is essential for compounds which are erroneously described as metals within LSDA, such as the hexagonal YMnO 3 , treated in section IV E.
A second major difference between our approach and Ref. 29 lies in the way in which relaxation corrections to the SIC potential are accounted for. In Ref. 29 they are incorporated as an additional, orbital-dependent number which is imported from atomic energy difference calculations (i.e. the so-called ∆SCF method). In contrast, in our case the relaxation contribution is included directly as an analytical correction in the SIC potential projector, and no further atomic calculations are needed.
In this paper we test the pseudo-SIC on three classes of materials, including the wide-band gap insulators (e.g. ZnO and GaN), the transition-metal oxides (MnO and NiO) and the manganese oxides (YMnO 3 ). Our main focus is on the electronic properties (i.e. the band energy structure) of these materials, since we want to verify the capability of our single-electron Hamiltonian to reproduce the main features of photoemission spectra. For some of these compounds we also calculated the theoretical structure by total energy minimization. In general we obtain very encouraging results and clear systematic improvements over the LSDA description, without considerable increase of computing effort.
Note that we have implemented the pseudo-SIC within the ultrasoft pseudopotential method 34 (USPP). This is instrumental in keeping a moderate computational cost even for large unit cells containing atoms with 2p and 3d electrons.
As is customary when a new scheme is introduced, we compare our results with those of other common beyond-LDA approaches. According to the set of results presented here, the pseudo-SIC ranks among the most accurate. In particular, it seems to perform equally well (or even better) than the very popular LDA+U. 32 Although a detailed comparison between pseudo-SIC and LDA+U is not the aim of this paper, we can point out some potential advantages for the pseudo-SIC. First, it does not require parameters imported from an external theory. Second, it can be applied to both magnetic and non-magnetic compounds, whereas the LDA+U is constructed to correct spin-polarized and/or orbital-ordered band structures. Finally, within LDA+U or even within the exact SIC-LSDA a choice of which orbitals are localized in space and thus, which orbitals are to be corrected, has to be made before the calculation, and the final result depends on this assumption. Instead, in pseudo-SIC the correction is applied to all the orbitals indiscriminately, and no choice of orbital localization is required. This is an advantage in case where a discrimination of the electron charge localization cannot be established a priori. This happens in non-bulk systems such as surfaces and interfaces, when localized surface states or resonance states are present, or even in bulk materials whenever the charge localization is not, or not only, due to the chemical nature of the involved electrons. Examples of this kind are the superconductor cuprates, in which the d orbital charges have different localizations corresponding to different physical solutions.
The remainder of this paper is organized as follows: in Section II we review the main features of earlier work on SIC implementations in LSDA. In Section III we describe the pseudo-SIC formulation applied in combination with norm-conserving pseudopotentials (NCPP). The extension of the pseudo-SIC to USPP is given in the Appendix. In Section IV we report our results obtained within the pseudo-SIC approach, and finally in Section V we present the conclusions.

II. OVERVIEW OF PREVIOUS WORK
The SI effects in atomic calculations are extensively explained in the seminal paper by Perdew and Zunger. 17 In the following we briefly summarize the most fundamental features. Due to the SI, the tail of the KS potential does not recover its physical long-range limit, -1/r. Thus negative ions which are stable experimentally cannot be described at all in LSDA since the outermost eigenstates are unbound. The atomic total energies are severely overestimated with respect to the experimental values, as a consequence of the reduction in binding energy caused by the spurious self-screening in the KS potential.
Consistently, the LSDA KS eigenvalues strongly overestimate the experimental electron removal energies. Although the KS eigenvalues represent in general Lagrange multipliers and the comparison with the removal energies may be questioned, calculations 17 clearly show that in atoms the SI represents, by far, the largest source of mismatch. Furthermore, it can be proved that the highest occupied KS eigenvalue of the exact density functional theory equals the atomic ionization potential. 35 Thus, at least for this quantity, the discrepancy must be attributed to the LSDA inaccuracy.
The SIC-LSDA proposed by Perdew and Zunger is based on the straightforward subtraction of the SI contribution from the LSDA KS potential: where n and m are the total charge and magnetization densities, and n σ i is the spin-charge density of the i th orbital. It is understood that in V σ HXC [n σ i ] the magnetization density is set equal to 1, i.e. the single electron charge is fully spin-polarized. Despite its conceptual simplicity, this modification is able to systematically improve total energies, ionization potentials and electron affinities, giving a much better agreement between KS eigenvalues and removal or addition energies.
Notwithstanding this initial success, the SIC-LSDA has not enjoyed a large popularity in the first-principles community, mainly because the correction severely compromises the feasibility of the LSDA, introducing an orbital-dependent, spatially localized (i.e. non-periodic) contribution in the KS potential. Thus different wavefunctions experience different Hamiltonians, and are no longer orthogonal. This is not a serious drawback in atomic calculations, since the non-orthogonality is a small effect and the atomic orbitals can be easily forced to be orthogonal during the self-consistent cycle. 17 In contrast, the direct application of Eq.1 to extended systems (where the solutions of the KS equations are Bloch states) is very challenging. First of all the SI of a Bloch state is an ill-defined quantity that depends on the normalization of the wavefunction, and vanishes as Ω −1/3 in the thermodynamic limit, where Ω is the system volume. Thus, if we assume that the electron charges are properly described as Bloch states the SIC becomes immaterial and discardable.
Clearly this is not necessarily true in general. In cases like simple metals or semiconductors with mostly covalent interactions (e.g. bulk Si) we can safely assume that the SI is discardable. However, in many systems the electron charges retain atomic-like features such as small band dispersion and a spatial distribution strongly localized around their ion-cores. A clear example of the variable influence of the SI is given by the energy gap of semiconductors: for bulk Si the discrepancy between the theoretical and experimental energy gaps (0.7 eV and 1.17 eV, respectively) must be attributed to non-locality and many-body effects. Instead, in LiF the SIC accounts for almost 95% of the LDA gap error. 36 Thus, in order to have a non trivial SIC, a description in terms of localized orbitals must be adopted, which causes violation of the translational invariance and consequently the inapplicability of the Bloch theorem. Furthermore, the SIC-LSDA eigenstates are not invariant under unitary rotations within the subspace of the occupied orbitals, and the SIC-LSDA solutions strictly depend on the assumption of choosing each orbital as extended (thus self-interaction free) or localized (therefore subject to SIC). We will show an example of this in Section IV D.
To our knowledge the implementation of a fully selfconsistent SIC-LSDA approach for extended systems was pioneered by Svane and co-workers. 18 They carried out a series of applications to a remarkable range of materials including the family of transition-metal monoxides, 18,24 the high-T c superconductor parent compounds La 2 CuO 4 19,20 and YBa 2 Cu 3 O 6 , 27 the rare-earth materials γ-Ce and α-Ce 22,23,25 , and Yb 26 and Pu 28 monopnictides and monochalcogenides, obtaining systematic improvements over the LSDA results. Other SIC-LSDA applications to transition-metal oxides, 21,33 implemented within different computational methodologies, are also present in the literature. The major drawbacks of these SIC-LSDA implementations are the rather complicated formulation with respect to the LSDA and the increase of computing cost which makes the SIC-LSDA almost impractical for large systems. One reason for this increase is the use of big supercells needed to describe the localized orbitals (e.g. ∼ 500 atoms for bulk MnO 19,33 ). Indeed, in the works previously cited, the SIC-LSDA is implemented using a basis of linear muffin tin orbitals within the atomic sphere approximation (LMTO-ASA), whereas, to our knowledge, there are no examples of implementations using the more expensive plane-wave basis and pseudopotentials.
An important step towards a practical (albeit approximate) expression for the SIC-LDA was accomplished by Vogel and co-workers 29 by incorporating the atomic SIC within the non-local pseudopotential projectors (SIC-PP) generated from the free atom. The idea underlying this approach is that the SI potential of a localized electron in the crystal can be well approximated by the SI which the same electron experiences in the free atom. The SIC-PP turned out to be quite efficient for describing the properties of some highly ionic compounds with atomic-like, poorly hybridized bands, such as II-VI semiconductors, 29 III-V nitrides, 30 and silver halides. 31 For these materials, the energy band structures calculated with the SIC-PP show a much better agreement with photoemission spectroscopy measurements than the LDA band structure calculations.
Our pseudo-SIC can be considered a generalization of the SIC-PP approach. It is still based on the idea of replacing the SIC potential with a non-local projector, but now the projector depends on the orbital occupation numbers calculated self-consistently within the crystal. The details are described in the next section.

A. Kohn-Sham equations within pseudo-SIC
Although the calculations reported in this paper have been performed within the USPP method, 34 here we describe the pseudo-SIC formalism adapted to NCPP. This allows a simpler formulation and an easier understanding of the logic behind this approach. (The generalization to USPP is given in the Appendix.) Within pseudo-SIC, the SIC potential is cast in terms of a non-local projector, which resembles the non-local part of the pseudopotential: Here n and m are the periodic charge and magnetization densities of the crystal and the SI is written in terms of atomic quantities only: i = [(l i , m i ), R i ] is a cumulative index for angular momentum quantum numbers and atomic coordinates, Y i are projector functions (e.g. spherical harmonics), and n σ i the charge densities of the (pseudo) atomic orbitals φ i : where p σ i are orbital occupation numbers. Through Eq.2, the Bloch wavefunctions are projected onto the basis of the atomic orbital charges n σ i . For each projection, the Bloch state is corrected by an amount corresponding to the atomic SIC potential V σ HXC [n σ i ]. Thus the SIC is applied without really introducing a dependence of the KS Hamiltonian on the individual Bloch wavefunctions, and the difficulties of the SIC-LSDA approach are overcome.
The presence of p σ i in Eq. 3 is a major novelty of our approach. In Ref. 29 the occupation numbers are implicitly set to 1, i.e. the SIC potentials are generated from fully occupied atomic states. Clearly, this cannot be a correct choice in general, since the occupation numbers may become fractions whenever hybridization, degeneracy, or spin-polarization effects arise. Furthermore, moving from the free atom to the crystal, the p σ i can change a great deal. For instance atomic orbitals which are fully occupied in the atomic ground state may become conduction states in the crystal (e.g. Zn 4s in ZnO).
Thus the p σ i must be allowed to be fractional and must be recalculated self-consistently within the crystal. Within a plane-wave basis set it is straightforward to calculate p σ i as atomic orbital projections onto the manifold of the occupied Bloch states: These quantities are analogous to the local orbital occupations calculated within the FLAPW implementation of the LDA+U. 52,53 In the limit of an isolated atom the p σ i recover the atomic values, whereas in the case of hybridized bonds, they rescale the atomic SIC by the amount of charge that actually occupies the atomic orbital. This fact has a fundamental consequence: since in most cases empty bands have dominant characters from orbitals whose occupation numbers are close to zero, the SIC potential will not affect these bands. This is consistent with the general assumption that conduction states are itinerant, thus SI-free. (In fact, even if the conduction states are not itinerant, it is theoretically justified to apply the SIC only to the occupied states, since these are the states which actually see their own charge. 48 ) Notice that partially occupied bands can be treated on the same footing as filled bands, thus the scheme can be applied to both insulators and metals.
Following the suggestion of Ref. 29, for the purpose of numerical efficiency we recast the pseudo-SIC as a fully non-local, Kleinman-Bylander projector: where and The pseudo-SIC KS equations are: whereV P P is the pseudopotential projector, and ǫ σ nk are the KS eigenvalues. The recalculation of V σ HXC [n σ i ] at each iteration of the self-consistency for each atom and angular component would result in a major increase of computing cost. A large saving of time can be achieved by assuming a linear dependence of the SI potential on the occupation numbers: so that V σ HXC [n σ i ; 1] (i.e. the SI potential for the fullyoccupied orbital) is set in the inizialization, and only the p σ i needs to be updated during the self-consistency. Eq.9 is exact for the Hartree term, which is the dominant contribution for large occupation numbers, whereas it introduces a non-linearity error of O(p 1/3 i − p i ) in the exchange-correlation part. 17 The time saving provided by the linear scaling argument is instrumental for calculating structural relaxations and electronic properties of large-sized systems. Furthermore, the assumption of linear scaling allows us to introduce relaxation effects in a very simple way into the pseudo-SIC scheme. Indeed if an electron state is localized its energy will change with the orbital occupation. Thus, in order to compare the calculated eigenvalue with the photoemission spectroscopy data (i.e. with the electron removal energy) the effects of the electron relaxation must be subtracted out of the one-electron potential. In DFT the energy required to remove a fraction p of an electron from a one-electron localized state is: [37][38][39] where ǫ is the corresponding KS eigenvalue. The leading dependence of the LSDA KS potential (and eigenvalues) on the orbital occupations is indeed contained in its SI part. Thus, if δǫ is the SI part of the LDA eigenvalue, within linear scaling we have ∆E(p) = 1/2 p 2 δǫ(1). It follows that the SIC relaxation energy, 1/2 p 2 δǫ(1), can be subtracted out by rescaling the SIC potential as follows: Through Eq.11 we directly incorporate in the pseudo-SIC KS Equations the electron removal energy due to the SI contribution. We point out that the eigenvalue relaxation is very important in order to match the photoemission spectroscopy results. Discarding this contribution, the SIC eigenvalues would strongly overestimate the electron removal energies.
Notice that, if the KS eigenvalues did not depend on the occupation numbers, according to Eq.10 they would be equal to the electron removal energies. This is in fact the case in Hartree-Fock calculations and in any other theory which obeys Koopman's theorem. This property does not hold in LDA or any LDA-related scheme (such as GGA, SIC, etc.) where the potential explicitly depends on the orbital occupation numbers.
As is customary in atomic calculations, we assume the radial approximation for the atomic orbital charges, so that the SIC projectors can be written: (12) where ν labels the atom type, and n σ l,ν (r) is the radial pseudo-charge density of orbital (l,m l ) (in radial approximation this does not depend on m l ). Finally, the normalization coefficients are: In summary, except for the electron occupation numbers, all the other atomic-like ingredients can be imported from the code used for the pseudopotential generation, and the Kleinman-Bylander projectors for the SIC are set during the initialization process. Since the calculation of the p σ l,m l ,ν through Eq.4 is not very demanding, the global computational cost of the pseudo-SIC for each self-consistent iteration is roughly equal to that of the basic LSDA. However within pseudo-SIC the number of self-consistent iterations required to reach the self-consistency can be larger than in LSDA due to oscillations in the values of the occupation numbers.

B. Total Energy within pseudo-SIC
In the previous section we estabilished the form of the pseudo-SIC KS Equations. Here we formulate a suitable expression for the total energy functional. We point out that, within our scheme, a physically meaningful energy functional which is also related to Eqs. 8 by a variational principle is not available. (It is, by construction, within LSDA and SIC-LSDA.) In SIC-LSDA the energy functional is: 17 where E[n, m] is the LSDA energy functional and E HXC [n σ i ] the Hartree exchange-correlation energy of the i th fully spin-polarized electron charge: where V H is the Hartree potential and E XC is the local exchange-correlation energy density. For the pseudo-SIC total energy we adopt the same expression as Eqs.14 and 15, with the orbital charges n σ i given by Eqs.3 and 4. Notice that, Eq.16 represents the Janak thoerem 37,39 applied to the SIC contribution since C σ i (see Eq.7) is the SI part of the atomic eigenvalue. In terms of the pseudo-SIC eigenvalues (Eq. 8) the total energy can be rewritten: where E ion is the usual Ewald term.
i,σ E HXC [n σ i ] produces a gentle modification to the LSDA energy functional whereas the fifth term in Eq.17 is a strongly varying contribution which compensates the same contribution present in the pseudo-SIC eigenvalues. Without this compensation, Eq.17 would give very inaccurate total energies which would be unphysically far from the LSDA values. A numerical example of this behavior will be shown in Section IV C.

A. Technical details
In this work we compare results from our plane wave USPP 34 pseudo-SIC implementation with results from conventional LSDA USPP method. The local exchange-correlation energy functional is modeled using the the Perdew-Zunger interpolation formula. 17 The use of USPP 34 allows us to obtain well-converged results with moderate cut-off energies (35 Ry for MnO and YMnO 3 , 40 Ry for ZnO and 45 Ry for GaN). In order to have highly transferable USPP, two projectors per angular channel are included for all the atoms. The 3d 10 electrons in ZnO and GaN are treated as valence states, while the semicore Y s and p electrons are placed in the valence for the LSDA calculations and in the core for the pseudo-SIC calculations (the reason will be explained in the following section). For total energy calculations we use up to 8×8×8 grids of special k-points. 40

B. Atomic ingredients for the pseudo-SIC
The atomic quantities necessary to build the pseudo-SIC hamiltonian are the pseudo atomic orbitals φ l,m l ,ν , required in the Kleinman-Bylander projector and for the calculation of p σ i , and the SIC potentials V σ HXC [n σ l,ν (r); 1] for the respective pseudo orbital charges calculated at full electron occupation (we use pseudopotential and not allelectron atomic functions for obvious reasons of smoothness.) To illustrate the impact of the atomic SIC, in Figure  1 we compare the Zn s, p, and d (unpolarized) pseudopotentials (V l ) to the same quantities minus the corresponding corrections ∆V SIC l = V HXC [n l,ν ; 1]. As expected, the SIC makes the electron potentials more attractive and recovers the physical long-range limit -2/r (in Ry).
Note that, before transfering V HXC [n l,ν ; 1] to the crystal, the long-range tails must be cut off, since the SIC should act as a local correction: each SIC potential must be directly applied only to the electron states localized on the same atom, otherwise the overlapping Coulomb tails would give rise to an unphysical SIC overestimation. In Figure 2 we report V HXC [n l,ν ; 1] (times r) for both pseudo and all-electron orbital charges (the latter show the cusps corresponding to the wavefunction nodes), and the pseudo orbital charges ρ l . Since only the products V HXC [n l,ν ] ρ l contribute in the Kleinmann-Bylander construction, the SIC potentials can be cut off as soon as the orbital charges vanish, without losing accuracy. The correction lowers the one-electron potential and recovers the physical long-range limit -2/r which is missing in LDA.
It is important to notice that V σ HXC [n σ l,ν ; 1] is not very sensitive to the atomic valence configuration from which it is actually generated. For example the integrated value of V σ HXC [n 3d ; 1] · φ 2 3d (i.e. the atomic SI) changes by only ∼ 1% if calculated in the 3d 10 4s 2 , the 3d 10 4s 1 , or the 3d 10 4s 0 configuration. In other words, despite being valence dependent, V σ HXC [n l,ν ; 1] is fairly transferable. Notice also that the pseudo-SIC does not interfere with the pseudopotential construction, and the same pseudopotentials used for the LSDA calculations can be used for the pseudo-SIC as well, due to the presence of the occupation numbers. Indeed, in the limit p i =0 the SIC vanishes, and the hamiltonian must continuously recover its LSDA value.
However for pseudo-SIC calculations we prefer to build pseudopotentials which take into account the SIC for the core states. Although this correction does not significantly modify the valence properties, 41 there are reasons related to pseudopotential transferability which suggest that the SIC should be included in the core states. First, the SIC could push strongly localized valence states (e.g. the semicore states) down in energy and unphysically close to the core states if the latter are not shifted down consistently. Second, the application of SIC to the semicore states may better justify their inclusion in the core even if in LSDA they need to be treated as valence states. This is the case for the Y 4s and 4p electrons, and the Ga 3d electrons. A third reason is even more fundamental: increasing the binding energy and the space localization of the core states makes them less sensitive to charge relaxation or chemical activity. This reinforces the hypothesis that core states are insensitive to the chemical environment, which is at the basis of the pseudopotential approximation.
The transferability of the atomic SIC to the extended system can be tested with the same procedure used for the pseudopotentials: the atomic pseudo-SIC hamiltonian should recover the all-electron, fully self-consistent SIC-LSDA eigenvalues. We generally find that pseudo-SIC and SIC-LSDA atomic eigenvalues agree within ∼ 1% Finally we stress a fundamental point: our procedure is not in any way restricted to the use of pseudopotentials or any other particular choice of basis functions. For example it can be equally well implemented within the context of all-electron methods. The advantage of using the pseudopotential approach relies on the fact that the SIC projectors can be easily assembed from quantities (such as the projections of Bloch states onto atomic orbitals) which are usually calculated in any pseudopotential-based computing code.

C. Wide-gap semiconductors
We begin our series of pseudo-SIC applications with two of the most prototypical wide-gap semiconductors, wurtzite ZnO and GaN. These materials are ideally suited as test cases since the LDA predictions for the band energies show some easily recognizable discrepancies with the abundantly available photoemission spectroscopy data. Furthermore, these compounds are of technological interest due to their applications in optoelectronic and piezoelectric devices. Nowadays the electric polarization and related properties can be efficiently calculated by first-principles techniques. 42 Thus a theory able to repair the LDA description of the electronic structure will also be important for the accurate evaluation of dielectric and piezoelectric response.
Since ZnO and GaN show similar structural and electronic behavior, we describe them together. In Table I we report our calculations within LDA and pseudo-SIC for the equilibrium structural parameters. These have been evaluated by energy minimization within the full space of parameters a, c and u (u is the anion-cation distance in units of c). In general, the SIC volume is larger than the LDA volume, since the stronger electron localization due to the SIC enhances the electron screening and reduces the electron-ion interaction and the charge hybridization. This is generally a favorable correction since it is well known that within LDA (or LSDA) the lattice parameters are underestimated compared with the experimental values. In the case of ZnO and GaN, however, our LDA values are already in excellent agreement (within less than ∼ 1%) with the experiments and no correction would be required. Happily, the pseudo-SIC calculations do not overcorrect the LDA values. (the reason will be explained later on the basis of the band structure results). This is an important aspect since typical beyond-LDA methodologies, also not derived by variational principles, often improve the LDA description of the electron excitation spectra but are less accurate than LDA in the determination of the equilibrium structure. In Figure 3 our calculated LDA energy band structure of wurtzite ZnO is shown. For such strongly ionic compounds each group of bands can be clearly labeled according to a single, dominant orbital character as indicated in the Figure. With respect to the photoemission spectroscopy results, the 3d bands calculated within LDA are too high in energy and overlap with the sp 3 valence band manifold. This produces a spurious p-d hybridization which shrinks the energy range of the sp 3 bands. Furthermore, there is the notorious problem of the fundamental energy gap which is underestimated by ∼ 40 % (see Table 3). All these undesirable features are largely overcome in our calculated pseudo-SIC band structure (Fig. 4). The largest SIC effect is a downshift of the d band energies, which are now placed ∼ 3 eV below the center of the sp 3 -band manifold, in excellent agreement with the experiments. Accordingly, the sp 3 band manifold is ∼ 1 eV broader than in LDA.
Furthermore, the bands with O p orbital character, which are the main contributors to the valence band top (VBT), are also corrected for an amount of SI corresponding to their average electron occupation (∼ 0.8 for the O p orbitals). In contrast the conduction band bottom (CBB), which is mainly Zn s in character, is almost unchanged, since the Zn s orbital occupation is ∼ 0.1. As a consequence the pseudo-SIC energy gap opens up to a value fairly close to the experimental band gap.
Our results show a substantive agreement with the calculated SIC ZnO band structure of Ref. 29 (they obtain E g = 3.5 eV). However the two methods act in quite different ways: in the approach of Ref. 29 all the bands (occupied and unoccupied) are corrected by the atomic SIC potential of the corresponding fully occupied orbital charges. Thus in Ref. 29 the increase of the energy gap with respect to the LDA can be roughly quantified as the difference between the SIC of the fully occupied Zn 4s and O 2p atomic eigenvalues. Instead, in our scheme only the occupied bands are corrected. We belive that this is a more conceptually sound approach for the reasons explained in Section III A.
In Figure 5 we show the LDA band structure of wurtzite GaN. As is usual in LDA calculations, the Ga d band energies fall within the same energy range as the N s bands, in disagreement with the experiments which place the d bands ∼ 3 eV below the N s bands. 49 This gives rise to a spurious s-d hybridization which causes an in-crease of lattice constant due to the closed-shell repulsion enhanced by the resonances between states on neighboring sites [44][45][46] (i.e. N s and Ga d). In our calculation this effect compensates the LDA tendency of underestimating the lattice parameters, and causes a better-thanusual agreement with the experimental data (Tab. I). However, the description of the band structure appears grossly inappropriate when compared to the X-ray photoelectron spectra. 49 The s-d hybridization splits the N s bands into two sections, one placed above and the other below the d bands. The s-d manifold is positioned only ∼ 6 eV below the bottom of the valence sp 3 manifold, and the energy gap is 40% smaller than the experimental value. The pseudo-SIC corrects in large part the LDA description. In Figure 6 we show the GaN band structure calculated within pseudo-SIC. Here the Ga d bands (E d ∼ -18 eV) are correctly located ∼ 3 eV below the center of the N 2s bands, thus the spurious s-d hybridization is avoided and the closed-shell repulsion reduced. This explains why the equilibrium lattice parameters calculated within pseudo-SIC are close to the LDA values although the pseudo-SIC is generally expected to enhance the electron-ion screening.
Also, the binding energies of Ga d and N s states are increased by ∼ 4 eV and 2 eV respectively, thus correcting in large part the faults of the LDA description. However within pseudo-SIC the energy gap is somewhat overcorrected (it is ∼ 0.9 eV larger than the experimental value). To conclude this Section, in Figure 7 we analyze the behavior of the pseudo-SIC total energy (Eq.17) for zincblende GaN as a function of the lattice constant in the region around its experimental value a 0 =8.505 Bohr. E 1 is the functional given by the first four terms of Eq.17. Clearly, its minimum value is far out of the physically meaningful region. E 2 is obtained by adding to E 1 the fifth term of Eq.17, thus E 2 equals the LSDA energy functional (except for the difference in the wavefunctions). This functional has a minimum value in excellent agreement with the experiments. Finally, E SIC is E 2 + i,σ E HXC [n σ i ]. This last term is almost independent of the lattice parameter, since it only changes through the p σ i , which are rather local quantities. As a consequence, E SIC and E 2 have similar behavior.
The LSDA misrepresentations of the electronic properties of the transition-metal oxides originate mainly from the description of the d electron states. At variance with the non-magnetic semiconductors considered in the previous section, in the transition-metal oxides the d states lie higher in energy and closer to the fundamental gap, often overlapping with the oxygen p states. Since in LSDA the d electron binding energies are severely underestimated, the d character is generally dominant at the VBT, so these compounds are described as Mott-Hubbard insulators. This is in striking contrast with the photoemission spectroscopy data which describe these materials as charge-transfer insulators (or eventually in the intermediate charge-transfer Mott-Hubbard regime) with a majority or even dominant O p character at the VBT. 7 Furthermore the LSDA energy gaps for these materials are even more severely underestimated than those of the III-V or II-VI semiconductors and in some cases the gap can be vanishing. This is because the major driving force leading to the formation of a gap separating bands of equal orbital character should be the on-site Coulomb energy U, but in LSDA the energy gap can only open due to Hund's rule and the crystal field splitting. These are both of order 1 eV, i.e. almost one order of magnitude smaller than U, and therefore can be easily overcome by the the band dispersion.
Finally, since the LSDA tends to broaden the space distribution of the electron charge and emphasize the interatomic charge hybridizations, the local magnetic moments are generally underestimated with respect to the experimental values. As applications of the pseudo-SIC to magnetic compounds we choose two of the most widely studied transition-metal monoxides, MnO and NiO. For these materials we can compare our results with both experimental data and the output of some common methodologies used for repairing the faults of the LSDA description. 5,7 These schemes include the SIC-LSDA, 18,24,33 the LDA+U, 32,24,51-53 and the model GW. [54][55][56][57] The ground state of MnO and NiO is characterized by A-type antiferromagnetic ordering, which consists of (111) ferromagnetic layers of Mn (or Ni) with alternating spin directions, while the oxygens are non-magnetic. The symmetry is rhombohedral, with 4 atoms per primitive unit cell. In our calculations we fix the lattice constant to its experimental value (a 0 = 8.37 Bohr and 7.93 Bohr for MnO and NiO, respectively).
In Figure 8 the LSDA band structure for MnO is shown. Notice that the orbital characters are referred not to the cubic, but to the hexagonal coordinates, i.e. x and y lie on the hexagonal (111) plane, and z is parallel to the [111] axis. Also, within the hexagonal (or rhombohedral) crystal field, the 3d orbitals are split into two doublets (d xy , d x 2 −y 2 ) and (d xz , d yz ), and one singlet (d z 2 ).
The nominal ionic configuration Mn 2+ O 2− suggests a Hund's rule-induced splitting between filled spin up and empty spin down 3d bands, resulting in an energy gap. From the analysis of the band character we see that the VBT is composed of 60% (d ↑ xy , d ↑ x 2 −y 2 ) and (d ↑ xz , d ↑ yz ) orbitals, with the remaining 40% coming evenly from the (p x , p y ) doublets of the two equivalent oxygens. In contrast, the conduction band bottom (CBB) consists of ∼ 80% d ↓ z 2 centered on the same Mn, which is equivalent to the d ↑ z 2 orbital localized on the other Mn. Thus, from the LSDA calculation, a picture of the MnO as a smallgap Mott-Hubbard insulator emerges, which is not what should be expected according to the experiments.
Furthermore, while the value of the Mn magnetic moments is satisfactory, the size of the LSDA energy gap severely underestimates the measured value (see Table  II). Also, notice that the occupied 3d bands are rather flat and well separated by a ∼ 1 eV gap from the underlying O 2p band manifold which is ∼ 6 eV wide.
In Figure 9 we see how the application of the pseudo-SIC modifies the results. First of all, there is no longer a gap between 3d and 2p bands. The SIC causes a downward shift of the Mn 3d bands with respect to the more dispersed O 2p bands, thus the 3d character is now spread across the whole (8 eV wide) valence band manifold, and is heavily mixed with the O 2p character.
This shift significantly increases the energy gap which is now well within the experimental uncertainty (see Table II). It also changes the VBT character, which is now composed of 40% 3d doublets and 60% 2p orbitals. 85% of the CBB still comes from the d ↓ z 2 orbital. Thus the agreement with experiments, which locate the MnO in the intermediate charge-transfer Mott-Hubbard regime, is restored. Furthermore, the enhanced on-site localization of the 3d orbitals due to the SIC increases the magnetic moment which is now in excellent agreement with the experimental value.

−10
A number of beyond-LSDA approaches have been attempted in the past and tried on the series of transitionmetal oxides. In Table II we listed the results of some of these calculations.
Whithin full SIC-LSDA, the results for the transitionmetal oxides are critically sensitive to the choice of orbitals taken as localized. In Ref. 33 (a LMTO-ASA calculation) the SIC-LSDA band structures are calculated in two ways: The first calculation assumes that 3d and the 2p orbitals are both localized and therefore affected by the SI, and the second is performed with only the 3d orbitals taken as localized. (In the Table we report the results of both the calculations; the d-only SIC-LSDA results are in parentheses). It is shown that the first option gives a band structure in better agreement with photoemission spectroscopy, with the exception of the band gap which is strongly overestimated. Instead, the choice of correcting only the d orbitals gives a better energy gap but also leads to a too large downshift of the d bands with respect to the p bands (e.g. for MnO and NiO the d bands lie ∼ 5-6 eV below the center of the p bands 33 ).
As a consequence, the calculation exaggerates the chargetransfer character of the band gap. The strategy of correcting both d and p orbitals represents our preferred point of view, 50 and is consistent with the pseudo-SIC which assumes that all the electron states should be corrected, without any a-priori discrimination. Indeed, the pseudo-SIC band structures are similar to those obtained in Ref. 33 with the first choice of orbitals, with the important exception that in our pseudo-SIC the energy gap is in much better agreement with experiments. We speculate that the reason for this difference is the absence in SIC-LSDA of the relaxation energy contribution, which shifts the occupied d and p bands up in energy.
The LDA+U corrects only the d states, and it shows similarities with the d-only SIC-LSDA in describing the properties of the transition-metal oxides. 32,24,51,53 Indeed in both these approaches, the transition-metal oxides are all described as charge-transfer insulators. However in LDA+U the amplitude of the d band shift is critically controlled by the parameter U. For example, in Ref. 53 (a projector augmented-wave (PAW) calculation) it is shown that for NiO the widely used value U=8 eV gives a good energy gap but a poor description of the occupied band manifold due to the exaggerated downward shift of the d bands. In contrast, the value U=5 eV gives a smaller energy gap, but optical properties and magnetic moments in better agreement with experiments. (Results for U=5 eV are reported in parentheses in Table II).
Finally, we report the results of the model GW, 54 an approach radically different from both the SIC or the LDA+U schemes, based on a model self-energy correction to the LSDA KS equations. The model GW is capable of giving an accurate description of both transition-metal oxides 55,56 and other non-magnetic semiconductors, 54,55,57 thus it should probably be considered as the most reliable reference. The fact that the pseudo-SIC values for the energy gaps are close to the results obtained within GW is an indication that, at least for this family of compounds, the SI, and not many-body effects, really represents the main source of the LSDA error.

E. Hexagonal manganites
As a final application of the pseudo-SIC, we consider hexagonal YMnO 3 , which has recently attracted the attention of both the experimental 63-66 and the theoretical 67,11 communities since it shows the uncommon characteristic of being both magnetically and ferroelectrically ordered within the same bulk phase. The study of YMnO 3 is motivated by the possible applications of ferroelectromagnetic materials as building blocks for spintronic devices 68 and, from a more fundamental standpoint, by the necessity to understand the interaction between magnetic and ferroelectric polarization and the conditions which favor this coexistence. 69,11 YMnO 3 can be grown in both orthorhombic 70 and hexagonal [71][72][73] structures, although the latter is the most stable. The orthorhombic phase, typical of manganese perovskites (e.g. LaMnO 3 ), is antiferromagnetic but not ferroelectric. The hexagonal phase shows a spontaneous ferroelectric polarization P ∼ 5.5 µC/cm 2 parallel to the c axis 74,72 below a critical temperature T c =900 K, and is A-type antiferromagnetic 75,76 at temperatures below T N =80 K. In this section we will consider only hexagonal YMnO 3 . . LSDA band structure for the A-type AFM paraelectric YMnO3. The orbital character is given in terms of the hexagonal cartesian coordinates. According to this calculation the system is a metal.
The structure of the paraelectric YMnO 3 (P6 3 /mmc symmetry) is shown in Figure 12. The Mn ions, sited on close-packed hexagonal positions, are surrounded by corner-sharing bipyramidal cages of oxygens. The unit cell is made of eight 1×1 hexagonal planes enclosing a total of 10 atoms. In the ferroelectric phase (P6 3 cm symmetry) the MnO 5 bipyramids are tilted around the axis passing through the Mn and parallel to one of the triangular base sides, thus the hexagonal planes are √ 3 × √ 3 and the unit cell has 30 atoms. However, for the purposes of illustrating the effects of the pseudo-SIC, the small oxygen rotations of the ferroelectric structure are not significant, thus in the following we only consider the paraelectric phase.
The technical features of our calculations are the same as those used in Ref. 11, where a detailed study of density of states, band structure and orbital charges within LSDA can be found. Here we focus on the remarkable changes in the band structure of this material produced by the pseudo-SIC. An obvious requirement for a ferroelectric material is that it must be insulating, but the LSDA calculations describes YMnO 3 as a metal, as can be seen in Figure 13. (The hexagonal crystal field splitting and the cartesian coordinates are the same as those described previously for MnO and NiO.) In hexagonal symmetry the four d electrons of the Mn 3+ ion entirely occupy the two orbital doublets (d ↑ xy , d ↑ x 2 −y 2 ) and (d ↑ xz , d ↑ yz ), leaving the d ↑ z 2 orbital, which is the highest in energy, empty. This ordering causes a magnetic moment on Mn equal to 3.8 µ B , slightly lower than its nominal value 4 µ B because of Mn d-O p hybridization.
However the LSDA crystal field splitting is too small to open a gap, and one band crosses E F in the direction Γ-Q-P which is parallel to the k z =0 plane. This band comes from a mix of Mn d ↑ xy , d ↑ x 2 −y 2 orbitals and p ↑ x , p ↑ y orbitals from the oxygens lying in-plane with the Mn (i.e. it is a pdσ hybridization). In contrast, the band energies are flat along [0001], as is typical for strongly layered compounds. (Notice that according to our LSDA calculations the system is metallic even within the ferroelectric phase.) The metallicity is a fatal shortcoming since it precludes the possibility of accessing any ferroelectric properties, such as spontaneous polarization, Born effective charges, and piezoelectric tensors. (according to our calculations the system is still metallic within the ferroelectric structure) The pseudo-SIC ( Figure 14) repairs the fault: a gap, E g = 1. 40  On the basis of the pseudo-SIC calculation which correctly describe the YMnO 3 as an insulator we are now able to study the ferroelectric structural displacements and other properties related to the electric polarization. This will be the subject of future work.

V. DISCUSSION AND SUMMARY
In this paper we have proposed and implemented an innovative first-principles approach which is able to improve the LDA and LSDA description of the electronic properties for a vast range of compounds, including strongly-correlated and magnetic systems, and requires only a minor increase of computing cost. Our strategy, based on the approximation of the true SIC potential with a pseudopotential-like projector, was built primarily on the work of Ref. 29, where it was shown that the inclusion of the atomic SIC within the pseudopotential construction can consistently improve the agreement of band structures with photoemission spectra for a range of strongly ionic II-VI and III-V insulators. Working on this original idea, we defined a pseudo-SIC functional explicitly dependent on the orbital occupation numbers which are self-consistently calculated from the Bloch wavefunctions and represent the natural extension to periodic systems of the atomic occupation numbers. They can be fractional due to charge hybridization, eigenvalue degeneracy, or Fermi-Dirac distribution.
We tested the pseudo-SIC on three different classes of compounds: non-magnetic semiconductors, magnetic insulators and ferroelectrics, and found a generally good (sometimes very good) agreement with experimental data. In particular, the size and the orbital characterization of the fundamental energy gap and the local magnetic moments are much better described than in LSDA.
We expect the same kind of improvement in any system where the electron localization plays a major role. Examples of systems which could benefit from the pseudo-SIC application are innumerable, including bulk systems with d and/or f electrons, Jahn-Teller and orbitally ordered systems, and non-bulk systems like defects and impurities, surface and interface states, bulk resonances and core states. The SIC may affect a vast range of phenomena, encompassing, among others, surface reconstructions, adsorbtion and diffusion of atoms and molecules on surfaces, doping in semiconductors, alloys, and homoand hetero-junctions.
Of course, the discrepancies between theoretical and experimental results which are not attributable to the SI should not be expected to disappear within pseudo-SIC. A typical example is the fundamental energy gap of bulk Si, whose LDA value is ∼ 0.5 eV lower than the experimental gap, due to genuine self-energy effects which are outside the realm of the DFT itself. As a consistency test, we calculated the band structure of bulk Si within pseudo-SIC. Happily, it comes out to be very similar to the LDA band structure. This is a consequence of the fact that each band (either occupied or empty) shares the same orbital character (sp 3 -hybridization), thus the pseudo-SIC mostly produces a trivial shift of the whole band manifold.
Comparing it to other corrective methodologies available in the literature, the pseudo-SIC is closest to the spirit of the LDA+U: in both cases the LSDA energy functional is augmented by a term depending on the orbital occupation numbers. In the pseudo-SIC the modification is milder, and consists of subtracting out a contribution (i.e. the SI) which is already present in the LSDA energy functional. The LDA+U is a more radical departure from the LSDA, since it rewrites the whole Coulomb electron-electron interaction in terms of local orbitals according to the multiband Hubbard expression (which is, in itself, self-interaction free). Based on the undeniably limited set of results presented in this paper, the pseudo-SIC seems to be at least as accurate as the LDA+U, but a fair comparison will need a much larger body of work.
Finally, we point out that our method suffers a drawback, that is the non-variationality of the energy functional. This is especially annoying when force or stress calculation is required, since the familiar Hellman-Feynman formulations cannot be applied and additional contributions arise due to the first-order change of the wavefunctions. We will address these problems and the studies of forces and stresses in future publications.

ACKNOWLEDGMENTS
We are indebted to Gerhard Theurich for many thoughtful discussions, and to David Vanderbilt for his critical reading of the manuscript. We acknowledge financial support from the National Science Foundation's Division of Materials Research under grant number DMR 9973076. Also, this work made use of MRL Central Facilities supported by the National Science Foundation under award No. DMR00-80034. Most of the calculations have been carried out on the IBM SP2 machine of the MHPCC Supercomputing Center in Maui, HI.

A. Review of the USPP formulation
In the USPP approach 34 the constraint of normconservation on the pseudo wavefunctions is relieved, so that the pseudo wavefunctions and the associated pseu-dopotentials can be extremely smooth (i.e. "ultrasoft"). In order to restore the correctly normalized electron charge, the charge density is written as: Here α = [n, l, m, R] and α ′ = [n ′ , l ′ , m ′ , R] are sets of orbital quantum numbers and atomic positions R, and β α (r) are the USPP projector functions. 34 The first term within the square brackets is the ultrasoft charge, and the second term the "augmented" charge, that is the portion of the valence charge which is localized within the atomic core radii and restores the normalization of the total charge. The atomic charges Q αα ′ are, by construction, where φ AE α and φ P S β are atomic all-electron and pseudo wavefunctions, respectively. The release of norm-conservation leads to the generalized Kohn-Sham equations: where V LOC (r) is the local part of the pseudopotential, D σ αα ′ is the non-local part andŜ is the overlap matrix which generalizes the orthonormality condition: Here q αα ′ are the integrals of the augmented charges Q αα ′ (r), and < ψ nk |Ŝ|ψ n ′ k ′ >= δ n,n ′ δ k,k ′ .
Finally, the non-local pseudopotential projector is made up of two contributions: D σ αα ′ =D αα ′ + dr (V LOC (r) + V σ HXC (r)) Q αα ′ (r). (22) The first term on the right side of Eq.22 is the usual Kleinman-Bylander projector, and contributes to the 'bare' pseudopotential (i.e. it is calculated within the atomic reference configuration). The second term is specific to the USPP formalism, and represents the action which the local and screening potentials exert on the augmented charges. Since this term depends on V σ HXC , it has to be updated during the self-consistency cycle.
In order to ensure better transferability, two atomic reference states, corresponding to different energy values, are usually included in the projector for each angular quantum number. As a consequence, the USPP projector contains non-diagonal terms (α,α ′ ), where α = (ν l , l, m), and α ′ = (ν ′ l , l, m), and ν l , ν l '= 1,2. The atomic reference eigenstates do not need to correspond to bound, normalized solutions of the free atom, but may be unphysical eigenstates of the Schrödinger equation, useful to extend the pseudopotential transferability into a larger energy range. Thus the atomic pseudo wavefunctions may diverge at large r, but the projector functions β α (r) and the matrixD αα ′ are always short-ranged and well-defined by construction.

B. Pseudo-SIC within USPP
The USPP implementation of the pseudo-SIC requires some generalization of the formalism described in Sections III A and III B. The charge densities n σ i of the (pseudo) atomic orbitals φ i are: and the occupation numbers p σ i become: Furthermore, at variance with the norm-conserving case, the non-local part of the USPP depends selfconsistently on the screening potential itself (see Eq.22). As a result it also must be self-interaction corrected. The SI part of the non-local USPP is given by: Thus, the pseudo-SIC KS equations finally are: whereV σ SIC is given by Eqs.5,12, and 13. The total energy is the same as that given in Eq.17, except for the fifth term which now is: