Reformulation of DFT þ U as a Pseudohybrid Hubbard Density Functional for Accelerated Materials Discovery

Luis A. Agapito, Stefano Curtarolo, and Marco Buongiorno Nardelli Department of Physics, University of North Texas, Denton, Texas 76203, USA Center for Materials Genomics, Duke University, Durham, North Carolina 27708, USA Materials Science, Electrical Engineering, Physics and Chemistry, Duke University, Durham, North Carolina 27708, USA Department of Physics and Department of Chemistry, University of North Texas, Denton, Texas 76203, USA (Received 10 June 2014; published 28 January 2015)


I. INTRODUCTION
High-throughput quantum-mechanics (HTQM) computation of materials properties by ab initio methods has become the foundation for an effective approach to materials design, discovery, and characterization [1].This datadriven approach to materials science currently presents the most promising path to the development of advanced technological materials that could solve or mitigate important social and economic challenges of the 21st century [1][2][3][4][5][6][7][8][9][10][11].In order for this approach to be successful, however, one needs the confluence of three key factors: (i) improved computational methods and tools, (ii) greater computational power, and (iii) heightened awareness of the power of extensive databases in science [3].While the last two are driven by technological advances in computing and research and development needs, the development of improved computational tools appropriate for HTQM frameworks is a grand challenge that is still in need of considerable advances.
Most HTQM infrastructures rely on the predictive capability of density-functional theory (DFT), the method of choice for the first-principles study of materials properties.However, despite the enormous success of DFT in describing many physical properties of real systems, its limitations in correctly describing the electronic band structure of insulators are well known.The method is limited by the presence of an unknown correlation term that represents the difference between the true energy of the many-body system of the electrons and the approximate energy that we can compute.The common approximations based on a nearly homogeneous electron-gas treatment of the electron density, the local-density approximation (LDA), and the generalized-gradient approximation (GGA) are extremely successful in the description of many physical properties of materials but dramatically underestimate the electron energy gap in insulators and semiconductors and thus fail to satisfactorily describe the electronic properties of these systems.Higher-order levels of theory exist that are able to predict, with great accuracy, the energy gaps (the GW approximation [12] and dynamical mean-field theory [13][14][15], among others), but they are computationally expensive and unsuitable for extensive high-throughput materials characterization [16,17], even when machine learning methods are employed to simplify the complexity of the task, especially when seeking for new materials systems [18,19].In order to address the energy-gap problem at a lower computational cost, the two most common corrections to traditional local and nonlocal approximations to DFT are "hybrid functionals" and DFT þ U.Both approaches aim to reduce, at some level, the self-interaction error [20] and introduce the derivative discontinuity in the exchange-correlation functional [21][22][23] responsible for the underestimation of the energy gap [24].
Hybrid functionals are based on the idea of computing the exact exchange energy from the Kohn-Sham wave functions and to mix it with the (semi)local approximation of exchange energy of DFT [25].The method is very successful in predicting the energy gap with good accuracy, and it points to the importance of introducing some degree of exact exchange in traditional DFT for a proper description of the electronic structure.However, the level of mixing is not determined from first principles, so the method suffers from some level of empiricism and, even after the introduction of range-separated functionals suitable for periodic system calculations [26,27], it is more computationally demanding than LDA or GGA.
The treatment of systems with strongly localized (correlated) electrons is another outstanding issue that limits the predictive power of (semi)local approximations to DFT for systems where localization is important, such as transition-metal insulators, and it directly impacts the energy-gap problem.The DFT þ U method introduced by Liechtenstein and Anisimov [28,29] aims at preserving the information of orbital localization from being averaged out as in LDA or GGA in order to improve the description of the band structure with a very modest increase in computational effort.Since, in this paper, we are mainly concerned with an extension of this approach, let us discuss its guiding principles in some detail.
Within the DFT þ U ansatz, localized states φ i largely retain their atomic nature and, therefore, can be expanded in term of an atomic-orbital basis set fϕ m g ≡ fmg.The Coulomb and exchange energy associated with these states is explicitly evaluated using the Hartree-Fock (HF) framework via electron repulsion integrals (ERI, also know as two-electron integrals), with a screened (renormalized) Coulomb interaction V ee , and added to the original DFT energy after the removal of double-counting terms in the energy expansion, in a spirit similar to the hybrid functionals approach.
The HF Coulomb and exchange energy of the localized states is given by [28] þ ðhmm 00 jV ee jm 0 m 000 i − hmm 00 jV ee jm 000 m 0 iÞ where n σ is the spin-density matrix n σ of the atomic orbitals ϕ m .This equation can be simplified via the introduction of the phenomenological parameters Ū and J that describe the on-site Hubbard-like interactions as expressed by Dudarev et al. [30]: Here, N σ m is the spin occupation number of the atomic orbital ϕ m .
From the equations above, it clearly follows that the new parameters, Ū and J, contain the information of all the ERIs in an averaged scenario.In physical terms, Ū is the strong correlation experienced between localized electronsonly subtly coupled to the sea of extended states in which they live.Thus, the most akin definition of Ū (for the nonspin-polarized case) is the average [31] where 2l þ 1 is the total number of localized states φ i , and l ¼ 2; 3 for d; f orbitals, respectively.The exchange contribution J is given by a similar average [31].
Although the physical picture is clear, an unambiguous procedure for computing f Ū; Jg from ab initio does not exist.Two factors need to be further clarified in Eq. ( 3): (i) the screened (renormalized) Coulomb interaction V ee arising from the "subtle coupling" to the background extended states and (ii) the actual orbitals φ i used to represent the "localized electrons." Amongst the most common ab initio methods to compute Ū are the constrained random-phase approximation (cRPA) [32] and the linear-response constrained DFT (or cLDA) [33,34].The former computes the screened Coulomb interaction as the bare Coulomb interaction renormalized by the inverse dielectric function, which is calculated using the random phase approximation.The latter circumvents the ambiguity of V ee by indirectly determining Ū as the second derivative of the total energy with respect to constrained variations of the atomic charge q I of the chosen Hubbard center I, Ū ¼ ∂ 2 E=∂q I 2 .E is the total energy of a supercell large enough to converge to the bulk environment for the atom I.It is assumed that the charge perturbation on atom I does not disturb the local environment.In DFT with a linear-combination-of-atomic-orbital (LCAO) basis, this is enforced by suppressing the hopping integrals to prevent charge rehybridization or transfer with its environment [33] and, in the case of plane-wave DFT, by subtracting a correcting term from ∂ 2 E=∂q I 2 as given in Ref. [34].This method has been widely used for open-shell systems; nonetheless, the numerical reliability becomes challenging for closed-shell systems where the localized bands are completely full, thus exhibiting very small response to the linear perturbation [35,36].
Regarding the representation of the φ i states, e.g., d or f electrons of transition metals, localized orbitals-obtained either from the linear-muffin-tin-orbital (LMTO) method [37] or from the Nth-order muffin-tin-orbital (NMTO) method [38]-can be used with both the cLDA and cRPA to obtain Ū.Recently, maximally localized Wannier functions (MLWF), an invariant choice suitable for plane-wave calculations, have also been employed [31,39].By construction, these functions are associated with a given angular momentum (l; m), and the direct correspondence makes them convenient to represent the localized d or f electrons.Ultimately, pinpointing a single localized state within the solid is arbitrary-any of these options are equally valid.In principle, the options should be equivalent for very localized states; nonetheless, the physical significance and construction become more ambiguous when bands corresponding to localized states are not fully disentangled [39].Despite attempting the computation of the same physical entity, the cLDA and cRPA methods do not yield the same value of Ū [40].Considering the number of assumptions taken in numerical implementations, the outcome is unsurprising.
In this article, we introduce an alternative ab initio method to compute Ū and J, which parallels the calculation of the HF energy for molecules and solids and follows closely the original definition of Anisimov et al. [Eq.(1)]: the Agapito-Curtarolo-Buongiorno Nardelli (ACBN0) pseudohybrid Hubbard density functional.In ACBN0, the Hubbard energy of DFT þ U is calculated via the direct evaluation of the local Coulomb ( Ū) and exchange ( J) integrals in which the screening of the bare Coulomb interaction is replaced by a renormalization of the density matrix.Through this procedure, the values of Ū and J are thus functionals of the electron density and depend directly on the chemical environment and crystalline field, introducing an effective procedure of giving the proper description of Mott insulators and other strongly correlated transition-metal oxides (TMOs).As a first application, we discuss the electronic properties of a series of transition-metal oxides that show good agreement with hybrid functionals, the GW approximation, and experimental results at a fraction of the computational cost.In particular, we demonstrate that the ACBN0 functional satisfies the rather ambitious criteria outlined by Pickett et al. in one of the first seminal articles on LDA þ U [41]: (i) ACBN0 reduces to (LDA)PBE when (LDA)PBE is known to be good, (ii) the energy is given as a functional of the density, (iii) the method specifies how to obtain the local orbital in question, (iv) the definition of Ū and J is provided unambiguously, and (v) the method predicts antiferromagnetic insulators when appropriate.
The article is organized as follows: The methodology is discussed in Sec.II.The application of the method for four prototypical transition-metal oxides is presented in Sec.III, and the results are compared against available experimental and theoretical data.Section IV discusses the important features of the method and suggests extensions of potential significance to the goal of discovering novel functional materials.Conclusions are summarized in Sec.V.

II. METHODOLOGY
The foundations of the approach for evaluating the on-site Coulomb and exchange parameters are as follows.
(i) Ū and J are on-site quantities derived from the energy in the Hartree-Fock method.The HF theory considers pairwise interactions of only two electrons at the time, and therefore, it misses the concept of screening.Because of this approximation, the HF theory is known to be inappropriate in describing delocalized metallic systems; however, it is qualitatively sound for molecules and insulators (especially in the strong localization regime) [42][43][44][45].The on-site energies are derived from the HF energies following the ansatz of Mosey and Carter [46,47] in which the occupied molecular orbitals (MO), which are needed for the computation of the HF energy, are considered populated only in the subspace fmg.(ii) No localized orbitals φ i need to be explicitly computed.As in the Hartree-Fock method, all the MOs, or crystalline wave functions for the case of solids, are used.This eliminates the indeterminacy in finding the subset of MOs that better corresponds to the localized states, which can lead to wide fluctuations of the calculated Ū [48].During the calculation of the on-site HF energies, the localized orbitals are implicitly taken as a linear combination of the basis functions of interest, fmg, with the expansion coefficients included in the renormalized density matrix coming directly from the solution of the Kohn-Sham equations projected onto the localized basis of choice (see below).(iii) A plane-wave basis set is the natural choice for DFT calculations of periodic systems, but on-site HF energies are more efficiently computed in a localized basis set.Electron-repulsion integrals are evaluated using pseudo-atomic orbitals (PAO) expressed as a linear combination of Gaussian-type functions, which we define as the PAO-3G minimal basis set.This is possible by the projection procedure that we have recently developed [11], which seamlessly maps the plane-wave electronic structure onto a localized atomic-orbital basis set (see Appendix A).However, it is important to note that the construction of Ū and J outlined below is completely general and can be applied to any choice of basis, localized or otherwise.(iv) E fmg HF is a true functional of the electron density in the spirit of the Hohenberg-Kohn theorems.This leads to the definition of the ACBN0 pseudohybrid Hubbard density functional.

A. Calculation of the electron repulsion integrals
The enormous quantity of ERIs needed in the calculation of the HF exchange energy is the fundamental bottleneck in the use of hybrid DFT functionals.In DFT calculations based on LCAO (PAO) basis sets, the problem is made more tractable when the PAOs are expressed as linear combinations of Gaussian-type functions, as is commonly done in commercial packages such as Gaussian09 [49] and Crystal06 [50].
The electron repulsion integrals used in Eq. ( 1) are defined as four PAOs interacting under the bare Coulomb interaction V ¼ jr 1 − r 2 j −1 as ERI≡ðmm 0 jm 00 m 000 Þ ≡ hmm 00 jVjm 0 m 000 i The real-space evaluation of these integrals is not directly possible when using a plane-wave basis set, which is the preferred choice for periodic systems.For this reason, we employ the auxiliary space of PAOs naturally included in the definition of the pseudopotentials.Given that the radial and angular parts of the PAO basis functions, ϕ lm ðrÞ ≡ ðR l ðrÞ=rÞY mfc;sg l ðθ; φÞ, are separable, they can be directly fitted using linear combinations of sphericalharmonic Gaussian functions.For efficiency, the latter functions are then further expanded as a linear combination of Cartesian Gaussians defining the PAO-3G minimal basis set.(See Appendix A for more technical details on these transformations).Once expressed in the PAO-3G basis set, the ERIs can be efficiently evaluated using any optimized quantum-chemistry library.We use the C routines included in the open-source quantum-chemistry package PyQuante [51].

B. Hartree-Fock Coulomb and exchange energies
The knowledge of the ERIs and the molecular (or crystal) orbitals allows the calculation of the HF Coulomb and exchange energies E HF .For isolated systems (molecules or clusters) and in the restricted case, Here, ψ σ i ðrÞ ¼ P i;μ c σ μi ϕ μ ðrÞ are occupied molecular orbitals expanded in the PAO basis; N σ ψ i ≡ i , and S μν is the overlap integral between the PAOs ϕ μ and ϕ ν .The last line of Eq. ( 5) is expressed in the basis of atomic orbitals ϕ μ with the density matrix P σ μν ¼ The expression of the Coulomb and exchange HF energies for a periodic system is analogous to the molecular case (see Pisani et al. [52]): where g, m, and l are lattice vectors and 0 refers to the primitive unit cell.However, the mapping of the crystalline wave functions ψ kσ i in a local basis (i.e., the expansion coefficients c kσ μi ) is not readily available when using a plane-wave basis to solve for the electronic structure of the material as is common for solids.We circumvent this problem by projecting the plane-wave solution into the chosen auxiliary space of PAOs following the method described in Ref. [11].This projection procedure is a noniterative scheme to represent the electronic ground state of a periodic system using an atomic-orbital basis, up to a predictable number of electronic states, and with controllable accuracy by filtering out high-kinetic-energy planewave components.See Appendix B for a summary of this procedure to calculate the expansion coefficients c kσ μi and the real-space density matrices of the solid, P σ;R μν .
C. Ū and J as functionals of the density: The ACBN0 functional The energy functional for the DFT þ U method is given by where E DFT is the DFT energy calculated using a LDA or GGA functional.The energy correction E U is given either in the original Anisimov-Liechtenstein formulation [28,29] or in the simplified Dudarev [30] formulation as with E fmg;I HF defined in Eq. (1) for a given atom I. E DC corrects for a possible double counting of the localizedstate interaction energy already captured (in an averaged way) in E DFT .The second formulation defines an effective on-site Coulomb interaction U eff ¼ Ū − J (henceforth referred to simply as U).It should be noticed that numerical implementations of the Anisimov DFT þ U functional [Eq.(1)], for instance, in QUANTUM ESPRESSO [53] or VASP [54], do not compute the ERIs explicitly.They are evaluated from tabulated Slater integrals, which ultimately depend on the provided values of Ū and J, or from phenomenological considerations (e.g., Refs.[9,55]).
On the contrary, we evaluate Ū and J by directly computing the on-site Coulomb and exchange energies on the chosen Hubbard center, from the Coulomb and exchange Hartree-Fock energies of the solid.The following assumptions are used.
(i) We follow a central ansatz, introduced by Mosey et al. [46,47] for the case of cluster calculations, that defines a "renormalized" occupation number which is the Mulliken charge of the basis f mg.The set f mg includes all the atomic orbitals in the unit cell that have the same quantum numbers as the orbitals fmg of the Hubbard center of interest.Correspondingly, we define a renormalized density matrix as The renormalized occupations can be interpreted as weighting factors that specify the on-site occupation of each electronic state.The expressions in Eqs. ( 8) and ( 9) are applicable to isolated systems (molecules and clusters).Solidstate systems are periodic and calculated in reciprocal k space that inherently applies periodic-boundary conditions; thus, surface effects are avoided.Periodic plane waves are the basis of choice for the solution of solid-state systems.The k-space electronic-structure information of the material can be projected into a PAO basis directly from the plane-wave DFT solution.We take advantage of this to compute the reduced density matrices and occupation number, which are needed to calculate Ū and J, as follows Here, N k is the total number of k vectors in the first Brillouin zone.See Appendix B for more details.
(ii) E fmg HF , the on-site HF energy associated with the basis fmg is obtained from Eq. ( 5) by restricting the summation indices to fmg.In the periodic case, it is reduced from Eq. ( 6) considering the central unit cell only, i.e., lattice vectors R ¼ g ¼ l ¼ m ¼ 0. Combining (i) and (ii), we obtain, in the general spin-unrestricted case, In this notation, each primed and unprimed index m runs over the whole fmg set.Clearly, the above equation is equivalent to Anisimov's original DFT þ U functional [Eq.(1)] after we replace n σ mm 0 with the renormalized density matrix Pσ mm 0 .σ ¼ fα; βg.However, while Eq.(1) requires the knowledge of a subjective screened Coulomb interaction V ee , Eq. ( 11) uses the bare Coulomb interaction.Although screening is not considered in the HF theory, the direct parallel between both equations shows that the renormalization of the density matrix effectively introduces a degree of screening.
The comparison of Eq. ( 2) with Eq. ( 11) leads to the definitions of Ū and J as density-dependent quantities in the ACBN0 functional: There are essential and relevant differences between Eqs. ( 12) and ( 13) and the similar framework of Mosey et al. [47]: (i) Our on-site energy requires the computation of a smaller number of electron-repulsion integrals, namely, only those involving the fmg set (i.e., 5 4 integrals for the d shell).This directly parallels the original definition of Anisimov, in sharp contrast with the methodology of Mosey et al. requiring ERIs between all basis functions contained in the cluster; (ii) we consider the larger set f mg of localized orbitals (instead of fmg) in the calculation of the reduced occupation number given in Eq. ( 8); (iii) in the approach of Mosey et al. [47], the atom of interest is embedded in a cluster of volume large enough to yield local convergence to bulk conditions, thus requiring calculations on clusters of increasing volume to ensure convergence, which can become extremely computationally expensive; (iv) finally, we do not solve the full Hartree-Fock problem for the solid (Roothaan's equations) but project the DFT Kohn-Sham wave functions on the minimal PAO-3G basis set, thus implicitly including, in the renormalized density matrix, all the local screening effects that come from the mean-field solution on the local set.

III. RESULTS
We selected four prototypical examples to benchmark the ACBN0 density functional: TiO 2 (rutile), MnO, NiO, and ZnO (wurtzite) are technologically important TMOs that have been extensively studied both theoretically and experimentally.These materials pose a methodological challenge to traditional energy functionals (LDA or GGA) because of the strong localization of the TM-3d electrons that results in significant errors in the description of their electronic structure [56].The set of chosen TMOs covers a wide range of 3d shell fillings, namely, 3d 2 , 3d 5 , 3d 8 , and 3d 10 in Ti, Mn, Ni, and Zn, respectively.
All our calculations use the Perdew-Burke-Ernzerhoff (PBE) [57] functional as a starting point and a planewave energy cutoff of 350 Ry with a dense Monkhorst-Pack mesh to ensure good convergence of all quantities.All DFT þ U calculations use the simplified rotationalinvariant scheme of Dudarev [30] and Cococcioni [34] as implemented in the QUANTUM ESPRESSO package [53].Notice that both the original DFT þ U formulation of Anisimov and the simplified rotationally invariant scheme are equivalent within numerical errors for the same values of Ū and J.For all elements, we used a scalar-relativistic norm-conserving pseudopotential from the PSlibrary 1.0.0 [58].
Although the calculation of the effective values of Ū and J should be performed concurrently within the general Kohn-Sham self-consistent loop for electronic convergence, in this work we follow a simplified scheme: The initial and objective guess for the first DFT þ U calculation is U 2p ¼ 0 eV.Then, the resulting electronic structure is used to compute U ð1Þ for the next DFT þ U step [from Eqs. ( 12) and ( 13)].The process is iterated simultaneously for both transition metal and oxygen atoms until the difference between two subsequent iterations is jU ðnÞ − U ðn−1Þ j < 10 −4 eV.This self-consistent scheme ensures the internal consistency of the results, while the true variational solution using the ACBN0 functional will be implemented in the near future.The converged effective values of U for the transition-metal oxides under study are reported in Table I.All the band structures presented follow the AFLOW standard integration paths [55].
In what follows, the ACBN0 results will be benchmarked against experimental measurements and a higher level of theory, whenever possible.Clearly, meaningful comparisons between theory and photoemission spectra require GW quasiparticle energies (and more for excitonic effects), which go beyond DFT.Even if the DFT þ U formalism has been shown to be a first-order approximation to the GW method for localized states in the static limit [29,59], in practice ACBN0 is clearly not nearly as inclusive and accurate as GW, and a comparison between the results of the two approaches provides us with the most stringent test of the validity of our method.
The valence manifold is predominantly of O-2p character with small Ti-3d hybridization except at the top of the manifold at Γ, where it takes almost exclusively an O-2p character.Conversely, the unoccupied manifold is predominantly of Ti-3d character; the conduction-band minimum (CBM) is at Γ, but in practice, it is degenerate with the minima at R and M. Two regions are distinguishable in the 3d projected density of states (PDOS) in the unoccupied manifold (Fig. 1), and they correspond to the octahedraltype crystal-field splitting of e g (higher energy) and t 2g states (lower energy).
The use of an increasing on-site Coulomb potential U 3d on Ti alone (without correcting the oxygen) has been shown to monotonically open the gap, which reaches satisfactory accordance with the experimental value of 3.03 eV [61] only at values of U 3d ∼ 10 eV [62,63].However, Park et al. [64] have found that large values of the Ti on-site Coulomb interaction (> 7) introduce unphysical defect states in the study of vacancies, and they suggested that a concomitant use of U 2p ¼ 7 eV on oxygen is necessary to achieve both the experimental band gap and a good treatment of vacancy states.
With small Löwdin charges of 0.29e-0.45e(out of 2e) per orbital [65], the Ti-3d states cannot be considered localized, and therefore, the use of large values of U 3d is  II), which improves the DFT prediction by 0.9 eV.Contrary to the predominant focus on the Ti-3d states, our results show that a correction on the oxygen 2p states can alone yield an equally satisfactory band gap.More generally, this suggests that a correct treatment of oxygen 2p states may be more relevant to the correct modeling of TiO 2 vacancies.
The occupied O-2p PDOS [red line in Fig. 1(a)] shows a split into two regions, upper 0-3 eV and lower 3.5-6 eV, which originates in the oxygen 2p crystal-field splitting.The 2p x and 2p y states form sp 2 -like σ bonds contained in the planar Y-shaped OTi 3 subunits, whereas the 2p z states remain as lone pairs perpendicular to the Y-shaped planes.The higher-energy 2p z states correspond to the upper PDOS region.These nonbonding lone pairs have been explained with a simple empirical molecular-orbital model [91,92], whereby the octahedral O h symmetry of the local environment of each Ti coordinated to six oxygen ligands ðTiO 6 Þ 8− frustrates the hybridization of the highest occupied orbitals [93].
Applying the on-site Coulomb potential U 2p on oxygen increases the localization of the 2p z lone pairs, thus increasing the splitting between the two 2p PDOS regions, cf.Figs.1(a) and 1(b).The main peak of the lower PDOS region downshifts by 1 eV, to around 5.4 eV, consistent with the value of 5 eV reported by the accurate full-frequency-dependent GW calculation of Khan and Hybertsen GW [94] and x-ray photoelectron spectroscopy (XPS) measurements [95].
Examining the G 0 W 0 @GGA band structure reported by Malashevich et al. [73], it is interesting to notice that besides a scissor-shift operation, the main correction with respect to the DFT bands is a downward energy shift of the 2p x 2p y bands (lower region of the occupied manifold, −6 to −4 eV), whereas the upper region (−4 to 0 eV) remains mostly unchanged.A possible mechanism is that the GW approach implicitly applies a self-interaction correction that increases the splitting between the 2p z and 2p x 2p y states  by further localizing the 2p z states, which can be captured in the ACBN0 calculation.In this regard, ACBN0 bands closely follow the G 0 W 0 @GGA bands of Ref. [73] in the range from −4 to 8 eV.The most significant difference happens in the remaining range of −7 to −4 eV, where our DFT þ U bands are overdownshifted with respect to the G 0 W 0 @GGA results.This can be expected from the explicit emphasis of the on-site Coulomb interaction on oxygen in the ACBN0 approach.
Our converged values for MnO are U 3d ¼ 4.67 eV for Mn and U 2p ¼ 2.68 eV for O.This value is in the range of other ab initio Us reported for Mn (3.6-6.04 eV) [41,59,[98][99][100].Note that because of the different assumptions for the physical quantities (i.e., screening, localized states), ab initio values of U are not expected to be unique.A close empirical value of U 3d ¼ 4.0 eV (albeit with no correction on oxygen) has been reported to reproduce well the experimental energy of formations of several manganese oxides [101,102].
Both PBE and ACBN0 band structures are shown in Fig. 2. The bottom of the conduction manifold is an itinerant sp band with noticeable parabolic dispersion, and it is predominantly of Mn-4s character at Γ.The set of low-dispersion bands located above the CBM are predominantly Mn-3d t 2g states.These bands are more narrowly resolved in the case of NiO.
In the occupied manifold, the PBE results (gray lines) distinctly show dispersionless Mn e g bands (separated from the rest at the top of the manifold) centered at about −0.5 eV and t 2g bands at about −1.5 eV.These bands have minor oxygen hybridization, whereas the bands below them (≲ − 1.6 eV) have a strong O-2p character.With the on-site repulsion correction in the ACBN0 results (black lines), the hybridization between Mn-3d and O-2p states increases.As a result, the 3d bands are pushed down in energy while increasing their dispersion (more noticeably on the t 2g bands).This leads to (i) an increase of the bandwidth of the occupied manifold to about 7.5 eV, in good agreement with results from higher levels of theory (7.6 eV GW@LDA þ U [103] and 8 eV sX-LDA [67]), (ii) an increase of the t 2g bonding-antibonding splitting across the band gap and, indirectly, the resolution of the 4s parabolic band (i.e., the energy difference between the parabolic CBM at Γ and the occupied t 2e bands) to 2.21 eV, which compares favorably with the more rigorous selfconsistent GW [77] and GW@LDA þ U [103] predictions of 2.42 eV and 2.27 eV, and (iii) an increase of the energy difference between the CBM at Γ and the unoccupied e g bands, i.e., the band gap.
The indirect band gap improves from the PBE value of 0.98 eV to 2.31 eV; however, the experimental value of 3.6-4.1 eV (Table II) is still underestimated.On the other hand, the magnetic moment is evaluated to be 4.79 μ B (Table III), which matches the experimental value (Fender et al. [104]).
Franchini and coauthors [102] have performed calculations with a larger value of U 3d ¼ 6.0 eV (without U on oxygen) yielding, however, a band gap and magnetization (2.1 eV, 4.67μ B ) smaller than our results.This evinces the importance of having the on-site Coulomb interaction not only on the TM but also on oxygen.Sakuma et al. [105] have shown that O 2p orbitals are considerably localized (as measured by the spread of their MLWFs) in TMOs; correspondingly, cRPA ab initio calculations find the Coulomb repulsion in oxygen Ū2p ≳ 4 eV [99], independent of the TM.Moreover, multideterminant correlated methods (complete active space self-consistent field, CASSCF) have recently shown that, contrary to conventional wisdom, the valence-band edge in NiO is a localized O 2p state [106].
ACBN0's MnO band gap is closer to the predictions of hybrid functionals (sX-LDA 2.5 eV, HSE03 2.6 eV); nonetheless, the more accurate G 0 W 0 @HSE03 and selfconsistent GW methods yield results of 3.4 and 3.5 eV, much closer to the experimental range.This underperformance of the hybrid functionals suggests that correlation effects may be particularly predominant in the case of MnO and, therefore, beyond the Hartree-exchange correction of hybrid functionals.
For NiO, PBE incorrectly locates the parabolic 4s band above the 3d ones, as seen in Fig. 3.With our values of U, 7.63 and 3.0 eV for Ni and O, the 4s CBM is correctly positioned at the Γ point [107], yielding an indirect band FIG. 2. Band structure (spin up) of manganese oxide.All energies are relative to the valence-band maximum E V .Effective values of U ¼ 4.67 eV for the Mn-3d states and 2.68 eV for the O-2p states are used in the ACBN0 calculation.gap (Z-Γ) of 3.8 eV and a direct gap of 4.29 eV.Considering that the 4s CBM has low spectral weight, the direct gap can well account for the dominant first peak at 4.3 eV observed with bremsstrahlung isochromat spectroscopy (BIS) [83].Our value of Ni U 3d ¼ 7.63 eV is in the same range as the effective value reported by Anisimov et al. of 7.1 eV [108], obtained with the cLDA method, widely used for ab initio calculation of U. Similar to the TiO 2 case, the bands around the bottom region of the NiO-occupied manifold are strongly of O-2p character and are noticeably downshifted by the use of O U 2p with respect to the PBE bands.The manifold bandwidth increases by 1.5 eV to 9 eV as seen in Fig. 3.Most band structures reported in the literature do not find an increase of the bandwidth; nonetheless, the same value of 9 eV is obtained with the self-consistent GW reported by Li et al. [107] who argued that such broader bandwidth accounts well for the presence of strong satellite structures observed experimentally in that range of energy [83,109].Admittedly, these satellites are not captured in other approximations such as the model GW of Massidda et al. [110].As pointed out by Gillen and Robertson [67], a bandwidth of 9 eV in NiO is in agreement with experimental measurements of 8-8.5 eV (x-ray emission spectroscopy XES [79]) and 8.5-9.5 eV (ultraviolet photoemission spectroscopy UPS [111]).
In the strongly ionic ZnO, the bands can be readily identified by their dominant orbital character.The bands in Fig. 4 (black lines) from 0 to −6 eV are mostly of O-2p character.The low-dispersion bands around −9 eV correspond to the Zn-3d states.The conduction bands are predominantly of Zn-4s character.
Within PBE, the 3d bands incorrectly overlap with the 2p manifold introducing spurious hybridizations, as shown in Fig. 4 with gray lines, which in turn leads to a strong underestimation of the band gap.The PBE gap is 0.85 eV, while the experimental gap is 3.3 eV [80].Zinc oxide highlights the underlying failure of LDA or GGA in treating materials with localized electrons and thus constitutes a case study for the application of the DFT þ U method.
Our converged values are Zn U 3d ¼ 12.8 eV and O U 2p ¼ 5.29 eV, and they yield a band gap of 2.91 eV, which compares favorably to the experimental value.The bandwidth of the O-2p manifold shown in Fig. 4 is about 6 eV, in accordance with the angle-resolved photoemission spectroscopy (ARPES) value of around 6.05 eV [78].
Although seemingly high, our parameters agree with values reported by Calzolari et al. [112] (U 3d ¼ 12.0, U 2p ¼ 6.5 eV) and Ma et al. [113] (U 3d ¼ 10, U 2p ¼ 7 eV), both of which were found by fitting to reproduce the experimental band gap and position of the 3d bands.
It is established that the 3d bands downshift monotonically with increasing values of U 3d .As the 3d bands downshift, the p-d repulsion with the O-2p bands is decreased, which in turn lowers the energy of the valence-band maximum (VBM) and, thus, monotonically increases the gap [114].After the 3d bands have been fully disentangled from the 2p manifold, however, the 2p bands are well resolved and remain mostly insensitive to a further increase of U 3d .Consequently, the band gap becomes progressively independent of U 3d , and after the 3d bands are fully disentangled, the application of an onsite Coulomb interaction on oxygen becomes necessary to further reach the experimental band gap [113].
For illustration, Fig. 5 shows a comparison of the band structure with different values of U 3d (12.8 and 9 eV), while which is equivalent to Eq. (7b) [34], and the corresponding Hubbard potential is where 0 ≤ λ Iσ m ≤ 1 is the occupation of the orbital ϕ m .For fully occupied orbitals such as Zn-3d in ZnO (Löwdin charge 9.97e out of 10e), i.e., λ m ≈ 1, the Hubbard energy reduces to E U ≈ 0, and the Hubbard potential becomes a rigid shift V U ≈ −U=2 applied to the localized orbitals.In principle, at this limit, the value of U 3d does not change the energy of the material, and it becomes irrelevant in pinning the position of the 3d bands.
The experimental position of the center of the 3d bands is at −7.5 eV [78,115], measured with respect to the VBM (E V ).Other experimental values have also been reported, −8.81 eV [116], −8.6 eV [117], and −7.8 eV [118].Our value of U 3d underestimates the position of the 3d bands at around −9 eV.As discussed above, for fully disentangled and occupied 3d states, the energy of the system becomes almost independent of U 3d because of a singularity of the DFT þ U energy functional.
Similarly, the GW method and hybrid functionals, while correcting the band gap and fully disentangling the 2p and 3p manifolds, consistently miss the position of the 3d bands by about 1 eV [78].Recently, Lim et al. [78] proposed an assisted GW þ V d approach in which the 3d bands are shifted by an ad hoc on-site potential V d ¼ 1.5 eV in the GW self-energy operator.
Analogously, the cLDA method fails in the case of Zn in ZnO.Because of the full occupancy of the 3d states, they become rather insensitive to small linear perturbations, yielding unreliable numerical values of U [36,119].Lee and Kim [35] have proposed an extension to cLDA method for systems with closed-shell localized electrons.They found U 3d ¼ 5.4 eV for Zn by applying a large perturbation potential and correcting for the excess potential needed to reach the onset of the electron-density response.

IV. DISCUSSION
Indeed, the ACBN0 functional satisfies the rather ambitious criteria outlined by Pickett et al. [41]: (I) ACBN0 reduces to (LDA)PBE when (LDA)PBE is known to be good.-Thereduction of U → 0 arises with the introduction of the "renormalized" density matrix P (instead of the regular density matrix P), which makes U dependent on the degree of localization of the Bloch states.
A toy model with two basis functions m and m 0 reveals the scaling of Eq. ( 12) as Ū ∼ 1 4 N m N m 0 ðmmjm 0 m 0 Þ, in contrast to Ū ∼ ðmmjm 0 m 0 Þ when using the regular density matrix instead.Delocalized Bloch states are assumed to be properly described at the LDA(PBE) level.The more delocalized a state, the lower the charge projected inside the atomic sphere (N m ≈ N m 0 → 0), and thereby U vanishes quadratically.See, for instance, the case of Ti U 3d in TiO 2 , or silicon where ACBN0 yields U 2p ≈ 0 eV.(II) The energy is given as a functional of the density.- The value of U in ACBN0 depends only on the electron density, with a set of fixed PAOs fmg chosen a priori to define U.The ACBN0 functional can be considered a zeroth-order pseudohybrid density functional in the sense that it includes an on-site form of the Hartree-Fock exchange.The results for the test systems studied here follow closely the more established and far computationally more expensive sX-LDA hybrid functional.Moreover, the methodology presented in this work can be immediately generalized to evaluate the nonlocal exchange energy for solids by computing the full set of required two-electron integrals.Thus, one could have a hybrid-functional plane-wave DFT calculation that is as fast as LCAO hybrid functionals while still benefitting from the robust parallel fast-Fouriertransform algorithms and systematic basis-set convergence of plane waves.(III) The method specifies how to obtain the local orbital in question.-ACBN0directly parallels the original orbital-dependent DFT þ U functional of Anisimov that uses atomic orbitals fmg.Conceptually, the localized states φ are linear combinations of fmg with the expansion coefficients obtained selfconsistently.Thus, they reflect the chemical environment of the site; however, the expansion coefficients need not be explicitly known.Even though the information of the coefficients is conceptually included in the renormalized density matrix, they are not individually resolved.Such LCAO expansion is more general and suitable for cases where the localized bands are not readily disentangled from other bands, which happens when the orbitals φ are not thoroughly well localized.(IV) The definition of Ū and J is provided unambiguously.-SeeEqs. ( 12) and ( 13).(V) The method predicts antiferromagnetic insulators when appropriate.-This is demonstrated by the results presented for TiO 2 , MnO, NiO, and ZnO.The flexibility of ACBN0 is that it allows the calculation of Ū and J for any atom in the system of interest, yielding, for instance, non-negligible values for the 2p lone pair of oxygen in transition-metal oxides or for the p states of the anion in transition-metal chalcogenides.Through the inclusion of these terms, ACBN0 corrects both the band gap and the relative position of the different bands, in particular, the ones deriving from the d orbitals of transition-metal atoms.This characteristic of ACBN0 is crucial for improving the agreement with experimental results.Generally, the experimental band gap of TMOs cannot be reached if considering only the TM.Paudel and Lambrecht [120] suggested the simultaneous use of U on both the 3d and 4s orbitals of Zn to reach the experimental gap; however, a large value on Zn U 4s ¼ 43.5 eV was needed.Finally, our results predict the stability of the antiferromagnetic phases of both MnO and NiO.However, a more thorough discussion on the relative stability of different magnetic phases will be the subject of a forthcoming publication [121].

V. CONCLUSIONS
In conclusion, we have introduced ACBN0, a pseudohybrid density functional that incorporates the Hubbard correction of DFT þ U as a natural function of the electron density and chemical environment.The values of Ū and J are functionals of the electron density, and they provide a variational way of obtaining the proper description of insulators such as transition-metal oxides.Although a more extensive validation of this functional is needed, the first results of our tests show improved agreement to higher levels of theory (hybrid functionals or the GW approximation) and to experimental measurements for the electronic properties of TMOs at a fraction of the computational cost.This is an essential requirement for the designefficient algorithms for electronic structure simulations of realistic material systems and massive high-throughput investigations [1].calculations.We follow the procedure by Mathar [122,126] to further convert each spherical-harmonic Gaussian into a linear combination of Cartesian Gaussians.Then, the Cartesian expansion of the PAOs is where f l;m ðx; y; zÞ is given in Table IV.An example of this fitting procedure is shown in Fig. 6 for Zn-3d and O-2p PAO.Having the PAOs expressed as a linear combination of Gaussian-type orbitals in Eq. (A2) is largely advantageous since it allows computation of the ERIs in a straightforward and analytical way.Furthermore, Gaussians allow the filtering out of ERIs with negligible energy contribution [128] further speeding up calculations, as implemented in the Heyd-Scuseria-Ernzerhof HSE03 [26] hybrid functional.

APPENDIX B: CALCULATION OF THE REAL-SPACE HAMILTONIAN AND DENSITY MATRICES FROM THE PLANE-WAVE ELECTRONIC STRUCTURE
The solid is efficiently calculated using plane-wave DFT on a unit cell with periodic boundary conditions.The planewave basis allows a systematic convergence of the basis-set energy error, which is controlled by a single energy-cutoff parameter.Periodic-boundary conditions are implicit to the plane-wave basis, thus avoiding the presence of surface effects intrinsic to molecular cluster calculations.Moreover, plane waves allow the use of robust and scalable Fourier-transform algorithms.We follow the method described in Ref. [11] to project the k-space electronic structure of the solid onto an atomic-orbital space by filtering out high-kinetic-energy plane waves.The resulting reciprocal-space Hamiltonian H σ;k and overlap S k matrices are then Fourier transformed into real space, resulting in The parameters κ and N, defined in Ref. [11], determine the shifting and filtering for the projection procedure.The overlap integral between a basis function ϕ μ located inside the primitive unit cell (lattice vector 0) and the periodic translation of ϕ ν to lattice vector R is the matrix element S 0R μν ¼ hμ 0 jν R i.The real-space density matrix is then computed as where N kσ ψ i ¼ 1 for all occupied states ψ kσ i .The expansion coefficients c kσ μi are the components of the generalized eigenvectors of H σ;k and S k .[96], 4.79 [104] 1.77 [104], 1.90 [96,127]

FIG. 1 .
FIG. 1.Comparison between the band structures and projected density of states of TiO 2 without on-site interactions U ¼ 0 (a) and with the converged effective values of U ¼ 0.15 eV for Ti and 7.34 for O (b).

FIG. 3 .
FIG. 3. Band structure (spin up) of nickel oxide.All energies are relative to the valence-band maximum E V .The effective values of U 3d ¼ 7.63 eV for nickel and U 2p ¼ 3.0 eV for oxygen are used in the ACBN0 calculation.

FIG. 5 .
FIG. 5. Effect of the zinc on-site Coulomb repulsion U 3d on the occupied 3d bands.An ad hoc value of U 3d ¼ 9.0 eV (a) is compared against the converged value of 12.8 eV (b).Both cases use the converged values of U 2d ¼ 5.29 eV for oxygen.

FIG. 6 .
FIG. 6. Fitting the radial component of the PAOs [r −1 R l ðrÞ, yellow line] as a linear combination of primitive Gaussians (black line).Each contracted Gaussian is the sum of the three primitive Gaussians shown in blue.The Zn d (a) and O p PAOs (b) are taken from the PSlibrary 1.0.0.

TABLE I .
Converged values of the effective on-site Coulomb parameter U (in eV) for the transition metal (TM) 3d and the oxygen 2p states.Our converged values for the rutile environment are Ti U 3d ¼ 0.15 eV and O U 2p ¼ 7.34 eV.These values yield a band gap of 2.83 eV close to the experimental range of 2.8-3.8eV (Table

TABLE II .
Minimum direct and indirect energy band gaps (in eV).

TABLE III .
Local magnetic moments (in μ B ) for the antiferromagnetic states of MnO and NiO.

TABLE IV .
Cartesian expansion of f l;m .