Extensive Benchmarking of DFT+U Calculations for Predicting Band Gaps

Accurate computational predictions of band gaps are of practical importance to the modeling and development of semiconductor technologies, such as (opto)electronic devices and photoelectrochemical cells. Among available electronic-structure methods, density-functional theory (DFT) with the Hubbard U correction (DFT+U) is a computationally tractable approach to improve the accuracy of band gap predictions beyond that of DFT calculations based on (semi)local functionals. At variance with DFT approximations, which are not intended to describe optical band gaps and other excited-state properties, DFT+U can be interpreted as an approximate spectral-potential method when U is determined by imposing the piecewise linearity of the total energy with respect to electronic occupations in the Hubbard manifold, thereby providing a (heuristic) justification for using DFT+U to predict band gaps. However, it is still frequent in the literature to determine the Hubbard U parameters semiempirically by tuning their values to reproduce experimental band gaps, which ultimately alters the description of other total-energy characteristics. Here, we present an extensive assessment of DFT+$U$ band gaps computed using self-consistent ab initio U parameters obtained from density-functional perturbation theory to impose the aforementioned piecewise linearity of the total energy. The study is carried out on 20 compounds containing transition-metal or p-block (group III-IV) elements, including oxides, nitrides, sulfides, oxynitrides, and oxysulfides. By comparing DFT+U results obtained using nonorthogonalized and orthogonalized atomic orbitals as Hubbard projectors, we find that the predicted band gaps are extremely sensitive to the type of projector functions and that the orthogonalized projectors give the most accurate band gaps, in satisfactory agreement with experimental data.


Introduction
Density-functional theory (DFT) [1,2] with approximate exchange-correlation (xc) functionals-e.g., local-density approximation (LDA) or generalized-gradient approximation (GGA)-has been remarkably successful in predicting ground-state properties of a large variety of systems, such as crystal structure and thermodynamic stability. However, these DFT calculations have known limitations, including the underestimation of the band gap (on the order of ∼40% in semiconductors and insulators [3]) due to self-interaction errors (SIE) inherent to approximate xc functionals [4,5]. Various xc functionals and methods were developed to alleviate this problem. One practical solution is to use hybrid functionals, which employ a fraction of the exact (Fock) exchange energy to reduce SIE [6][7][8][9][10][11]. The drawback of this method is its very high computational cost (which prevents its wider application to high-throughput computational studies [12]), and difficulties in determining the magnitude of the screening parameter for solids (that controls the amount of the exact exchange), although there has been considerable progress in this direction in recent years (see, e.g., Refs. [13][14][15][16][17][18]). In parallel, a range of computationally efficient meta-GGA functionals have emerged [19], including the strongly constrained and appropriately normed semilocal density functional (SCAN) [20]. The SCAN functional uses a xc potential which depends on Kohn-Sham (KS) wavefunctions via a kinetic energy density; however, SCAN still contains significant SIE (especially when applied to transition-metal compounds [21,22]) and exhibits some potential limitations in describing magnetic systems [23,24]. An alternative approach is provided by Koopmans-compliant (spectral) functionals [25]. The main idea of these is to enforce the piecewise linearity of the total energy with respect to all orbitals occupations; this class of functionals has proven to be effective in improving band gaps in finite and extended systems [26][27][28][29]; however, this method is not straightforward to implement in practice, though recent progress in this direction has been achieved [30]. A very popular extension of DFT is the DFT+U approach [31][32][33], in which the Hubbard U correction acts selectively on a subset of states in the system (typically, of d or f character) by imposing the piecewise linearity in the energy functional as a function occupations of this subset [34]. DFT+U is only marginally more expensive than DFT within the LDA or GGA approximations, while significantly improving various properties of materials such as transition-metal compounds. A variant of DFT+U to take into account inter-site Hubbard V interactions (DFT+U+V [35]) was introduced, showing success in describing materials with strong inter-site electronic hybridizations (i.e. covalent interactions) [35][36][37][38][39][40][41]. Finally, there are various methods beyond DFT, including many-body perturbation theory (MBPT) [42] (e.g., GW approximation [43][44][45][46]) and dynamical mean field theory (DMFT) [47][48][49][50][51], which are widely used for predicting optical properties of (strongly) correlated systems.
In contrast to MBPT and DMFT, the scope of DFT and DFT+U is limited to ground-state properties, including electronic band gaps but excluding optical band gaps [28,52]. The optical band gap of a material is the minimal energy required for absorbing an incident photon (corresponding to a neutral excitation), whereas the electronic band gap is the difference between the occupied and unoccupied states (corresponding to charged excitations) [52]. In most inorganic semiconductors, the electronic band gap and the optical band gap are approximately equal, but for organic semiconductors, the difference between these two energies (the exciton binding energy) may be substantial. It is thus expected that DFT+U would be more accurate at predicting optical band gaps for inorganic semiconductors than for organic semiconductors [53,54]. To improve the precision of the computed electronic band gaps, DFT+U incorporates a corrective U term designed to restore piecewise linearity of the total energy with respect to occupation of orbitals within the Hubbard manifold, which acts to approximate derivative discontinuities [28,35]. Often the U parameter is chosen semiempirically by fitting it to reproduce experimental band gaps [12,55,56], experimental oxidation enthalpies [57,58], or other properties. However, fitting U to reproduce a subset of experimental data is not a predictive approach, especially in cases when no prior data is available. If instead the Hubbard U parameter is determined using ab initio methods, the band gap prediction might not exactly match that of the experimental gap, but this approach typically improves band gap predictions drastically (over semilocal DFT) and provides an more accurate, overall description of the material [12].
In practical terms, the Hubbard U parameter can be calculated using ab initio methods, such as constrained DFT (cDFT) [34,[59][60][61][62][63][64][65][66][67], constrained random phase approximation (cRPA) [68][69][70][71][72][73][74], and Hartree-Fock based methods [75][76][77][78]. In particular, a linear-response formulation of cDFT was introduced in Ref. [34]: its main idea is based on the fact that the Hubbard corrections are meant to remove SIE from approximate energy functionals that manifests itself through a curvature of the total energy as a function of atomic occupations. Recently, this formulation was recast in the framework of density-functional perturbation theory (DFPT) [38,79], allowing the replacement of computationally expensive supercells by primitive cells at the cost of having multiple monochromatic perturbations instead of just one (localized). In all these ab initio methods for computing U, key is the choice of the Hubbard projector functions. There are different types of Hubbard projectors [80]; in particular, we mention here nonorthogonalized atomic orbitals [34,[81][82][83], orthogonalized atomic orbitals [36,37,39,84], nonorthogonalized Wannier functions [85], orthogonalized Wannier functions [86], linearized augmented plane-wave approaches [87], and projector-augmented-wave projector functions [88,89]. It is well known that the values of the computed Hubbard U parameters strongly depend on the type of Hubbard projectors that are used [79,84]. Moreover, the final properties of interest (e.g., band gaps, oxidation reaction energy, etc.) obtained using DFT+U and different projectors are not always the same [90]. Therefore, special attention must be given to the choice of Hubbard projectors when computing band gaps using DFT+U.
Historically, DFT+U was introduced to accurately describe the dand f -type electrons of transition-metal (TM) and rare-earth compounds, as it was thought that the U parameter would correct these strongly correlated, localized states (i.e., dand f -type states). It was later discovered that the Hubbard U parameter corrects for SIE, which is present in not only d and f , but also s and p states [91,92]. While light elements (e.g., O, N, S) contain only s and p states and thus do not have as large SIE as TM, U has been increasingly applied to light elements in recent years to improve the computational prediction of properties, such as the band gap [78,93,94]. Additionally, in the literature, U has been applied to d 10 elements from the p-block (specifically group III-IV) of the periodic table, such as Ga 2+ [35]. While elements, like Ga 2+ , are not TM, they have electronic configurations that match that of a d 10 TM under certain ionizations. Applying U to light elements and/or d 10 group III-IV elements (e.g., Ge, In, Sn, Pb) is still in debate in the literature [35,78,93,94]. This issue will be addressed in this work using DFPT [79], and in particular we will explore the sensitivity of band gaps to the application of the Hubbard U correction to p states of light elements.
In this paper we present a benchmark of DFT+U using ab initio Hubbard U [79] for predicting band gaps in 20 compounds containing transition-metal or p-block (group III-IV) elements, including oxides, nitrides, sulfides, oxynitrides, and oxysulfides. This is done using two types of Hubbard projector functions, namely, nonorthogonalized atomic orbitals (that we will refer to as "atomic" orbitals, for simplicity) which are provided with pseudopotentials and which are orthonormal within each atom, and orthogonalized atomic orbitals which are obtained by orthogonalizing the atomic orbitals from different sites using the Löwdin method (that we will call "ortho-atomic" orbitals) [95]. By inspecting the dominant character of electronic states at the top of the valence bands and at the bottom of the conduction bands we select the most relevant states and apply the Hubbard correction to them and analyze the subsequent changes in the band gaps. Our study reveals that the DFT+U band gaps are very sensitive to the type of Hubbard projectors, and that the most accurate results are obtained when using "ortho-atomic" projectors. Moreover, this work shows that DFT+U with Hubbard parameters determined nonempirically can be used to predict reliable band gaps at lower computational cost than that of hybrids functionals (such as HSE [96]) that are popular for predicting band gap values. This is important for high-throughput studies that require efficient and reliable estimates of the band gaps at moderate computational costs for various technological applications such as photocatalysis [97].
The remainder of the paper is organized as follows. Section 2 briefly recalls the formulation of the DFT+U approach and how Hubbard U is computed using DFPT, as well as describes which materials are studied in this work. Section 3 contains the technical details of our calculations. In Sec. 4 we present the projected density of states for all materials of this work, computed values of Hubbard parameters, and the comparison of the computed band gaps with those measured in the experiments. Finally, Sec. 5 presents our concluding remarks. Appendix A presents our convergence tests when computing U, and Appendix B contains a comparison of the computational costs of HSE and linear-response U calculations.

DFT+U
In this section, we briefly review DFT+U in the simplified, rotationally invariant formulation [33] which is used in this work. For the sake of simplicity, the formalism is presented in the framework of norm-conserving pseudopotentials in the collinear spin-polarized case. The generalization to the ultrasoft pseudopotentials and the project augmented wave method are described in Ref. [38]. In this section we use Hartree atomic units.
DFT+U is based on an additive correction to the approximate DFT energy functional [33]: where E DFT is the approximate DFT energy (constructed, e.g., within the spin-polarized LDA or GGA), while E U contains the additional Hubbard term: where I is the atomic site index, m 1 and m 2 are the magnetic quantum numbers associated with a specific angular momentum, and U I are the on-site Hubbard parameters. The atomic occupation matrices n Iσ m 1 m 2 are based on a projection of the lattice-periodic parts of KS wavefunctions, u σ v,k (r), on the Hubbard manifold: where v and σ represent, respectively, the band and spin labels of the KS wavefunctions, k indicate points in the first Brillouin zone (BZ) and N k is the number of k-points, f σ v,k are the occupations of the KS states, andP I m 2 m 1 k is the projector on the Hubbard manifold: where ϕ I m 1 k (r) is the Bloch sum computed using ϕ I m 1 (r − R I ) which are the localized orbitals centered on the Ith atom at the position R I . The Hubbard manifold can be constructed using different types of projector functions, as was discussed in Sec. 1; here we consider two types of localized functions, "atomic" and "orhto-atomic".
In DFT+U, the values of Hubbard parameters are not known a priori, and often their values are adjusted semiempirically, by matching the value of properties of interest, which is fairly arbitrary. In this work we instead determine Hubbard parameters from a piecewise linearity condition implemented through linear-response theory [34], based on DFPT [38,79]. Within this framework the Hubbard parameters are the elements of an effective interaction matrix computed as the difference between bare and screened inverse susceptibilities [34]: where χ 0 and χ are the susceptibilities which measure the response of atomic occupations to shifts in the potential acting on individual Hubbard manifolds. In particular, χ is defined as where α J is the strength of the perturbation on the Jth site. While χ is evaluated at self-consistency of the linear-response KS calculation, χ 0 (which has a similar definition) is computed before the self-consistent re-adjustment of the Hartree and exchange-correlation potentials [34]. The main goal of the DFPT implementation is to recast the response to such isolated perturbations in supercells as a sum over a regular grid of N q q-points in the Brillouin zone [79] dn Iσ It is important to remark that the monochromatic responses in Eq. (7) are calculated in the primitive unit cell, thus avoiding the use of the computationally expensive supercells. In Eq. (7), the atomic indices have been replaced by atomic (s and s ) and unit cell (l and l ) labels [I ≡ (l, s) and J ≡ (l , s )]; R l and R l are the Bravais lattice vectors, and the grid of q-points is chosen fine enough to make the resulting atomic perturbations effectively decoupled from their periodic replicas [38,79]. Since ∆ s q n sσ mm are the lattice-periodic responses of atomic occupations to a monochromatic perturbation of wavevector q, they can be obtained by solving the DFPT equations independently for every q, which allows us to parallelize and speed up the calculations of Hubbard parameters. To cover a sufficiently representative set of materials, we have investigated binary oxides, binary nitrides, binary sulfides, oxynitrides, and oxysulfides. Metals in these compounds were limited to d 0 (Sc 3+ , Y 3+ , Ti 4+ , Zr 4+ , Hf 4+ , V 5+ , Nb 5+ , Ta 5+ , Mo 6+ , W 6+ ) and d 10 (In 3+ , Ge 4+ , Sn 4+ , Pb 4+ ) cations to tailor the study towards photocatalytic materials [98]. Using the data from the Materials Project [99], we created permutations of O, N, S, and these d 0 and d 10 cations. This resulted in an extensive list of materials which was paired down using criteria of unit cell size and experimental band gap: the number of atoms per unit cell was limited to <20 and an experimental band gap had to be available in the literature for comparison. Table 1 shows the resulting list, containing 8 oxides, 2 nitrides, 6 sulfides, 3 oxynitrides, and 1 oxysulfide. These materials are all d 0 or d 10 closed-shell systems due to the oxidation states of the TM and group III-IV elements, as listed above.

Technical details
All calculations were performed using the plane-wave (PW) pseudopotential method as implemented in the QUANTUM ESPRESSO distribution [100][101][102]. We used GGA for the exchange-correlation functional constructed with the PBE prescription [103], and we have used the ultrasoft pseudopotentials from the GBRV library [104,105].
KS wavefunctions and charge density were expanded in PWs up to a kinetic-energy cutoff of 60 and 480 Ry, respectively. Electronic ground states were computed by sampling the BZ with uniform Γ-centered k-points meshes with a spacing of 0.04 Å −1 .
DFT+U calculations were performed using the simplified rotationally-invariant formulation [33]. As projectors for the Hubbard manifold, we have used "atomic" and "ortho-atomic" orbitals [84]; for the latter we have used the Löwdin orthogonalization method [95]. Hubbard U parameters were computed using the linear-response approach [34] as implemented using DFPT [79] (see Sec. 2.1). We have used the self-consistent procedure for the calculation of U as described in detail in Ref. [38], which consists of cyclic calculations containing structural optimizations and recalculations of Hubbard parameters for each new geometry. We have performed convergence tests with respect to the q-points sampling for 3 materials (TiO 2 , NbNO, Y 2 SO 2 ); the convergence criteria was ∆U < 0.1 eV. The details can be found in the Appendix A. After selecting the q-points sampling, hybrid functional (HSE) calculations were performed for two representative materials (TiO 2 and NbNO) to compare the computational costs of HSE and U calculations. The results of this analysis are detailed in Appendix B.
All materials were tested for possible finite magnetization by performing collinear spin-polarized calculations within standard DFT. In order to allow for fractional occupations we have used the Marzari-Vanderbilt smearing [106] with a broadening parameter of 5 × 10 −3 Ry. In this case we also sampled the BZ using the uniform Γ-centered k-point meshes with a spacing of 0.04 Å −1 . Materials were only tested for ferromagnetic ordering; thus, this study is not exhaustive and more testing would be needed to check for other types of possible magnetic orderings. Based on our computational testing for ferromagnetic ordering, all of the materials in this study showed zero magnetization. Therefore, all materials were modeled as nonmagnetic, which could lead to errors in predicted band gap values for materials that do show magnetic ordering other than ferromagnetic. In fact, some of the materials are known to be antiferromagnetic (V 2 O 5 ) or ferromagnetic (SnO 2 , SnS 2 ) in their defective forms [107][108][109], but the primary focus of this work is to investigate the pristine bulk form.
We have computed the projected density of states (PDOS) using both DFT and DFT+U with the Gaussian smearing and a broadening parameter of 5 × 10 −3 Ry. The PDOS was calculated by summing up states from all equivalent atoms in the unit cell.
The data used to produce the results of this work are available in the Materials Cloud Archive [110].

Determination of the states requiring Hubbard corrections
One can assess the necessity of applying the Hubbard U correction to certain states of a given element by looking at the PDOS of a given material. The PDOS plots show, in particular, which specific states have the highest density of states (DOS) near the valence band maximum (VBM) and conduction band minimum (CBM) in semiconductors and insulators. By applying U corrections to elements with states at the VBM and CBM, the value of the band gap is expected to change with respect to standard DFT predictions. In contrast, applying U values to elements with low DOS or lacking states near band edges is not expected to significantly affect the predicted band gap value. Using PDOS, we can also assess whether a given material is a band insulator, a Mott-Hubbard insulator, or a charge-transfer insulator [111]. In Mott-Hubbard insulators, the VBM and CBM are of the same kind (e.g., of d character), in charge-transfer insulators instead the VBM and CBM are of different kinds (e.g., the VBM is mainly of the p character and the CBM is of the d character), while in band insulators the band gap is not determined by strongly localized electrons of d or f character. Therefore, by determining what is the dominant character of the VBM and CBM in the PDOS, it is possible to identify the insulator type of a given material. Thus, prior to DFT+U calculations, PDOS calculations were performed in order to determine the dominant character of states at the VBM and CBM, and thus to know which states might require the Hubbard U correction. The PDOS for all materials studied here are shown in Fig. 1. It can be seen that all 14 TM-containing materials can be classified as charge-transfer insulators, while all of the remaining materials that contain group III-IV elements (except InN that comes out to be metallic at the DFT level of theory) are band insulators. Also we note that in Fig. 1 those cases that have multiple lines of the same color indicate the PDOS for non-equivalent atoms (e.g., multiple red lines for V 2 O 5 indicate 2p states of non-equivalent oxygen atoms).
All the materials of this study are oxides, nitrides, and sulfides containing d 0 /d 10 cations. The PDOS of d 0 TM-containing materials (TiO 2 , ZrO 2 , HfO 2 , V 2 O 5 , Ta 2 O 5 , WO 3 , TiS 2 , ZrS 2 , HfS 2 , MoS 3 , Sc 2 S 3 , NbNO, TaNO, Y 2 SO 2 ) in Fig. 1 show dominant d states from the TM d 0 cations at the CBM and dominant p states from the light elements (O, N, S) at the VBM. Hence, one might expect that the application of the Hubbard U correction to these states would allow us to obtain the most accurate prediction of the band gap at the DFT+U level of theory. The PDOS of d 10 p-block (group III-IV) containing materials (SnO 2 , PbO 2 , InN, Sn 3 N 4 , SnS 2 , Ge 2 N 2 O) show dominant p states from the light elements (O, N, S) at the VBM and a mixture of p and s states from the group III-IV elements at the CBM, indicating that the U correction can be applied to the p states of the light and group III-IV elements to obtain the most accurate prediction of the band gap [112]. In addition, it might be useful to investigate in these materials whether the application of the U correction only to the d states of the TM elements in d 0 TM-containing materials or only to the p states of the light elements (i.e. at the VBM) in d 10 p-block containing materials would be sufficient to improve band gaps with respect to standard DFT. Moreover, in this latter case we will show that the application of the U correction to the p states at the CBM turns out to be not so important, since this leads only to minor changes in the band gaps at the DFT+U level of theory. A similar test for d 0 TM-containing materials will not be performed because the d states of TM elements appear not only at the CBM but also they appear in the valence region and show significant hybridization with the p states of light elements, which means that the application of U to the d states of TM elements is expected to be important. Finally, we note that InN is metallic at the DFT level of theory (in Fig. 1 the DOS above the Fermi level is very small and is not well visible). In the next section we present the values of Hubbard parameters computed for the states discussed above that appear at the VBM and CBM for all the materials.

Hubbard U parameters from first principles
In this section we present the Hubbard parameters for all systems studied here, computed using the DFPT approach that was briefly discussed in Sec. 2.1. Based on the PDOS computed at the DFT level of theory (see Sec. 4.1), we calculated U for those states that appear at the VBM and CBM, and the results are shown in Table 2. More specifically, for the d 0 TM-containing materials we computed the Hubbard parameters for the partially occupied d states of TM elements (Sc, Ti, V, Y, Zr, Nb, Mo, Hf, Ta, W) that appear at the CBM and also for the p states of light elements (O, N, S) that appear at the VBM. In the case of d 10 p-block (group III-IV) containing materials we computed the Hubbard U parameters for the p states of light elements that appear at the VBM and for the p states of the group III-IV elements that appear at the CBM.
First of all it is important to discuss the dependence of the computed values of the U parameters on the type of Hubbard projectors that are used. In this work we consider two types of projector functions, namely "atomic" and "ortho-atomic". The difference between the two is the inter-site orthogonalization of atomic orbitals that is present in the latter and is absent in the former. As can be seen in Table 2, the values of Hubbard parameters depend strongly on the type of projector functions that are used. In particular, we can see that U for the d states of TM elements is larger when "ortho-atomic" projectors are used, while U for the p states of light elements are much smaller with respect to the case when "atomic" projectors are used. Notably, U for the p states of light elements in some materials are extremely large, and in fact this finding is not surprising and is well known in the literature [113]: linear-response theory for computing Hubbard parameters according to Ref. [34] is not suited for application to fully occupied states (closed-shell systems), which was observed when using "atomic" Hubbard projectors. Indeed, it can be seen from Table 2 that the values of U for the p states of light elements (O, N, S) are unphysically large when computed using "atomic" projectors. As will be seen in Sec. 4.3, this leads to dramatic worsening of the band gaps when computed using DFT+U with these values of U and "atomic" projectors. However, interestingly we find that when "ortho-atomic" Hubbard projectors are used, the values of Hubbard U for the p states of light elements are much smaller and they fall in the range of physically meaningful values (though still quite large). As will be seen in the following section, the band gaps computed using DFT+U and using these U values with the "ortho-atomic" projectors are in good overall agreement with the experimental band gaps. Therefore, this observation highlights the crucial role played by the Hubbard projectors [90] when computing U and using it for predicting various materials' properties, and in particular band gaps.
Next, in the case of d 0 TM-containing materials, the U values computed for the partially occupied d states were applied to either 3d, 4d, or 5d states depending on the TM period in the periodic table. In general, the smaller the principal quantum number the larger the localization of the d orbitals (because they are closer to the nucleus); thus, the value of U is expected to be larger for 3d states than for 4d and 5d states [114]. Indeed, we see this trend for the U values of the TM d 0 elements in Table 2. In the case of the d 10 p-block containing materials with "ortho-atomic" Hubbard projectors we find that the U values of the p states of the group III-IV elements are very small (< 1.5 eV), which is not surprising since these states are almost fully empty. It will be shown in Sec. 4.3 that the application of the U correction to the p states of the group III-IV elements has only a minor effect on the computed band gaps.

Band gaps from DFT+U calculations
In this section we present our findings for the band gaps of all 20 materials computed using standard DFT and using DFT+U with ab initio values of Hubbard U (see Sec. 4.2) using "atomic" and "ortho-atomic" Hubbard projectors. The resulting predicted band gap values are plotted in Fig. 2 (a), for the materials containing TM elements, and in Fig. 2 (b), for the materials containing p-block (group III-IV) elements; this data is also summarized in Tables 3 and 4.
We use the following short-hand notations: (i) in the TM-containing materials we use "U-TM" for DFT+U calculations with the U correction applied only to the d states of TM elements, and "U-All" for DFT+U calculations with the U correction applied both to the d states of the TM elements and to the p states of the light elements (O, N, S); (ii) in the materials containing the p-block (group III-IV) elements we use "U-Light" for DFT+U calculations with the U correction applied only to the p states of the light elements, and "U-All" for DFT+U calculations with the U correction applied both to the p states of the light and group III-IV elements (Ge, In, Sn, Pb).
Upon inspection of the results for the TM-containing materials [see Fig. 2 (a)] we can see large variations in the computed band gaps depending on the type of projector functions that are used (e.g. compare the slope of the lines obtained by fitting the data points). Overall, very good agreement with the experimental band gap values is obtained for DFT+U calculations using "ortho-atomic" projectors with respective U values. In contrast, DFT+U calculations with the "atomic" projectors result in large, unrealistic values (> 8 eV) for the band gaps in some materials (Ta 2 O 5 , Y 2 SO 2 , ZrO 2 , and HfO 2 ) when Hubbard Data was fit to first order polynomials so that data points from DFT and different DFT+U calculations could be more easily compared using the resulting solid trendlines (note that all data points were included in these fits). DFT+U calculations were performed using "atomic" and "ortho-atomic" projectors. On panel (a), the legend also indicates whether U corrections are applied to only d states of the TM elements (U-TM), or to the d states of TM and p states of light elements (U-All). On panel (b), the legend also indicates whether U corrections are applied to p states of light elements (U-Light), or to p states of light and p-block (group III-IV) elements (U-All). Note that in panel (b), the "DFT+U Ortho-Atomic U-Light" trendline (blue) is dashed so that the "DFT+U Ortho-Atomic U-All" trendline (orange) can also be seen.
corrections are applied to TM and light elements (indicated by U-All) and underestimated band gaps when Hubbard corrections are applied to only TM (indicated by U-TM). In the case of "ortho-atomic" projectors, both U-TM and U-All results are in good agreement (on average) with the experimental band gaps, with U-All fitting line (mean absolute error of 0.55) being better than the U-TM one (mean absolute error of 0.66). This latter finding can be further clarified by looking at Table 3: we can see that when the "ortho-atomic" projectors are used, both U-TM and U-All corrections perform equally well (which one is better depends on the material). It is worth highlighting the outliers from these trends: V 2 O 5 shows comparable accuracy between band gaps computed using "ortho-atomic" projectors with Hubbard corrrections applied to only TM and standard DFT, and TiS 2 shows higher accuracy for the band gap with DFT+U calculations using "atomic" projectors. For V 2 O 5 , it is important to remark that the antiferromagnetic ordering seen Table 3. Comparison of the band gaps (in eV) as computed using DFT and DFT+U (using "ortho-atomic" and "atomic" projectors) and as measured in experiments for compounds containing TM elements, with the mean absolute error (MAE) given for each method. DFT+U results obtained with "otho-atomic" and "atomic" projectors are split into two columns each: (1) band gaps predicted with U corrections applied only to the d states of TM elements (U-TM), and (2) band gaps predicted with U corrections applied to the d states TM and p states of light elements (U-All). The computationally predicted band gaps with the closest values to the experimental band gaps are highlighted in bold.

Formula
Expt. DFT DFT+U (ortho-atomic) DFT+U (atomic) U-TM U-All U-TM U-All Materials containing p-block (group III-IV) elements show even more striking differences between the band gaps computed using "atomic" and "ortho-atomic" projectors with their corresponding values of U parameters [see Fig. 2 (b)]. Notably, it can be seen that the DFT+U band gaps computed using "atomic" projectors are much larger than the experimental gaps, and, surprisingly, they are even worse than DFT band gaps. However, DFT+U band gaps computed using "ortho-atomic" projectors with respective U values are in good overall agreement with the experimental band gaps, and importantly these are better than the DFT band gaps. In addition we find that in the majority of cases U-All corrections give slightly better band gaps than U-Light, though these differences are very small. This latter finding suggests that the application of the U correction to the p states of the p-block (group III-IV) elements is not very important. Thus, U values were not applied to the p states of group III-IV elements for the case of DFT+U band gaps computed using "atomic" projectors.
Therefore, from our analysis of the band gaps computed using DFT+U for all of the 20 materials considered here, we can conclude that the Hubbard projectors play a fundamental role in the accuracy of band gap prediction. More specifically, depending on the type of projector functions, the Hubbard U parameters vary dramatically (see Sec. 4.2) and as a consequence the DFT+U band gaps also vary largely. Therefore, the choice of projector functions must be done with care prior to DFT+U calculations. In particular, in our study we find that the "ortho-atomic" Hubbard projectors give band gaps that are in much better agreement with the experimental values, compared to DFT+U with the "atomic" projectors. This result should be taken into account in future DFT+U studies. In addition, it is important to mention that from Tables 3 and 4 we can see that still our best DFT+U band gaps are not perfect and differ from the experimental ones by as much as almost ∼40% in some cases (e.g., SnO 2 ). This finding suggests that DFT+U with "ortho-atomic" projectors is a valuable tool for improving band gaps over standard DFT and Table 4. Comparison of the band gaps (in eV) as computed using DFT and DFT+U (using "ortho-atomic" and "atomic" projectors) and as measured in experiments for compounds containing group III-IV elements (Ge, In, Sn, Pb), with the mean absolute error (MAE) given for each method. DFT+U calculations using "ortho-atomic" projectors are split into two columns: (1) band gaps predicted with U corrections applied only to the p states of light elements O, N, and S (U-Light), and (2) band gaps predicted with U corrections applied to the p states of group III-IV and light elements (U-All). For "atomic" calculations, band gaps were predicted with U corrections applied only to the p states of light elements (U-Light). The computationally predicted band gaps with the closest values to the experimental band gaps are highlighted in bold.

Formula
Expt. thus can be useful for high-throughput screening of thousands of materials, but whenever one is interested in very accurate predictions of band gaps then it is important to resort to more advances and more accurate methods (see Sec. 1). Finally, for the sake of completeness, we present the PDOS for all 20 materials using DFT+U with "ortho-atomic" projectors in Fig. 3. The U correction was applied both to the d states of TM elements and to the p states of the light and p-block (group III-IV) elements (i.e., U-All). By comparing Fig. 3 with Fig. 1 we can see not only an increase in the value of the band gap thanks to the U correction but also "sharpening" of the PDOS due to the higher level of localization of states that are corrected. Although a comparative analysis of the PDOS on the basis of experimental spectroscopic data is beyond the scope of the present work, it would be an instructive topic for future studies.

Conclusions
We have presented a detailed investigation of the accuracy of DFT+U calculations for 20 materials containing transition-metal or p-block (group III-IV) elements, including oxides, nitrides, sulfides, oxynitrides, and oxysulfides. Hubbard U parameters were computed from first principles using density-functional perturbation theory [37,38], thus avoiding any fitting or tuning parameters. We have found that the accuracy of DFT+U band gaps depends strongly on the type of Hubbard projector functions used for applying the U corrections. More specifically, by comparing DFT+U results obtained using nonorthogonalized and orthogonalized atomic orbitals as Hubbard projectors, we found considerable deviations in the band gaps, with the orthogonalized projectors showing the highest accuracy, in overall close agreement with experimental data.
Additionally, we have analyzed systematically which electronic states contribute predominantly to the valence band maximum and to the conduction band minimum in order to determine the states that might require Hubbard corrections. We have found that in the 14 TM-containing materials the d states of TM elements are mostly found at the bottom of conduction bands while the p states of light elements (O, N, S) appear at the top of the valence bands. We have not found any clear trend as to whether one should apply the U correction to both groups of states or only to the d states of TM elements-either of these two options perform equally (for the band gaps) or one of them is better than the other depending on the specific material. We note that all materials were considered to be nonmagnetic in our simulations, and hence future studies with more detailed investigation of various magnetic orderings would be useful. As for the remaining 6 materials containing the p-block (group III-IV) elements, we found that the p states of light elements (O, N, S) contribute mainly to the top of the valence bands and the p states of the group III-IV elements appear mostly at the bottom of the conduction bands. Application of the U corrections to both groups of states or only to the p states of light elements give essentially the same band gaps, which means that the application of U corrections to the p states of group III-IV elements leads to minor improvements (if any).
Hence, this study conclusively shows that DFT+U with ab initio U may serve as a consistent and reliable electronic-structure method for improving band gaps beyond (semi)local functionals, provided that care is taken in choosing appropriate Hubbard projectors.

Conflicts of Interest:
The authors declare no conflicts of interest.

Appendix A. Convergence tests for the self-consistent calculation of Hubbard parameters
In this Appendix we discuss how we performed convergence tests when computing Hubbard parameters using the DFPT approach (see Sec. 2.1). We did these tests for TiO 2 , NbNO, and Y 2 SO 2 which represent the material families considered in this study. As explained in Ref. [79], Hubbard U must be converged with respect to the k-and q-points sampling of the BZ: the former is the Bloch wavevector that is used to describe Kohn-Sham wavefunctions and charge density, while the latter is used to describe modulations of monochromatic perturbations that are used in linear-response calculations when computing U (see Sec. 2.1). As the criterium to monitor the convergence of U we use the ratio N k /N q , where N k and N q are the number of k-and q-points in the BZ, respectively. It is important to remark that we use this criterium with the condition that the k-point mesh for each material remains unchanged (i.e. it is fixed such that the spacing between the k-points is 0.04 Å −1 ). The result is shown in Fig. A1.
In these plots, convergence occurs from right to left, i.e. the smaller the ratio N k /N q the denser the q-points mesh at a fixed k points mesh. For example, N k /N q =1 is the case where the q points are equivalent to the k-points for a given material and hence, the most computationally expensive option (in this case the spacing between the q points is also 0.04 Å −1 ). Since we enforce a consistent k points sampling density in our calculations, k-points are unique for each material in this study. Thus, when we choose the number of q-points for the sampling of the BZ corresponding to the primitive cell, we look at the N k /N q ratio by which we divide the k-points for each material and truncate down to the nearest integer value. For example, if the k-points mesh is 6 × 5 × 5 and we have a N k /N q ratio of 2, then the resulting q points mesh is 3 × 2 × 2. The N k /N q ratios tested, as well as their subsequent q meshes for each material are shown in Table A1. From the q-point mesh testing, a N k /N q ratio of 2 was selected with a convergence threshold of ≤ 0.1 eV. By selecting a higher N k /N q ratio than 1, we are able to speed up Hubbard parameter calculations and reduce the computational cost.
Once the convergence of U is reached with respect to the q-point sampling of the BZ at a fixed k-point sampling, the next step is to perform the self-consistency convergence test. This is explained in detail in Ref. [38]. In short, it is necessary to reach self-consistency between the computed value of Hubbard U and the crystal structure; to do so, we need to perform structural optimization at the DFT+U level of theory with the current value of U, then recompute U using DFPT on top of the new geometry, then perform new structural optimization with the updated value of U, and so on until convergence is reached [38]. The results for the selected systems are shown in Fig. A2. It can be seen that the accuracy of ∆U ≤ 0.1 eV is reached after 3 cycles for all systems. Table A1. Comparison of q-point meshes for each material (TiO 2 , NbNO, Y 2 SO 2 ), corresponding to N k /N q ratios of 1, 2, 3, and 4. For the case where N k /N q =1, the q-point mesh is equivalent to the k-point mesh. Figure A1. Dependence of the Hubbard U parameter on the N k /N q ratio showing the results of N k /N q testing for 3 representative materials: (a) TiO 2 , (b) NbNO, and (c) Y 2 SO 2 . For each plot, the q points density increases from right to left, where at N k /N q =1, the q-point mesh matches the k-point mesh. In panel (c), Y 2 SO 2 does not have a N k /N q ratio of 4 plotted because its N k /N q ratio of 3 results in a q mesh equivalent to that produced using a N k /N q ratio of 4 (see Table A1).

Appendix B. Computational cost comparison of HSE and Hubbard parameter calculations
In this Appendix we present a comparison of the computational costs of hybrid functional (HSE) calculations and linear-response calculations of Hubbard U parameters [134]. This comparison is made for two representative materials, TiO 2 and NbNO. These materials were selected as examples of relatively small unit cells (rutile TiO 2 has 6 atoms per unit cell) and relatively large unit cells (NbNO has 12 atoms per unit cell) -in comparison to the diverse materials studied in this work. HSE self-consistent-field calculations were performed using the optimized geometry of the materials at the DFT+U level of theory. The computational setup is the same as described in Sec. 3. We recall that hybrid functional calculations are computationally very expensive because the calculation of the non-local exact-exchange potential contains double sums over electronic bands and double sums over k-points (one of these k-points meshes is replaced by a smaller q-points mesh to speed-up the calculations). As said above, it is common to lower the computational cost of HSE calculations (with, typically, minor losses in the accuracy) by reducing the number of q-points relative to the k-points of a given material by setting a ratio N k /N q higher than 1 [135]. For the current HSE calculations, N k /N q was set to 3 for each material [136]. The resulting k-and q-points meshes for the U linear-response and HSE calculations are listed in Table A2 (we stress that the q-points sampling has different physical meaning in these two types of calculations). The Hubbard parameters calculations reported in Table A2 were performed using a N k /N q ratio of 2, as discussed in Appendix A. It is clear from Table A2 that the HSE calculations are computationally more expensive than Hubbard U calculations, even in this particular case where a lower sampling density of q-points was used for HSE (N k /N q = 3) than for the Hubbard U calculations (N k /N q = 2) (even though these two q-points samplings are not really directly comparable). Additionally, it is seen from Table A2 that the relative computational cost increases rapidly with the system size, with the HSE calculations being 2.7 times more costly than Hubbard U calculations for TiO 2 and 8.2 times more computationally costly than Hubbard U calculations for NbNO. However, it is important to remark that in Table A2 we report the timing of only the one-shot calculation of the U parameters, while in practice we need to use the self-consistency procedure that requires typically 2-3 such calculations (see Appendix A). Thus, by considering the factor of 3 that was used in this work for the self-consistency loop, we obtain the ratios t HSE /t U of 0.9 for TiO 2 and 2.7 for NbNO. Hence, we can still see that DFT+U with self-consistent U values is more advantageous for systems with more than ∼ 10 atoms in the unit cell. In fact, in Ref. [36] it was shown that DFT+U with U computed using linear response theory can be used for systems as large as 320 atoms in the supercell, while this is beyond the reach for hybrid functionals with the current high-performance computing hardware.
We note that the comparison shown in Table A2 is very approximate as each of these calculations can be optimized and speed-up further (e.g. by lowering the kinetic-energy cutoff for the density and potentials for the exact-exchange term in the HSE calculations, by lowering the convergence threshold for the response matrices in the U calculations [see Eqs. (6)], by using more efficient parallelization strategies, etc.). Moreover, it is not only the computational cost that matters when comparing HSE with DFT+U calculations, but also the accuracy of the final results. In the fully first-principles DFT+U approach, U is a material-specific parameter and it is computed ab initio as a response property of the material. In HSE (and other popular hybrids), the amount of the exact-exchange (the screening parameter) is fixed (i.e. it is not material-specific), and in solids this often leads to unsatisfactory results and thus frequently this parameter is tuned until the experimental value of some property (e.g. the band gap) is reproduced, which makes this approach not fully first-principles. However, it is important to mention that this screening parameter in hybrids can be computed ab initio (see e.g. Refs. [13,15,18]), which is not always an easy task; in this case, an extra computational cost for HSE must be added in Table A2 which will make the ratio t HSE /t U even larger.
Therefore, the overall observations in this appendix further confirm that DFT+U can be considered as a reliable and efficient method to calculate band gaps at a fraction of the computational cost of hybrids functionals such as HSE (especially for large systems).