First-principles modeling of localized d states with the GW @ LDA + U approach

Hong Jiang,1 Ricardo I. Gomez-Abal,2 Patrick Rinke,2 and Matthias Scheffler2 1Beijing National Laboratory for Molecular Sciences, National Laboratory of Rare Earth Material Chemistry and Application, Institute of Theoretical and Computational Chemistry, College of Chemistry and Molecular Engineering, Peking University, 100871 Beijing, China 2Fritz-Haber-Institut der Max-Planck-Gesellschaft, 14195 Berlin, Germany Received 20 April 2010; revised manuscript received 19 June 2010; published 14 July 2010


I. INTRODUCTION
Kohn-Sham ͑KS͒ density-functional theory ͑DFT͒ ͑Refs. 1 and 2͒ in the local-density or generalized gradient approximation ͑LDA/GGA͒ to the exchange-correlation ͑xc͒ energy functional has become "the standard approach" for firstprinciples electronic-structure calculations of extended systems. 3 Although the Slater-Janak 4,5 transition state theorem relates KS DFT single-particle energies at half occupation to excitation energies, KS eigenvalues of a ground-state calculation are often used to interpret excited-state properties as, for example, probed by direct and inverse photoemission spectroscopy ͑PES/IPES͒ or optical absorption. This practice, however, should be exercised with caution. Even for weakly correlated systems such as sp semiconductors, KS-LDA/GGA band gaps are considerably underestimated ͑by about 20-50 % compared to experiment͒. 6 The problem can be more severe for systems with open d or f shells, often called strongly correlated systems, for which certain wide gap insulators are predicted to be metallic. For these systems, even ground-state properties such as the magnetic ordering might be qualitatively wrong. This is best illustrated in the well-known failure of LDA/GGA for the later transitionmetal oxides. 7 For the quantitative description of quasiparticle excitations in solids as measured by PES/IPES many-body perturbation theory in Hedin's GW approximation 8 has become the method of choice. Having earned its merits for sp bonded systems, GW's application to d-or f-electron systems is relatively recent and not as successful as anticipated.  Applied in the standard way as perturbation to an LDA ground state ͑G 0 W 0 @ LDA͒ the GW calculation suffers from the pathologies of the LDA starting point. [14][15][16][17][18]20,21,[24][25][26]29,[31][32][33][34][35] Following our previous work for f-electron systems, 29 we demonstrate in this paper that applying Hubbard U corrections to the LDA calculations ͑LDA+ U͒ provides an insightful way to systematically analyze the problem. We investigate examples from three common classes of semiconductors: empty d states ͑ScN͒, fully filled semicore d states ͑ZnS͒, and partially filled d states ͑transition-metal oxides NiO, MnO, FeO, and CoO͒. For physically meaningful values of U, G 0 W 0 calculations based on LDA+ U ground states ͑G 0 W 0 @ LDA+ U͒ give a balanced description of both the itinerant and localized states, which is consistent with our previous findings when applying the same approach to f-electron systems. 29 Compared to 4f-electron systems like lanthanide oxides, where highly localized 4f states barely participate in the chemical bonding, typical d-electron systems exhibit a much stronger interaction between localized d states and itinerant states. It will therefore be elucidating to apply the G 0 W 0 @ LDA+ U approach to these systems.
The GW approximation is the first-order term in a systematic expansion of the self-energy ͑⌺͒. The exact ⌺ could in principle be obtained by solving a set of integrodifferential equations. 8 In practice, this is still beyond the reach of today's computational machinery even for the simplest systems like the homogeneous electron gas. Like this exact set, the GW equations could also be solved self-consistently. However, the omission of the vertex function that would introduce higher order interactions with every iteration makes a self-consistent solution inconsistent and worsens the spectral properties of the Green's function. 6,35 In practical calcula-tions, a "best G best W" strategy is therefore adopted, where the GW quasiparticle energies are expressed as a first-order correction to a reference single-particle Hamiltonian H 0 . Both G and W are calculated using eigenenergies and eigenfunctions of H 0 , hence the abbreviation G 0 W 0 .
Early GW studies of late transition-metal oxides, NiO and MnO, showed that G 0 W 0 calculations based on LDA singleparticle energies and wave functions only slightly improve over the LDA description. However, introducing an approximate level of self-consistency leads to better agreement with experiment. 10,11,13 This indicates that the main difficulty for d / f-electron systems may not come from GW itself, but from the failure of the LDA/GGA as a starting point, although partially occupied d or f states may indeed require the inclusion of higher order correlation effects that go beyond the GW approach. Before venturing into the realm of higher order correlation effects ͓by, e.g., combining GW with dynamical mean-field theory ͑DMFT͒ ͑Refs. 36 and 37͒ or through vertex corrections 38,39 ͔ it is important to establish the limitations of the GW approach.
In the spirit of the best G best W approach, much work has been invested recently to develop better single-particle reference Hamiltonians on which to base G 0 W 0 quasiparticle energy calculations. 15,21,[24][25][26]32,[40][41][42] In this work we follow the same strategy. For d / f-electron systems, a simple and effective approach to overcome the major failure of LDA/ GGA is the LDA+ U method in which the LDA total energy is augmented by a local Hubbard correction, characterized by the on-site Coulomb interaction U. [43][44][45][46] Since the correction is only applied to a subset of states, LDA+ U itself is not expected to provide a quantitatively accurate description of the whole band structure. It can, however, serve as a reasonable starting point for G 0 W 0 calculations. 25,29,47 In a previous study, 29 we have applied the G 0 W 0 @ LDA+ U approach to f-electron systems, using lanthanide oxides as examples, and found that both localized and itinerant states are described quite accurately and in fact superior to state-of-the-art DMFT. To further assess the performance of the G 0 W 0 @ LDA+ U method in a critical manner, we investigate a series of prototypical 3d-electron systems: ͑1͒ ScN with empty d states, ͑2͒ ZnS with shallow semicore d states, and ͑3͒ transition-metal oxides ͑MnO, FeO, CoO, and NiO͒ with partially occupied d states. Although similar conclusion can be reached, 3d-electron systems also show some significant differences compared to 4f systems, and in some sense pose a more serious challenge for first-principles many-body theory calculations.
The paper is organized as follows. In the next section, the main theoretical framework is outlined and certain aspects of our implementation are discussed. In Sec. III, results for three type of representative systems are presented. Section IV summarizes our work.

A. GW approximation for quasiparticle excitations
Hedin's equations 8 for the Green's function G, the polarizability P, the screened ͑W͒ and bare ͑v͒ Coulomb interaction, the self-energy ⌺, and the vertex function ⌫ ⌫͑1,2,3͒ = ␦͑1,2͒␦͑1,3͒ G͑4,6͒G͑7,5͒⌫͑6,7,3͒d͑4,5,6,7͒ ͑4͒ form a closed set of integrodifferential equations. In Eqs. ͑1͒-͑4͒ we adopted the short-hand form 1 ϵ͑r 1 , 1 , t 1 ͒ to denote a triple of space, spin, and time variables. We will also use x = ͑r , ͒ to denote the collective space and spin coordinate. Accordingly ͐d͑1͒ is a shorthand notation for the integration in all three variables of the triple. By means of Dyson's equation which links the noninteracting system with Green's function G 0 to the fully interacting one ͑G͒ via the self-energy ⌺, Hedin's equations could, in principle, be solved selfconsistently starting from a given G 0 .
The GW approximation is formally obtained by retaining only the zeroth-order term in the vertex function ͓⌫͑1,2,3͒ = ␦͑1,2͒␦͑1,3͔͒. A few attempts to solve this reduced set of equations fully self-consistently have been made in the past. 37,[48][49][50][51][52][53][54][55][56][57][58][59] Although full or partially self-consistent GW calculations can give accurate ground-state total energies and are essential for particle number conservation, 60,61 the quasiparticle properties deteriorate. [49][50][51] This is a result of the successive introduction of higher order electron-electron interaction terms of certain type that are not balanced by other higher order terms contained in the vertex function.
The most widely used procedure to solve the GW equations in practice is the so-called G 0 W 0 approach in which G and W are calculated from the eigenenergies ͕⑀ nk ͖ and wave functions ͕ nk ͖ of a single-particle reference Hamiltonian, H 0 . The self-energy then takes the form where is an infinitesimal positive number. 62 follows from the polarizability

͑10͒
The QP energies E nk are then calculated by first-order perturbation theory, treating ␦⌺ ϵ ⌺ − V xc as the perturbation where V xc ͑r͒ is the exchange-correlation potential already included in H 0 . Equation ͑11͒ is often linearized in energy where the QP renormalization factor Z nk is given by For sp semiconductors it has been demonstrated that further improvement can be obtained without introducing too much computational overhead by partial self-consistency. 6,63 In the so-called energy-only self-consistent or GW 0 approach the energy denominator in the Green's function is updated by E nk 's but W remains unchanged.
Recognizing that LDA or GGA are not always the best reference Hamiltonian for a G 0 W 0 calculation several alternatives have been proposed. 22,29,32,34,41,42,47,64,65 They roughly fall into two categories. The first comprises different exchange-correlation functionals in H 0 , e.g., exact exchange in the optimized effective potential ͑OEPx͒ approach, 32,66 hybrid functionals, 34,65 and LDA+ U. 25,29,47 All these variants do not increase the computational cost of the G 0 W 0 calculation ͑although some xc functionals might be more computationally expensive than others in the ground-state calcula-tion͒ but their adequacy is not known a priori and has to be carefully investigated on a case by case basis. The second category includes several variants of approximate selfconsistency in GW. 15,22,30,41,42,64 These schemes appeal because they remove the dependence of G 0 W 0 on the starting point but their nonuniqueness remains a disconcerting factor.

B. LDA+ U method
For systems with d / f shells, the severe self-interaction error of LDA or GGA often results in an inadequate description. A simple and effective approach to correct for this is to introduce a local, Hubbard-type correction ͑LDA+ U͒, characterized by the on-site Coulomb ͑U͒ and the exchange in-teraction ͑J͒. 44,46,67 In its most general form, the LDA+ U total energy is written as where = ↑ , ↓ is the spin index ͑from now on we will write out the spin degree of freedom explicitly, assuming a collinear spin polarization͒. n is the local-density matrix defined as where f nk denotes the occupation number of the state nk and ͕͉m͘ϵ͉In p lm͖͘ denote a set of atomiclike local orbitals on the Ith atom with the principle, angular, and magnetic quantum numbers n p , l, and m, respectively. n ϵ Tr n is the local occupation number. E ee contains the electron-electron interaction of the localized electrons and the double-counting term Ē ee removes the part that was already included in the LDA Hamiltonian. The LDA+ U approach is obtained by treating E ee in a Hartree-Fock-type fashion 45 with an effective screened Coulomb interaction V ee . Using the angular expansion of V ee the Coulomb matrix element ͗m 1 , m 2 ͉V ee ͉m 3 , m 4 ͘ can be expanded as follows: ͗m 1 ,m 2 ͉V ee ͉m 3 ,m 4 ͘ = ͚ L F L C L ͑m 1 ,m 2 ,m 3 ,m 4 ͒. ͑18͒ F L are radial Coulomb integrals ͑Slater's integrals͒ F L ϵ ͵ dr ͵ drЈr 2 rЈ 2 ͉R n p l ͑r͉͒ 2 v L ͑r,rЈ͉͒R n p l ͑rЈ͉͒ 2 ͑19͒ and C L ͑m 1 , m 2 , m 3 , m 4 ͒ are angular integrals, which, using Wigner's three-j symbols, 68 read For d electrons, only F 0 , F 2 , and F 4 are nonvanishing, and they are related to U and J by U = F 0 and J = ͑F 2 + F 4 ͒ / 14. By fixing the ratio F 4 / F 2 , which is nearly constant ͑ϳ0.625͒ in free atoms, 69 one can use U and J as only parameters to determine E ee . 45,46 The double-counting correction term Ē ee ͓n ͔, on the other hand, is more arbitrary and is one of the largest problems in the LDA+ U approach. 44,45,67,[70][71][72][73] Most frequently Ē ee is taken as the following function of the local occupation numbers n : which can be obtained from Eq. ͑16͒ by neglecting orbital polarization effects, often called fully localized limit ͑FLL͒. 70,74 The single-particle Hamiltonian corresponding to E LDA+U ͓͑r͔͒ reads ͗mЈ͉ is a site-and orbitaldependent nonlocal potential arising from the LDA+ U correction term, Neglecting the anisotropy of the local Coulomb interaction, i.e., dropping all L Ͼ 0 terms in Eq. ͑17͒ and further using the identity C L=0 ͑m 1 , m 2 , m 3 , m 4 ͒ = ␦ m 1 ,m 3 ␦ m 2 ,m 4 , 68 we obtain

͑24͒
Within this approximation we have

͑25͒
If we define the local projection so that the on-site density matrix is diagonal, ␦V takes the simple form

͑26͒
The main physical effect of ␦V is therefore to push occupied localized states down and unoccupied ones up in energy, which effectively opens a gap that might have been absent in the LDA description.

C. LDA+ U as an approximation to GW
For highly localized d / f states LDA+ U can be viewed as an approximate GW scheme, as first pointed out by Anisimov et al. 46 The original derivation is fairly involved and based on several specific assumptions that turn out to be not nec-essary. We present here a derivation that is simpler and more general, starting from the static Coulomb-hole and screened exchange ͑COHSEX͒ approximation 6,8 to the GW selfenergy. The COHSEX approximation is obtained by omitting the dynamic features of the screened Coulomb interaction Using the closure relation of the KS states we obtain The matrix elements of ⌺ ͑r , rЈ͒ projected onto the local subspace ͕͉m͖͘ can therefore be written as We now decompose the Kohn-Sham wave functions according to with C m;nk ϵ͗ m ͉ nk ͘. ͉ nk ͘ can be regarded as a "pure" itinerant state in which the contribution of localized states has been projected out. Inserting Eq. ͑31͒ into Eq. ͑30͒ and further assuming that any matrix element that involves the overlap of m ͑r͒ and nk ͑r͒ can be neglected, we obtain where we have used the relations

͑33͒
In the spirit of the LDA+ U approach, we consider only corrections to localized states where V H loc is the Hartree potential of the localized electron density

͑36͒
V LDA loc, ͑r͒ϵ␦Ē ee loc / ␦ ͑r͒ is the interaction potential among the localized states in the LDA. Using a similar argument as for the FLL approximation to the double-counting correction, i.e. Eq. ͑21͒, we obtain where F 0 ͑0͒ and J ͑0͒ are the first Slater integral and the on-site exchange term from the bare Coulomb interaction v͑r , rЈ͒, respectively.
We therefore have By neglecting the anisotropy in both the bare and screened Coulomb interaction, i.e., using the approximation in Eq. ͑24͒, we obtain exactly the same expression as in LDA+ U ͓Eq. ͑25͔͒ where U is identified as the first Slater integral F 0 arising from the static screened Coulomb interaction W͑r , rЈ ;0͒. 46 To summarize, LDA+ U follows from the GW approach under the assumption that: ͑1͒ The frequency dependence of the screened Coulomb interaction is neglected; ͑2͒ quasiparticle corrections are only applied to localized states, whereas itinerant states are still treated at the LDA level; ͑3͒ all Coulomb matrix elements that involve the overlap of a localized state and an itinerant state are neglected ͑which is equivalent to omitting the many-body exchange interaction between localized and itinerant electrons͒.
None of these assumptions are of course fully satisfied in realistic systems. The dynamic character of the screened Coulomb interaction is actually stronger for localized electrons than for itinerant ones, as demonstrated by the fact that the renormalization factor ͓Eq. ͑13͔͒ typically takes a value of ϳ0.5-0.6 for localized states, in contrast to the typical value of ϳ0.8-0.9 for itinerant states. 24 The LDA description of the itinerant states suffers from the band-gap problem and the coupling between localized and itinerant states is critical for the physical and chemical properties of d / f-electron systems.

D. GW @LDA+U method
The preceding discussion elucidates that LDA+ U is a relatively crude approximation to GW and is therefore not necessarily expected to provide an accurate description of d or f-electron systems on its own. However, since LDA+ U corrects the major failure of LDA for localized systems, it is likely to serve as a good starting point for G 0 W 0 and GW 0 ͑i.e. energy-only self-consistent GW with fixed W 0 ͒ calculations. Formally, the only difference between the LDA-and

͔͘. ͑40͒
This equation illustrates that the double-counting term contained in ⑀ nk is exactly canceled out in the GW @ LDA+ U approach, which is a significant advantage of the scheme. Before we introduce the specific values of U for the systems in this work, we will treat U as a parameter and analyze the U dependence of GW @ LDA+ U, which is quite different from that of LDA+ U. The ␦V contribution in the LDA+ U single-particle energies can be split off and ⑀ nk can be regarded as LDA single-particle energies calculated with LDA+ U wave functions. If we now expand the self-energy around ⑀ nk instead of ⑀ nk , under the assumption that ⌺ nk ͑E͒ can be approximated by a linear function around ⑀ nk we obtain, after some simple algebraic manipulations, has no explicit U dependence and the quasiparticle energies in GW @ LDA+ U therefore depend only implicitly on U. The two main factors are ͑1͒ the difference in wave functions between LDA+ U and LDA, and ͑2͒ the change in screening introduced by changes in the singleparticle spectrum, most notably the band gap. Regarding point ͑1͒, occupied states with strong d / f character will be pushed toward lower energies and become more localized, whereas unoccupied states with strong d / f character are pushed to higher energies and become more delocalized. This implies that ⑀ nk in Eq. ͑42͒ is larger for more localized states. The U dependence of the second term in Eq. ͑42͒ is determined by the change in screening, which can be characterized by the macroscopic static dielectric constant, Increasing U usually increases the band gap and therefore reduces ⑀ M . The limit of vanishing screening gives the Hartree-Fock self-energy, whose corresponding band gap is dramatically overestimated. The second term therefore has a tendency to increase the GW @ LDA+ U band gap with increasing U. In practice both terms can significantly shape the U dependence of the system under study, as we will demonstrate in the next section.

E. Computational details
All DFT calculations are performed using the WIEN2K package 75 in which the Kohn-Sham equations are solved in the full-potential ͑linearized͒ augmented plane wave plus local orbital ͓FP− ͑L͒APW+ lo͔ approach. 76 The following parameters for the FP− ͑L͒APW+ lo basis are used: muffin-tin ͑MT͒ radii R MT ͑in units of Bohr͒ are ͑2.13, 1.89͒ for ScN, ͑2.05, 1.95͒ for ZnS, ͑2.10, 1.86͒ for MnO, ͑2.10, 1.77͒ for FeO, ͑2.05, 1.75͒ for CoO, and ͑1.97, 1.75͒ for NiO; wave functions are expanded by spherical harmonics with l up to l max = 10 in the MT spheres, and by plane waves with the energy cutoff determined by min R MT ϫ K max = 7.0 in the interstitial ͑IS͒ region; the potential and electron density are expanded in cubic harmonics l up to l max = 4 within the MT spheres, and by plane waves in the IS region. For the LDA+ U calculations, the double-counting correction is treated in the FLL scheme ͓i.e., Eq. ͑21͒, LDA+ U, selfinteraction correction ͑SIC͒ in the WIEN2K notation͔. 44,77 GW calculations were performed using the FHI-gap ͑Green's function with augmented plane waves͒ package, a recently developed all-electron GW add on to WIEN2K. 29,78-80 All important convergence parameters have been monitored to achieve an overall numerical accuracy of Ӎ0.05 eV. The Brillioun zone was sampled with 4 ϫ 4 ϫ 4 k meshes; unoccupied states with energy up to 136 eV were taken into account. Densities of states ͑DOSs͒ are calculated using ϳ1000 k points in the Brillouin zone. In the GW case, DOS are obtained from GW quasiparticle energies, first calculated on the sparse k mesh and then interpolated to the fine k mesh ͑ϳ1000͒ using the Fourier interpolation technique. 81 When the calculated DOS is compared to the experimental spectral data, a Gaussian broadening of 0.6 eV was chosen to mimic typical experimental resolutions.
We use the experimental lattice constants for all materials considered in the next section. ScN has a NaCl structure with a͑exp͒ = 4.50 Å. 82 ZnS takes the zinc-blende structure with a͑exp͒ = 5.42 Å taken from Ref. 63. The transition-metal monoxides ͑MO, M = Mn, Fe, Co, and Ni͒ crystallize in the NaCl structure in the paramagnetic phase. In this work, however, we consider the type-II antiferromagnetic ͑AFM-II͒ phase, which is the most stable structure below the Neel temperature. 7

A. Determination of U
The parameters U and J describe the effective Coulomb and exchange interaction between localized electrons in a crystalline environment. The determination of U from first principles has attracted growing interest in recent years. 67,74,[84][85][86][87][88][89][90][91][92] The most widely used approach is the constrained DFT formalism, which was developed in a general form by Dederichs et al. 86 and later widely used to calculate parameters for effective Hamiltonians ͑Hubbard or Anderson͒. 93,94 The Hubbard U is obtained from the totalenergy variation for different integer occupations of a fictitious atom ͑usually denoted the impurity͒ embedded in a crystalline environment via the supercell technique. In practice, a constrained DFT calculation proceeds along the following lines ͑using NiO as an example͒: ͑1͒ in a supercell of NiO one Ni atom is treated as an "impurity;" ͑2͒ a suitable constraint ͑see below͒ is imposed that allows variations in the local occupation number of the impurity 3d states ͑n d ͒ while keeping the system as a whole neutral ͑by either changing the number of valence electrons correspondingly or adding a uniform compensating charge͒; ͑3͒ from the corresponding total energy E͑n d ↑ , n d ↓ ͒ for different occupations U and J are determined according to There are different ways to constrain n d , depending on the definition of the local projection and the implementation for the Kohn-Sham equations. In Anisimov and Gunnarsson's original work, 67 the 3d electrons were constrained by setting the hopping integrals between the impurity 3d orbitals and other orbitals to zero in their linear muffin-tin orbital ͑LMTO͒ code based on the atomic sphere approximation ͑ASA͒. In the linearized augmented plane wave ͑LAPW͒ framework as implemented in WIEN2K, the constraint is imposed by removing the d orbitals on the impurity from the LAPW basis and treating the impurity 3d states as core states with fixed occupation number. 95 Alternatively one can also constrain n d to the desired value by solving a constrained optimization problem 84,91 E͑n d ͒ ϵ min This approach is more general and applicable to any code as long as the local projection in Eq. ͑15͒ can be defined. The values of U obtained from constrained DFT can depend significantly on the implementation of the constraint and are additionally basis set specific. This explains the spread of U's reported in the literature for the same material by different groups. However, as long as a consistent approach is adopted at all stages, i.e., the same basis set and local projectors are used in the constrained DFT and subsequent LDA+ U and G 0 W 0 @ LDA+ U calculations, as in our work, the U values are meaningful.
For transition-metal oxides, the U's obtained in the LMTO-ASA based constrained DFT calculations by Anisimov, Zaanen, and Andersen ͑AZA͒ ͑Ref. 43͒ have been widely applied. 83 However, several reports in the literature indicate that these U's are too large. 71,96,97 More recent studies employing the linear-response approach 84,85 or the constrained random-phase approximation ͑RPA͒ ͑Refs. 90 and 98͒ yield significantly smaller U's. In this work, U's and J's were calculated with constrained DFT as implemented in the LAPW basis. 95 Our values of U are listed in Table I. They are significantly smaller than AZA's but very close to the linear response and constrained RPA ones. 84,85,92 The difference between our and AZA's constrained DFT results can be attributed to the different basis used in the calculations, i.e., full-potential vs atomic sphere approximation, and LAPW vs LMTO.
In practical LDA+ U calculations, U and J are often combined by redefining U as U eff = U − J and setting J =0. 71 We found that such a simplification can have noticeable effects on the resultant electronic properties, as also recently reported by Ylvisaker and Pickett. 73 We therefore only use J = 0 when investigating the U dependence but use the finite J's as obtained by Eq. ͑45͒ and listed in Table I when we compare our G 0 W 0 @ LDA+ U results to experiment.

B. ScN: Empty d states
The 3d states in ScN form the lowest conduction bands and apart from a small hybridization with the valence bands are nearly empty. 99,100 Interest in ScN has recently increased both experimentally and theoretically due to its possible application in optoelectronic devices. [100][101][102][103][104] LDA predicts ScN to be a semimetal with a negative band gap of Ϫ0.14 eV. Recent calculations based on the exact-exchange OEPx ͑Ref. 101͒ or the screened exchange approach 99 Figure 1 shows the indirect ͑E g ⌫-X ͒ and direct ͑E g X-X and E g ⌫-⌫ ͒ band gaps of ScN in LDA+ U, G 0 W 0 , and GW 0 @ LDA+ U as a function of U. In LDA+ U, ScN becomes insulating when U is larger than ϳ2.0 eV and the gap increases almost linearly with increasing U. To obtain the experimental band gap of about 0.9 eV, the Hubbard U would have to be as large as 10 eV, much larger than the U − J = 2.9 eV obtained from constrained DFT calculations and clearly unphysical. Conversely, in G 0 W 0 and GW 0 @ LDA+ U the band gaps decrease with increasing U.
With the exception of very large U͑Ͼ ϳ 9 eV͒ the GW 0 band gaps are always larger than G 0 W 0 ones, similar to what is usually observed for sp semiconductors. 63 This somewhat peculiar U dependence of GW @ LDA + U can be understood by analyzing the effects of U on the single-particle energies and wave functions from LDA+ U using Eq. ͑42͒, which is shown for the direct gap at the ⌫ point in Fig. 2. To cross-check that Eq. ͑42͒ is valid in this case we have verified that it gives nearly the same band gaps as Eq. ͑40͒ ͑circles and crosses in the right panel of Fig. 2͒. With increasing U, the two terms in Eq. ͑42͒ exhibit the opposite behavior. ⑀ nk LDA , the LDA single-particle energies calculated with LDA+ U wave functions, decreases, because the bottom of the conduction, which is mainly of Sc 3d character, is pushed upward strongly by the U-dependent correction term and becomes more delocalized, whereas the N 2p character leaves the top of the valence band largely independent of U. Conversely, the second term increases, because a larger LDA+ U band-gap results in weaker screening, as demonstrated by the decreasing dielectric constant ͑right panel of Fig. 2͒. As a result, the self-energy in the second term of Eq. ͑42͒ increases the band gap with increasing U.  , ͑b͒ E g X-X , and ͑c͒ E g ⌫-⌫ .
FIRST-PRINCIPLES MODELING OF LOCALIZED d… PHYSICAL REVIEW B 82, 045108 ͑2010͒

045108-7
We note that the U dependence in ScN is similar to that of the empty 4f-shell compound CeO 2 . 29 The band gap in CeO 2 is formed between valence-band states of predominantly oxygen 2p character and unoccupied 4f states that lie below conduction-band states of predominantly 5d character. The p-f gap depends strongly on U in LDA+ U but is nearly constant in G 0 W 0 @ LDA+ U. 29 In Table II different theoretical band gaps of ScN are compared to experiment. The G 0 W 0 @ LDA band gaps were obtained by extrapolating the G 0 W 0 @ LDA+ U results to the U − J → 0 limit via a quadratic fitting function. We verified that this gives nearly the same band gaps ͑0.94 vs 0.87 eV͒ as extrapolating the G 0 W 0 @ LDA band gaps for different lattice constants, following the same procedure as in Ref. 100. We further note that the G 0 W 0 @ LDA band gaps for ScN in this work differ from those in Ref. 100, where a pseudopotential plane-wave implementation was used. The discrepancy of ϳ0.2 eV is likely due to the core-valence exchange-correlation partitioning and pseudoization error of the pseudopotential GW approach. 78 In general our GW results are in good agreement with experiment, whereby the best agreement is achieved by GW 0 @ LDA+ U.

C. ZnS: Semicore d states
It has been recognized that the G 0 W 0 @ LDA approach performs worse for systems with shallow semicore d states, e.g., II B -VI compounds or group-III nitrides, than for other sp semiconductors. The actual reason is still a matter of debate but it is likely that it is associated to the absence of strong short-range correlations in the GW approximation. 35,42,63,[105][106][107][108][109][110] The error in the band gap can be considerable and the binding energies of semicore d electrons ͑⑀ d ͒ are significantly underestimated. Miyake et al. 47 have recently studied ZnS in G 0 W 0 @ LDA+ U with a particular emphasis on the d-band position. They found that ͑1͒ using U = 8 eV and J = 1 eV, LDA+ U reproduces the experimental d-band position but the band gap is still significantly underestimated; ͑2͒ G 0 W 0 @ LDA+ U pushes the d bands back to their G 0 W 0 @ LDA values. The same observation was made by Shishkin and Kresse. 63 In a more recent study by the same authors 42 it was shown that partial selfconsistency in GW improves the description of band gaps for II B -VI compounds but the error in ⑀ d remains large even when both self-consistency and an improved dielectric function in W are taken into account. 42 It was further argued that a more accurate description for the d-electron binding energies will require vertex corrections in both W and the selfenergy.
Our own LDA+ U based G 0 W 0 and GW 0 calculations for ZnS are summarized in Fig. 3 and confirm these observations. Since the semicore d states are fully occupied and fall between the S 2p and S 2s derived bands, the effect of U on the band gap is indirect and very weak. The valence and conduction bands change only slightly for U's from 0 to 11 eV, and the band gap increases by 0.42 eV, 0.28 eV, and 0.24 eV in LDA+ U, G 0 W 0 , and GW 0 , respectively. This is consistent with the fact that the orbital-dependent potential ␦V is applied only to the d states and influences the band gap mainly through the p-d coupling. 111 With increasing U − J, ⑀ d increases linearly in LDA+ U and reaches the experimental value at U ϳ 10 eV, which is considerably higher than our U − J value of ϳ6.0 eV obtained from constrained DFT. In addition, once the GW corrections are added to the LDA+ U energies, the G 0 W 0 and GW 0 results become very similar to the respective LDA-based GW calculations and exhibit only a very weak variation with U. The inset of Fig.  3 shows the RPA dielectric constant ⑀ M , which decreases as U increases. But even at U as large as 10.9 eV, ⑀ M is still significantly overestimated when compared to experiment, which is consistent with the underestimation of the band gap in G 0 W 0 and GW 0 .

D. Transition-metal monoxides
The transition-metal monoxides MnO, FeO, CoO, and NiO are frequently considered to be prototypical strongly  correlated electron systems and have become testbeds for many new first-principles approaches. 10,15,21,25,43,83,96,[112][113][114][115][116] The four compounds roughly fall into two classes: ͑1͒ MnO and NiO, for which standard band theory within LDA or GGA still predicts an insulating ground state when applied in the AFM-II phase but with significantly underestimated band gaps; ͑2͒ FeO and CoO, for which LDA and GGA give metallic ground states even in the AFM-II phase. 7 In this work we consider only the AFM-II phase. Further investigations concerning the relation between magnetic ordering and electronic band structures will be reported elsewhere.

MnO and NiO
The GW approach has been applied to study the quasiparticle properties of MnO and NiO by several authors. 10,12,13,[15][16][17]21,26 The earliest GW calculations for NiO were carried out with a LMTO-ASA implementation. 10 It was found that the G 0 W 0 @ LDA corrections open the band gap from 0.2 to 1.0 eV, which is still dramatically underestimated compared to the experimental gap of ϳ4 eV. By further applying a nonlocal potential to e g orbitals and updating G self-consistently, a gap of 5.5 eV was eventually obtained. Massidda et al. employed a self-consistent yet approximate model-GW scheme to study the quasiparticle band structure of MnO ͑Ref. 12͒ and NiO ͑Ref. 13͒ and obtained band gaps in good agreement with experiment. Furthermore, they observed a significant enhancement of the O 2p character in the highest valence-band states, in accord with the charge-transfer model of later transition-metal oxides. 117 These findings were confirmed by more recent GW studies. [15][16][17] Since many aspects of the electronic structure of MnO and NiO have already been extensively discussed in previous studies, 10,13,[15][16][17] we will mainly focus on the effects of U in LDA+ U and GW @ LDA+ U in this work. Figure 4 shows the band gaps of MnO and NiO in LDA+ U, G 0 W 0 , and GW 0 as a function of U. Unlike in ScN and ZnS, the LDA+ U band gaps in MnO and NiO increase nonlinearly with increasing U, as changes in the p-d hybridization counteract the scissor behavior of ␦V . The G 0 W 0 and GW 0 band gaps follow accordingly and exhibit a much stronger U dependence than in ZnS and ScN or the lanthanide sesquioxides. 29 For NiO the U dependence is even more pronounced than in LDA+ U, which saturates at large U. This, however, is an artifact of the LDA+ U approach resulting from the fact that at around U = 6 eV the d bands crossover with the Ni 4s derived states. Past this point the band gap changes little because the description of the Ni 4s states remains essentially at the LDA level. In GW @ LDA+ U, on the other hand, all bands are subject to the GW corrections, including Ni 4s derived states. Furthermore, we observe that the G 0 W 0 and GW 0 corrections to the LDA+ U band gap are always positive for MnO and NiO, unlike in ScN or, as we will demonstrate later, FeO. With regard to the order of G 0 W 0 and GW 0 MnO behaves more like a conventional sp semiconductor where GW 0 gives larger band gaps than G 0 W 0 . 63 In Fig. 5 we illustrate the evolution of the density of states with U in LDA+ U and GW 0 for MnO and NiO. Due to their similarity with GW 0 the G 0 W 0 DOSs are not shown here. The most distinct feature we observe apart from the continuous band gap opening with increasing U is a narrowing of the valence bandwidth. In MnO, the LDA+ U and GW 0 @ LDA+ U valence bandwidth decreases by ϳ1 eV over a U range from 0 to 8.2 eV. In NiO, the LDA+ U valence bandwidth first decreases and reaches a minimum at U ϳ 6 eV before it increases again, whereas GW 0 @ LDA + U exhibits a steady decrease by nearly 1.4 eV as U changes from 0 to 8.2 eV. In contrast to ScN and ZnS the transitionmetal d electrons hybridize strongly with the oxygen 2p states and their relative energetic position determines the degree of hybridization and the width of the valence band. As the occupied 3d states move toward lower energy the top of the valence band develops more O 2p character with increasing U. This was also observed in previous GW calculations for MnO and NiO. 10,12,13,15,16 To better understand the U dependence in later transitionmetal oxides, we apply the same analysis to NiO that we have previously applied to ScN. Figure 6͑a͒ shows the band  gap at the ⌫ point determined with Eq. ͑42͒ and its two components as a function of U. The second term increases again as the increasing LDA+ U gap decreases the screening strength. In contrast to ScN, however, the first term also grows with increasing U. This behavior can be explained by the different way that the U-correction term influences the character of the highest valence-and lowest conduction-band states. As already pointed out previously, the U correction significantly enhances the O 2p character in the highest valence-band states so that the latter become more delocalized as U increases. This is reflected in the quasiparticle renormalization factor Z nk , which increases significantly from 0.57 to 0.82 in the investigated U range ͑conventional sp semiconductors have renormalization factors around 0.8͒. 118 As a result, the first term in Eq. ͑42͒ for the valenceband maximum ͑VBM͒ state moves to lower energy as U increases. The lowest conduction-band states, on the other hand, are dominated by Ni 3d states at U = 0 but become itinerant band states ͑of mainly Ni 4s character͒ at finite U. This is again clearly seen in the Z nk factor, which levels out at a value of ϳ0.8 for U's larger than 3 eV. In that case, the term ⑀ nk LDA corresponding to the conduction-band maximum state becomes constant for U տ 3 eV. The combination of these two factors explains why the contribution of the first term in Eq. ͑42͒ increases as a function of U. We observe that the Hubbard U correction has a much stronger influence on the hybridization between localized and itinerant states in these open d-shell transition-metal oxides than in ScN, ZnS, or the 4f lanthanide oxides. 29 This gives rise to a much more pronounced and nonlinear U dependence in both LDA+ U and GW @ LDA+ U. This behavior is consistent with the different extent to which localized states participate in the chemical bonding in these compounds.

FeO and CoO
Compared to MnO and NiO, FeO and CoO are less well studied theoretically. They also pose a more serious challenge for first-principles approaches, and the first GW study for FeO and CoO has been reported only recently. 26 In a cubic crystal-field environment without symmetry breaking, the ground-state high-spin configurations of Fe 2+ ͑t 2g ↑3 e g ↑2 t 2g ↓1 ͒ and Co 2+ ͑t 2g ↑3 e g ↑2 t 2g ↓2 ͒ are degenerate, and therefore any band theory would predict FeO and CoO to be metallic. In the AFM-II phase a symmetry reduction that breaks the degeneracy is possible, but LDA and GGA, due to their strong self-interaction error, fail to induce such a break in symmetry.

Comparisons with experiment
Finally we compare our GW @ LDA+ U results for later transition-metal oxides, using U and J determined by con-strained DFT calculations, to experiment. Table III summarizes the fundamental band gaps of MnO, FeO, CoO, and NiO obtained in this work in relation to experimental and other theoretical results. Compared to widely cited experimental values, the band gaps from GW @ LDA+ U are considerably underestimated. We note, however, that a direct comparison of the fundamental band gap between theory and experiment has to be done with caution. An accurate experimental determination of fundamental band gaps is anything but trivial and the accuracy can be impaired by several factors, including sample quality, limited instrumental resolution, suitability of the chosen technique, or a mix of bulk and surface features. The conceptually most straightforward way to obtain band gaps is by direct and inverse PES. 126 While direct PES is now a routine technique used for probing electronic properties of occupied ͑including deep core and valence-band͒ states with ever-increasing resolution, inverse PES is not nearly as developed and the resolution is typically limited by charge accumulation on the sample surface in the course of the measurement. In practice, the band gap is therefore often measured with optical-absorption techniques. However, features in the absorption spectrum are mainly determined by dipole selection rules and the strength of the electron-hole attraction ͑excitonic effects͒. 127 In addition, the determination of the band edge can be obscured by  Table I 126 Further studies in this direction and more refined experiments are clearly needed for a more unambiguous comparison between experiment and theory.
In Fig. 9, we compare the theoretical DOS with photoemission spectra ͓x-ray photoemission spectrum ͑XPS͒ for occupied states and bremsstrahlung-isochromat spectroscopy ͑BIS͒ for unoccupied states͔. The theoretical and the experimental spectra are aligned in terms of the upper valenceband edge and/or the main VB peak position and not the Fermi level since this is not well defined at the 0 K at which the calculations are performed. Since the alignment is not unambiguous, caution has to be applied when comparing theory and experiment. In general, we observe that the DOSs from LDA+ U, in spite of its great improvement over LDA, are significantly different from experiment. For GW 0 @ LDA+ U, however, most of the peaks lie at the same position as in XPS+ BIS, with the notable exception of the lowest peak in the BIS spectrum of FeO. In NiO, for example, not only the main peak in the valence-band region is accurately described but also small peaks at −2.5, −4, and −5.5 eV are all well reproduced. Difference in the peak intensities are most likely due to final-state effects in the photoemission experiments that are not taken into account in our theoretical treatment but higher order correlation effects that go beyond the GW approach cannot be ruled out at present. Another feature that is absent in the GW 0 @ LDA+ U spectra of all four compounds is the spectral weight at higher binding energies ͑ϳ−9 eV͒, the so-called satellite structure. Its origin has been attributed to transitions from transition-metal d electrons 129 and it is widely believed that a many-body description beyond the GW approach is required for its description. 10,62 Our GW 0 DOSs do, however, exhibit peaks with strong d character 2 eV higher in energy. This energy difference is of the same order as the underestimation of d-electron binding energies in ZnS and other II-VI compounds and group-III nitrides and one could therefore surmise that the failure of GW has the same physical origin in both cases. This is consistent with recent observations by Shishkin et al. 42 that a correct description of ⑀ d in ZnS requires higher order vertex corrections in both the screened Coulomb interaction W and the self-energy.
Regarding the satellite structure in NiO, we note that the GW 0 DOSs reported in this work are obtained directly from quasiparticle energies. An alternative way to characterize electronic states of these systems is to calculate the spectral function directly from the interacting Green's function with the full energy dependence and nonhermiticity of the selfenergy taken into account. The latter may well exhibit features that could be related to the satellite structure. Some analysis along this line was already conducted by Aryasetiawan and Gunnarsson ͑AG͒, 10 who concluded that the satellite structure in NiO arises from higher order correlation beyond the GW approximation. Considering the approximations employed in AG's work in which the LMTO-ASA approach was used for solving Kohn-Sham and GW equations, and a nonlocal potential applied only to the e g band was introduced to overcome the failure of LDA, we consider this issue as still not fully settled, and will pursue it in future work.
In Fig. 10 we further compare theoretical projected density of states on the oxygen atom to experimental spectra obtained from oxygen K␣ x-ray emission spectroscopy ͑XES͒ and oxygen 1s x-ray absorption spectroscopy  Table I are compared to oxygen x-ray emission and absorption spectra ͑Ref. 123͒.
͑XAS͒. 123 In our calculations we do not explicitly create a core hole, whose effect is expected to be small for oxygen emission and absorption spectra. 123,130 We also do not consider the influences of absorption cross sections. These two facts combined mean that only a comparison of peak positions and not relative intensities is meaningful. Taking this into account, we again see an overall good agreement between our GW results and experiment and all main features are well reproduced, in particular, the peak spacing of the three peaks in the XAS spectra of NiO.

IV. CONCLUSIONS
In this work, we employed many-body perturbation theory in the G 0 W 0 and GW 0 approaches based on LDA + U for a series of prototypical d-electron systems. We provide a derivation of the LDA+ U approach as an approximation to GW. Applied to the empty d-state system ScN, we found that LDA+ U gaps ͑indirect and direct͒ depend strongly on U, and increase linearly with U, but G 0 W 0 and GW 0 band gaps exhibit only a weak U dependence. In ZnS with shallow semicore d states, we observe a weak band-gap dependence on U in both LDA and GW. d-electron binding energies can be tuned to their experimental value in LDA+ U but the GW corrections restore the binding energies to their value in LDA-based G 0 W 0 calculations, which is considerably underestimated. In late transition-metal oxides with partially occupied d states, we note a strong U depen-dence in LDA+ U and GW @ LDA+ U. For U's determined from constrained DFT reasonable agreement between GW 0 @ LDA+ U and direct and inverse photoemission data is observed. The fact that the different systems show a very different U dependence in LDA+ U and GW @ LDA+ U can be understood by analyzing how U changes single-particle wave functions and the screening strength in the screened Coulomb interaction W. For late transition-metal oxides the strong U dependence can then be traced back to the large U-induced change in hybridization between the oxygen 2p states and the transition-metal 3d orbitals. To be consistent the U's to be used in GW @ LDA+ U calculations should be determined self-consistently from constrained RPA calculations 90,98,131 which will be the subject of future work. With U determined in this fashion differences between experiment and GW @ LDA+ U would indicate where manybody effects beyond GW are required.