Efficient First-Principles Approach with a Pseudohybrid Density Functional for Extended Hubbard Interactions

For fast and accurate calculations of band gaps of solids, we present an {\it ab initio} method that extends the density functional theory plus on-site Hubbard interaction (DFT+$U$) to include inter-site Hubbard interaction ($V$). This formalism is appropriate for considering various interactions such as a local Coulomb repulsion, covalent hybridizations, and their coexistence in solids. To achieve self-consistent evaluations of $U$ and $V$, we adapt a recently proposed Agapito-Curtarolo-Buongiorno Nardelli pseudohybrid functional for DFT$+U$ to implement a density functional of $V$ and obtain band gaps of diverse bulk materials as accurate as those from $GW$ or hybrid functionals methods with a standard DFT computational cost. Moreover, we also show that computed band gaps of few layers black phosphorous and Si(111)-($2\times1$) surface agree with experiments very well, thus meriting the new method for large-scale as well as high throughput calculations with higher accuracy.

since other methods involve additional calculations. A recent proposal by Agapito-Curtarolo-Buongiorno Nardelli (ACBN0) [37] allows a direct self-consistent evaluation of U . They demonstrated improved agreements with experiments with a negligible increase in computational cost [37,38]. Second, inter-site Hubbard Coulomb interaction (V ) between the localized orbital of interest in DFT+U and its neighboring orbitals also need to consider properly because it could lead to better descriptions of electronic structures of some solids [39][40][41]. Moreover, DFT+U hardly improves LDA and GGA gaps of simple semiconductors such as Si while DFT+U with V does [40]. Therefore, combining these progress may realize a fast and accurate large-scale and high-throughput computational tool for materials with band gaps.
In this paper, we extend the ACBN0 functional for DFT+U [37] to implement a new density functional for the inter-site Coulomb interaction of V . With this, we achieve excellent agreements between self-consistent ab initio band gaps of diverse semiconductors and insulators and those from experiments. The band gaps comparable to those from GW and HSE methods [24][25][26][27] can be obtained within the standard DFT-GGA computational time. Moreover, for low dimensional systems in which the screening of Coulomb interaction varies significantly, the new method also can compute accurate band gaps of few layers black phosphorous and Si(111)-(2×1) surface, respectively, demonstrating its flexibility on structural and dimensional variations. We expect that this new approach could accelerate efficient high-throughput calculations with better accuracy for materials discovery.
Let us first consider mean-field (MF) energy of the Coulomb interaction between electrons in a pair of atoms I and J with the HF approximation, In the pairwise HF energy in Eq. (1), the general occu-arXiv:1911.05967v1 [cond-mat.mtrl-sci] 14 Nov 2019 pation matrix is written as where f nk is the Fermi-Dirac function of the Bloch state |ψ σ nk of the n-th band at a momentum k and w nk is the k-grid weight. The Löwdin orthonormalized atomic wavefunctions (|φ I i ) is used for projector of i-orbital in atom I. We note that the diagonal terms in Eq. (2) are the usual on-site occupations for DFT+U . In Eq. (1), we neglect other small pairwise interactions, e.g., the cross charge exchanges between neighbors [40] and discuss their effects in Supplementary Information (SI) [42].
Assuming the effective interactions of V ee in Eq. (1) are all equal to their atomic average [40], the rotationally invariant form of E MF can be written as E MF = E Hub = E U + E V where E U is for the case of I = J and E V for I = J. E U is a well-known density functional for U by Dudarev et al. [43]. For I = J case where atoms I and J locate at different positions, respectively, where the {I, J} indicates the summation for a pair of atoms I and J of which interatomic distance of d IJ is less than a given cutoff. In Eq. (3), V IJ is the inter-site Hubbard interaction for the pair and will be determined based on the method of ACBN0 [37].
To obtain a functional form of V IJ , as is discussed in ACBN0 [37], we also follow a central ansatz by Mosey et al. [35,36] that introduces a "renormalized" occupation number for the pair such as where the sum is divided by 2 if the pair of I and J is the same type of atoms. Corresponding to ACBN0 functional for U , we can replace n IJσ ij in Eq. (1) by a renormalized density matrix for the pair, Simultaneously, the bare Coulomb interaction between electrons belong to the pair can be expressed by the electron repulsion integral [37], where i and k are orbital indices belong to atom I and j and l to atom J. Using Eqs (5) and (6), the ACBN0-like energy expression (E V ACBN0 ) for the inter-site Hubbard interaction can be written as where the additional prefactor of 1/2 is for a double counting of same orbitals. Equating Eq. (3) to Eq. (7), then we can obtain a density functional for V IJ , ] .
(8) An energy functional for V can be constructed by subtracting a double counting term (E dc Following the discussion in Ref. [40], we use a fully localized limit and E dc With this, the final functional for inter-site V interaction can be written as where n IJσ is the matrix notation for the general occupation in Eq. (2), {n} = n IIσ , n IJσ and V IJ [{n}] in Eq. (8).
For the on-site repulsion, we used the ACBN0 functionals in Eqs (12) and (13) in Ref. [37] so that we complete a construction of pseudohybrid type functionals for the two important Hubbard interactions. As discussed in Ref. [40], the minus sign in E V [{n}] highlights the role of inter-site Hubbard interaction that localizes electron between atoms, thus improving description of bondings between the atoms or augmenting 'overlocalization' [41] caused by U in case that the bonding between neighboring d-and p-orbitals plays an important role.
We implemented E V [{n}] in Eq. (9), ACBN0 functionals and other related quantities in Quantum ESPRESSO package [44]. For the Kohn-Sham potential corresponding to Eq. (9), we used Eq. (13) in Ref. [40]. To compute V ERI in Eq. (6), we used the PAO-3G minimal basis set as in Ref. [37]. With the aid of PyQuante package [45], the integrals were done quickly in an analytical way. For all calculations here, the cut-off for d IJ sets within the second-nearest neighbors. The on-site interactions for s-orbitals were neglected. Fully converged U and V were obtained when the difference between the energies in the self-consistent loop is less than 10 -8 Ry. We used the GBRV ultrasoft pseudopotentials [46] for all materials considered. Regarding on pseudopotential dependence of ACBN0 functionals [38], we tested another set of the norm-conserving pseudopotentials provided by PseudoDoJo project [47] and discuss its effects in SI [42]. The kinetic energy cutoff was set to 160 Ry and the Brillouin zone (BZ) integration were performed with a Γ-centered k-point grid spacing of 0.2Å -1 . The lattice structures are from the experimentally available data.
We first tested our method for selected bulk solids with diverse characteristics. Table I and Fig. 1 summarize the calculated band gap of 23 solids. We also list the results from other methods and experiments. We select solids from group IV and III-V semiconductors, an ionic insulator and metal chalcogenides and oxides. As shown in Table I, our calculated band gaps are in excellent agreement with experiments and are as accurate as those from HSE and GW methods. We may infer from mean absolute relative error (MARE) with respect to experimental data that our method, HSE and GW methods are closer to experiments than PBE and ACBN0. We also check trends Ref. [25] except GaP, InP, AlAs from scGW in Ref. [54], ZnSe from a partially self-consistent GW (GW 0 ) in Ref. [55], BP from GW 0 in Ref. [56], ZrO 2 from scGW in Ref. [57], TiO 2 from scGW in Ref. [58], and Cu 2 O from scGW in Ref. [59]. d Monoclinic structure from mean relative error (MRE) that ACBN0 and PBE underestimate the gaps (minus sign) while GW methods overestimate. Hereafter, we mainly focus on calculated gap values and, for future references the band structures of all solids are displayed in SI [42].
For the group IV semiconductors, the effects of U on the band gaps are almost negligible as shown in Table 1 (see ACBN0 column) while the inter-site Hubbard terms improve the band gaps dramatically as was also discussed in a previous study using the linear response theory [40]. For the group III-V semiconductors, both U and V affect their electronic structures because of their mixed covalent and ionic bonding characters. So, ACBN0 improves PBE gaps and the inter-site V increases these further to match experiment values. Details of computations such as self-consistent U and V for Si and GaAs compared with Ref. [40] are presented in SI [42]. We note that PBE and ACBN0 wrongly describe Ge and InAs to be metallic and to have an inverted band gap (or topological insulator), respectively, while our method confirms them as semiconductors like HSE and GW results.
In case of an ionic compound LiF, the on-site U improves the PBE band gap significantly because of its strong local Coulomb repulsion. Nonetheless, the intersite V still increases the ACBN0 gap further to agree with its experiment value very well.
A similar trend also can be found for metal monochalcogenides (here Zn compounds only). For these compounds, U and V functionals play similar roles as they do for LiF so that our results with U and V are quite closer to experiment values than those with U only. We note that the calculated gaps depends on choice of pseudopotential of Zn while there is no such a dependence in cases of IV and III-V semiconductors. We will ). Because effects of on-site and inter-site interactions depend on degree of localization or cut-off in projector for localized orbital [38], it is important to select or generate pseudopotentials with care to obtain accurate results [61].
For Cu 2 O and Zr 2 O, our results are comparable to those from HSE and GW calculations. Considering limits in mBJLDA to obtain gaps for these compounds [29], our method could be a good alternative fast tool for studying zirconia and cuprite. We note here that for Cu 2 O the calculated energetic position of fully filled d-orbitals is lower than that of HSE value [53] by ∼ 0.9 eV [Fig. S7 in SI [42]]. So, it needs to improve the way to treat the completely filled d-orbitals with U and V . Now, we consider low dimensional systems where the screening of Coulomb interaction varies rapidly. The GW approximation can calculate quasiparticle gaps quite accurately but its convergence is very slow with respect to the k-points density and other parameters [62,63]. The hybrid functional methods do not suffer such a issue but they produce unreliable gap values with structural or dimensional variation [64]. The mBJLDA, another low-cost alternative for bulk solids, also suffer a similar problem with hybrid functionals [65]. Our method selfconsistently computes occupations of atoms at boundary and bulk and their influence on screening of Coulomb interaction through Eqs. (4), (5) and (8). So, we expect that the current method may overcome the aforementioned difficulties for low dimensional materials.
To test the new method, we first calculated electronic structures of Si(111)-(2 × 1) surface. Because of a unique surface reconstruction resulting in a quasi-onedimensional π-bonded chain of Si p z -orbitals [68] and because of large difference between screenings on surface and in bulk, it is a good test bed for a method to compute surface and bulk gap simultaneously [64,69]. A 24-layer slab with ∼15Å vacuum was optimized with GGA. The kinetic energy is set to 80 Ry and d IJ to the nearest neighbors. The surface has two degenerate and coexisting reconstructions [70,71] Fig. 2, the calculated averaged surface gap is 0.83 eV, agreeing well with experimental value of 0.75 eV [66,67], together with an accurate bulk gap. We note that the converged Hubbard interactions of Si atoms change spatially, reflecting local variation of screening such that converged U and V are confirmed to gradually increase from inside to the surface.
Next, few layers black phosphorus were chosen to test our method [Fig. 3]. We use fully relaxed crystal structures using a rev-vdW-DF2 functional [75,76]. All intersite interactions between valence s and p electrons of P atom within the plane are considered. Fig. 3 shows the calculated band gaps as a function of a number of layers, together with results from other methods and an FIG. 3. Band gaps of black phosphorus as a function of a number of layers. We also list other gaps from PBE, HSE06 [65], mBJLDA [65], GW0 [72], G0W0 [73] and GW -BSE [73]. A experimental band gap [74] is denoted with black dots. experiment. It is noticeable that without including V , all ACBN0 gaps are quite smaller than those from GW , and that the gap values from HSE [65], mBJLDA [65] and ACBN0 are close to optical gaps from GW -BSE method [73]. Considering qualitative difference in shape of optical spectrum between GW -BSE and HSE or other hybrid functionals [64,77], we conclude that they underestimate band gaps. As shown in Fig. 3, our results are consistent with GW results [72,73] and an available experiment [74]. Unlike bulk solids, nanostructures typically have a large number of atoms for which expensive calculations hardly apply so that our method will play an important role in calculating their accurate band gaps.
In conclusion, we report a new ab initio method for electronic structures of solids employing a pseudohybrid density functional for extended Hubbard Coulomb interactions. We demonstrate that the new method significantly improves the ACBN0 functional in obtaining band gaps of bulk and low dimensional materials. Its selfconsistent calculation for accurate band gap can be done with a computational time comparable to DFT-GGA. With further validations and improvements of the current method such as noncollinear spin and forces [38], our new method could fulfill requirements [6] for first-principles simulations suitable for massive database-driven materials research with improved accuracy.

Effects of cross exchange interactions between orbitals
The mean field expression of the electronic interaction energy in terms of atomic orbitals in its most general form can be written as , where, The generalized occupation matrix in which its diagonal terms are the on-site occupations is given by where w k is the k-point weight and f nk is the occupation of the Bloch state |ψ σ nk of the n-th band at momentum k. The |φ I i is an orthonormalized atomic projector of i-orbital in atom I. Considering E IJKL ijkl,σσ , there are many possible arrangements for IJKL and ijkl [40], respectively. Among them, here we consider the first three large contributions, E IIII ijij,σσ , E IJIJ ijij,σσ and E IJJI ijji,σσ where I = J. The first and second terms were discussed in the main manuscript and corresponds to the on-site and the inter-site Hubbard interactions, respectively. The last one is the cross charge exhanges between the neighboring atoms I and J [40]. The first case where all IJKL are equals will lead to the well known Hubbard density functional for LDA+U method and if we use a rotationally invariant on-site interaction (E I U ), the Eq. (S1) will become to be Dudarev U functional [43]. The second term becomes the inter-site Hubbard interaction as discussed in the main manuscript. Now, we consider the second and third ones simultaneously. If we use rotationally invariant forms for the interaction of φ I i φ J j |V ee |φ K k φ L l in Eq. (S2), then we can rewrite the second interaction using where N IJ is a number of degeneracy of angular momentum for atoms I and J [40]. Likewise, the third interactions, we can use With these considerations, the MF energy in Eq. (S1) can be written as where the {I, J} indicates the summation for a pair of atoms I and J of which interatomic distance of d IJ is less than a given cutoff.
Using the matrix notations of n IJ = σ n IJσ and n I = σ n Iσ = We assume a fully localized limit for double countings as was also discussed in a previous work [40] so that the final expression for Hubbard pairwise energy is given by, where we neglect K IJ σ =σ Tr[n IJσ n JIσ ] thanks to V IJ − K IJ K IJ . If we compare Eq. (S6) with the inter-site Hubbard functional shown in Eq. (9) in the main manuscript, we immediately notice that a replacement of V IJ by V IJ eff = V IJ − K IJ is enough for including the effects of cross charge exchange between the pair of atoms I and J. Now using Eqs. (4)-(6) in the main manuscript, a pseudohybrid functional or ACBN0 functional expression for K IJ can be obtained in a straightforward way and the final expression can be written as The effects of cross exchange interactions on band gaps are summarized in Table S1. As shown in Table S1, the calculated band gaps with and without effects of K are negligile. For solids considered here, all computed gaps with K are little bit smaller than the gap without K.

Effects of choice of pseudopotentials
In this work, we used two kinds of pseudopotentials: PseudoDoJo norm-conserving [47] and GBRV ultrasoft [46] pseudopotentials. Table S2 shows the effects of choice of pseudopotentials. Almost every s and p electron systems are not affected by the choice. However, particularly, in the case of Zn compound, the discrepancies are quite large.

Comparison of U and V for Si and GaAs
We present the calculated U and V values for Si and GaAs to compare with the values obtained by the linearresponse approach in Ref. [40]. Table S3 shows the calculated results. Despite of the differences in the Hubbard interaction terms, the calculated band gap of Si (1.34 eV) is in a good agreement with the value (1.36 eV) in Ref. [40]. In the case of GaAs, however, the on-site interactions are significantly lower than the previous results and the band gap (1.20 eV) is much close to the experimental gap rather than the value (0.90 eV) in Ref. [40].

Comparison of U and V for transition metal oxides
To compare with the previous studies [37,38], we present the obtained on-site Hubbard interaction terms U d (for d electrons of transition metal) and U p (for p electrons of oxygen) for TiO 2 , MnO, NiO and ZnO using ACBN0 method. We also provide those values and the first-nearest neighbor inter-site Hubbard interaction terms (the d-p interactions) calculated with our method. The results are summarized in Table S4. As mentioned in Ref. [38], the discrepancies of U d and U p is mainly attributed to the way to calculate Coulomb integrals and the treatment of the localized orbitals. Moreover, we show that pseudopotentials also play crucial role to make the discrepancies.

Calculated band structures of selected bulk solids
In Figs. S1-S8, We calculated the band structures of selected bulk solids given in the main manuscript. In order to show the effects of inter-site interactions V on the band structures, we present the results obtained by ACBN0 method (gray dashed lines) and our method (red solid lines).
Surface band gaps of two types of buckled Si(111)-(2 × 1) The cleaved Si(111) surface undergoes a 2 × 1 reconstruction resulting in a quasi-one-dimensional π-bonded chain which is buckled to decrease the total energy [68,70]. There are two energetically degenerate buckled structures as shown in Fig. S9 (a) and (b). The experimental observation clearly shows that the coexistence of positively and negatively buckling on Si(111)-(2 × 1) surface at broad range of temperatures (6 K ∼ room temperature) [70]. The previous density functional study also demonstrated that the total energy difference between the two buckled types is negligible [71]. In the main manuscript, therefore, we provide the band structures of the two buckled types all together. The Fermi level is set to the bulk valence band maximum at Γ in both cases. We found that the negative surface states (E g,negative = 0.72 eV) are in between the positive surface states (E g,positive = 0.94 eV). The experiment reported that the energetic differences between the surface valence bands of two buckled structures are 129 (10) meV and the surface conduction bands 99 (10) meV [70], which is in a good agreement with our calculations (103 meV and 117 meV, respectively). S1. Calculated band gaps (in eV) with and without the cross charge exchange in Eq. (S7). '+V ' column summarizes the band gap with V IJ and without K IJ while '+V eff 'column with V IJ eff = V IJ − K IJ . Experimental data for energy gaps are from Refs. [28,29,37,48,78,79]. a All data from self-consistent GW (scGW ) calculation results in Ref. [25] except GaP, InP, AlAs from scGW in Ref. [54], ZnSe from G 0 W 0 in Ref. [55], ZrO 2 from scGW in Ref. [57], TiO 2 from scGW in Ref. [58], and Cu 2 O from scGW in Ref. [59].