Wannier–Koopmans method calculations for transition metal oxide band gaps

The widely used density functional theory (DFT) has a major drawback of underestimating the band gaps of materials. Wannier–Koopmans method (WKM) was recently developed for band gap calculations with accuracy on a par with more complicated methods. WKM has been tested for main group covalent semiconductors, alkali halides, 2D materials, and organic crystals. Here we apply the WKM to another interesting type of material system: the transition metal (TM) oxides. TM oxides can be classified as either with d0 or d10 closed shell occupancy or partially occupied open shell configuration, and the latter is known to be strongly correlated Mott insulators. We found that, while WKM provides adequate band gaps for the d0 and d10 TM oxides, it fails to provide correct band gaps for the group with partially occupied d states. This issue is also found in other mean-field approaches like the GW calculations. We believe that the problem comes from a strong interaction between the occupied and unoccupied d-state Wannier functions in a partially occupied d-state system. We also found that, for pseudopotential calculations including deep core levels, it is necessary to remove the electron densities of these deep core levels in the Hartree and exchange–correlation energy functional when calculating the WKM correction parameters for the d-state Wannier functions.


INTRODUCTION
Transition metal (TM) oxide is one of the most important class of materials used in electronics, solar cells, catalysts, batteries, and many other devices. Band gap of the material is an essential characteristic for these applications. It is, thus, of paramount importance to be able to predict the band gaps of the material by using first principle calculations. However, it is well known that the density functional theory (DFT) approximated by local density approximation (LDA) 1 or generalized gradient approximation (GGA) 2 and implemented with Kohn-Sham equations significantly underestimates their band gaps due to the lack of discontinuity in their functional derivatives. There are methods utilizing empirical parameters to derive correct band gaps, such as meta-GGA 3 , HSE06 4 , B3LYP 5,6 , and Hubbard U methods 7 . However, the empirical nature of these methods leaves a lot to be desired. A more standard first principle approach requiring no empirical parameters to calculate the band gap is the many-body perturbation GW method. However, even for the GW method 8 , getting the correct TM oxide band gap is a challenge. Not only it can have different flavors in regard to the self-consistency, including G 0 W 0 , GW 0 , and GW (here 0 means no self-consistency), which can give different results, but also the initial input wave function and eigenenergies can also influence the final results for the not fully self-consistent calculations. Moreover, there are also challenges for numerical implementations (e.g., plane wave basis versus localized basis set) and other approximations (e.g., plasmon pole approximation for the dielectric function). Besides, some of the calculations require a large number of conduction bands to yield fully converged results, which can lead to debates about the accuracy of the GW method. Thus, even if fully converged GW method might be good, there could be many problems in practical use of its calculations either due to numerical implementation issues or convergence issues. For example, for ZnO, the experimentally measured band gap is 3.43 eV 9 , while the reported GW calculated band gap ranges from 2.1 to 3.9 eV 10 . The situation is even worse for TM oxides with partially occupied d states. For rock salt MnO system, the experimentally measured band gap is believed to be 4.0 eV 11 15 . Jiang reported a G 0 W 0 result with highenergy local orbitals and linearized augmented plane waves of 3.32 eV, and the GW 0 result is 3.69 eV in 2018 16 . Moreover, Carter tested G 0 W 0 , GW 0 , and GW, on PBE, LDA + U, PBE + U, HSE06, and PBE0 for Fe 2 O 3 in 2011, and the results vary from 1.3 to 4.8 eV compared with an experimentally measured band gap of around 2.6 eV 17 . On top of all these, GW calculation can be expensive, especially when a large number of conduction band states are needed or when the system size is large. Given this situation, it will be very interesting to find whether there are alternative parameter-free methods to predict the TM oxide band gaps.
There have been many prior works in using Koopmans condition to improve the accuracy of band gap of system especially for atoms and small molecules. For example, Dabo et al. 18 proposed a generalized Koopmans theory to remove not only the electrons from the highest valence states but also the electrons from other states to correct the energies of these states. They applied this method to atoms and small molecules and achieved excellent agreement with the experimental results. Cococcioni et al. 19 used Koopmans condition to fix the U parameter in DFT + U calculations. Kraisler et al. 20 showed that it is possible to use Koopmans condition to correct any exchangecorrelation functional, restore the discontinuity of its derivative, and used that to predict the electron affinities and ionization energies of atoms and small molecules. Over the past few years, we have developed a Wannier-Koopmans Method (WKM) which can yield accurate band gaps on main group covalent semiconductors 21 , alkali halides 22 , 2D materials 23 , and organic crystals 24 . For these tested materials, the accuracy of WKM is similar to that of GW method. In our approach, Wannier functions and the Koopmans condition are combined to calculate the band gap of bulk systems. Similar methods have also been adopted in recent years. For example, Yang et al. proposed a localized orbital scaling correction framework 25 . Instead of implementing a post-DFT correction method like ours, they applied it during the selfconsistent (SCF) calculations. However, the main purpose is similar: to use a localized function in compliance with the Koopmans condition to correct the DFT band gap. Marzari et al. reported a Koopmans-compliant spectral functional method 26 . For the conduction band, they used an approach similar to ours with localized Wannier functions. For valence band, they used the selfinteraction correction (SIC) function to derive the localized orbital. Since SIC localized orbital and maximally localized Wannier functions are very similar, their results are also very similar to WKM results. Given all the interest in this approach, it is thus important to test the limit of WKM. One class of important materials which have not yet been systematically tested is the TM oxides. The TM oxides can be classified into two groups. One group is with d 0 or d 10 configurations for the TM elements, and the other group is with partially occupied d states. The second group of TM oxides have long been regarded as strongly correlated Mott insulators, which are difficult to be described using mean field-like approaches (e.g., DFT, SIC, or even GW). It will, thus, be interesting to see whether there are qualitative differences between these two groups of TM oxides when the WKM method is applied.
In our following calculations, we found that the WKM method works well for the d 0 or d 10 close shell systems, whereas it performs badly for the partially occupied open shell configurations. The main cause of this problem for the open shell system is due to a strong interaction between the occupied and unoccupied d-Wannier functions as we will show later. This is perhaps compounded by the fact that the DFT ground-state properties (e.g., the magnetic moment) of the open shell systems can also be significantly wrong.
WKM is an extension of the ΔDFT method 27 . In the traditional ΔDFT method, the band gap is calculated by using the difference between the electron affinity (EA) energy and the ionization energy (IE). EA and IE can be calculated by using self-consistent ground-state energy E(N + 1), E(N), and E(N − 1). Here, N is the number of electrons in the neutral system, and N + 1 and N − 1 indicate that the system has one more or one less electron, respectively. EA can be expressed as E(N + 1) − E(N) and IE as E(N) − E(N − 1). Although this method works well for isolated molecules, the density change caused by adding or subtracting one electron is infinitesimally small in an extended system. As a result, according to Janak's theory 28 , the total energy difference is the same as Kohn-Sham orbital eigenenergy. To overcome this problem, we added an electron into a localized Wannier function instead of the extended Kohn-Sham orbitals. The resulting LDA total energy with partial occupation s w of the Wannier function ϕ w can be denoted as E LDA ({s w }). This function is not a linear function of s w between 0 and 1. However, in many-body quantum mechanics, the total energy of a system with a fractional number of electrons can be defined as a statistical mixture of the N electron and N ± 1 electron state. This leads to a linear segment total energy function of s w . This linear segment property is also called Koopmans condition. More deeply, this means that the system (or a relatively isolated part of the system) likes to occupy an integer number of electrons, instead of partial number of electrons. This would be a consequence of a grand canonical minimization using the linear segment total energy function with a given external Fermi energy (a consequence of a linear programming optimization). The integer number of electrons in a given system is also a property of many-body electron wave treatment of the system. Finally, to satisfy the Koopmans condition, one compensation term E w (s w ) can be added to the total energy expression. As a result, the WKM total energy can be expressed as Here w indicates a set of Wannier functions that are orthogonal to each other. s w (0 < s w < 1) indicates the occupation number of this Wannier function ϕ w . E w (s w ) is a simple function of s w obtained by requiring that E WKM ({s w }) is a linear function of s w . During the LDA calculation, we add or remove the electron from ϕ w of one spin channel and all the other orbitals in this spin channel should be orthogonal to this Wannier function ϕ w . All the other orbitals (except this one Wannier function) are variationally changed to minimize the total energy, which results in the ground-state energy E LDA ({s w }).
A simple analytical expression of E w (s w ) can be written as The λ w can be determined from E LDA ({s w }) (to make E WKM ({s w }) a straight line versus s w ). The E LDA ({s w }) calculation is done using a supercell that should be large enough to remove the image interactions between neighboring supercells. In our calculations, we used sufficiently large supercells to reduce the image charge problem. Based on a test in our previous work 23 , the image charge effect on band gap is likely to be in the order of 0.01 eV when the distance between Wannier functions from the neighboring supercell is >10 Å. Note, in the future, image interaction correction techniques can be used to reduce the needed supercell sizes, like that in a charged impurity defect calculation.
After doing the self-consistent calculation for E LDA ({s w }) with a Wannier function ϕ w , we can obtain λ w . For a Kohn-Sham orbital ψ i , we can have the modified eigenenergy ε i with More details of this method can be found in ref. 21 .  Fig. 3 and Supplementary Note 1). We used lattice constants from the Materials Project relaxed result for all these systems 29 . The details of these structures are listed in the Supplementary Materials.

RESULTS AND DISCUSSION
WKM calculations were performed in two different sets of norm conserving pseudopotentials: the FHI 30-34 pseudopotential and the recent SG15 pseudopotential 35,36 . The FHI pseudopotential does not include semicore states, while the SG15 includes the semicore states. Their detailed valence electron configurations are included in Supplementary Table 3. In general, the semicores are defined as the subshell levels that are not in the row of the periodic table where this element resides. For example, for the fourth-row elements (e.g., Fe), the 3s and 3p levels will be semicores. However, we have also taken O-2s level as the semicore as it is often isolated 10 eV below the other part of valence bands near the band gap.
For the LDA calculations, the SG15 and FHI pseudopotential yield similar band gaps. Here, for clarity, we will call the straight forward WKM procedure described in ref. 21 , as the "normal" WKM procedure. This is to distinguish it from a modified WKM procedure where the semicore electron charge density is excluded from the Hartree and exchange-correlation energy (E hxc ) calculations (as will be discussed later). The normal WKM calculation shows that, in most d 0 and d 10 cases, the SG15 gives a larger WKM band gap than FHI, and it can be 1 eV too big compared with the experiment. The results are shown in Fig. 1 labeled as the blue squares and SG15-WKM. In stark contrast, the FHI-WKM results are close to the experimental values as shown in red circle.
We believe that the difference between the SG15 and FHI-WKM band gaps comes from the interaction between the semicore-level electrons and the d-orbital Wannier function electrons via the Hartree and exchange-correlation energy. In a sense, if we believe that the band gap correction should come from the valence band properties, then the correction parameter like λ should also be able to be calculated by the valence bands alone. If the inclusion of the semicore level in the calculation makes a major difference, then we should trust the original calculation based on the valence bands alone. This argument leads us to the following procedure to exclude the semicore electron charge density from the E hxc calculation.
The original Hartree and exchange-correlation energy E hxc is calculated as Here ρ c , ρ v , and ρ w stand for semicore, valence, and Wannier orbital charge densities, respectively. Now, in a revised version, we have: This expression is used when calculating E LDA ({s w }) of Eq. (1) as well as for total energy minimization when ϕ w is partially occupied with s w in the WKM SG15 pseudopotential calculations. During the SCF iterations, we have also kept the semicore-level wave function ψ c unchanged from their original LDA wave function inputs. Physically, this procedure is equivalent to excluding the semicore levels from participating in the screening of the Wannier function charge. Usually, one might think such screening should reduce the charge self-interaction, thus reducing the amplitude of λ w . Surprisingly, in reality, we found that including such semicore can increase the amplitude of λ w , probably due to a highly nonlinear exchange-correlation energy contribution. This might come from the error of the LDA exchange-correlation function. It thus also partially justifies the removal of the semicore term from the E hxc term. At this stage, however, we do treat this procedure of removing the semicore level from exchange-correlation energy as an empirical procedure.
Let us now use anatase TiO 2 as an example to discuss this procedure in more detail. The LDA band structure and projected density of states of anatase TiO 2 is shown in Fig. 2a, b. The LDA band gap of anatase TiO 2 is 1.8 eV. The calculations are performed using SG15 pseudopotentials. In the SG15 pseudopotentials for Ti and O, we have several semicore states. Since Ti is a fourth-row element, the Ti-3s and Ti-3p levels are treated as the semicore, and their corresponding band structure and density of states are shown in Fig. 2a, b. We then generate Wannier functions for anatase TiO 2 by wannier90 code 37 . For valance band part, we used initial guess of p orbital on O atoms; for conduction band part, we used initial guess of d orbital of Ti atoms. The Wannier functions with the maximum project weights on conduction band minimum (CBM) and valence band maximum (VBM) are shown in Fig. 3 within the supercell used for the λ w calculation.
We then remove different numbers of semicore electrons in the E hxc calculation. The results are shown in Table 1. The experimental band gap for anatase TiO 2 is 3.2 eV 38 . The more semicore states are removed in the E hxc , and the smaller band gaps are obtained in WKM. We also notice that the main difference comes from the λ w from the conduction band.
To test whether this semicore removal from E hxc procedure is necessary for other systems where there is no d-state Wannier functions, we also calculated the AlAs. The experimental band gap for AlAs is 2.2 eV 39 . In the Wannier generation, we used As-p orbital for valance bands and Al-s for conduction bands. The results are shown in Table 2. As we can see, for this system, the effects of the semicore level are rather small. Thus, we conclude that, it is not necessarily to have this E hxc semicore-level removal procedure when calculating the λ w of non-d-state Wannier functions.
We also used GGA for E hxc in our calculation. A trend similar to LDA is shown in Supplementary Table 4. Lastly, we expect the same semicore removal procedure should also be used for allelectron calculations.
The results of using this E hxc semicore-level removal procedure for the d 0 and d 10 TM oxides using SG15 pseudopotentials are shown in Fig. 1 as green diamond (labeled as SG15-WKM-fix). Now, we can see that the SG15 and FHI results are very similar, and they are all close to the experimental band gaps. The band gap values of SG15 calculation, SG15-WKM-fix calculation, FHI calculation, and experimental result can be found in Supplementary Tables 5-8, respectively. We, thus, conclude that the WKM works well for the d 0 and d 10 TMs after this core-level removal procedure in E hxc . It is interesting to note that most of the WKM-predicted band gaps are still slightly larger than the experimental results. It will be interesting to test in the future, whether the electron-phonon coupling effect, which tends to reduce the band gap 40 , can bring the theoretical results even closer to the experimental values.
WKM calculations for partial occupied d-orbital TM oxide Having finished the calculations for d 0 and d 10 TM oxides, we now turn our attention to the partially occupied d-state oxides. In this      Fig. 4, while the original data can be found in Supplementary Table 8. Because FHI-LDA gaps are not consistent with SG15-LDA gaps for these systems, and their ground-state magnetic moments are also very different (see Supplementary Tables 13 and 14), here we have only considered the SG15 results.
During the E LDA ({s w }) calculation of these open shell systems, we found that the number of electrons in a given spin can change dramatically after s w electron is placed inside one Wannier function. Very often, the loss (or increase) of the charge from one spin will cause a charge transfer between spin channels. This might be caused by a strong on-site interaction between the occupied d-Wannier (constructed from the valence band) and unoccupied d-Wannier (constructed from the conduction band). For example, in the ground state of NiO (an antiferromagnetic system in LDA calculation), the number of electrons in spin-up and spin-down channels are the same (both are 1536 electrons in a supercell). However, after adding one electron to the unoccupied d ϕ w of spin-down channel in a WKM supercell SCF calculation, the converged ground state will have 1534.56 electrons in the spin-up channel and 1538.44 electrons in the spin-down channel (instead of the expected 1537 electron if there was no spin flow between these two channels). Corresponding to this spin flow, the system lost its band gap, and became metallic (judged by the eigenenergies for orbitals besides the Wannier function). Such spin flow between the spin-up and spin-down channels does not exist in the WKM calculation of d 0 or d 10 systems, like ZnO.
The above strong on-site interaction between the opposite spin Wannier functions is a manifestation of the strongly correlated system. To avoid this spin flow between opposite spins, for all the E LDA ({s w })calculation below, we have forced the number of spin-up and spin-down electrons to be unchanged in the self-consistent iteration (so spin up and spin down will have different Fermi energies).
The tested results are shown in Fig. 4. We have tested the normal WKM procedure (shown as SG15-WKM), and the procedure removing the semicore in the E hxc calculations (shown as SG15-WKM-fix1), both have forced the number of spin-up and spindown electrons unchanged during SCF calculations. As one can see, although removing semicore procedure does reduce the WKM gap as it did on the d 0 and d 10 materials, it can make the result even worse when compared with the experiments. For Cr 2 O 3 , Fe 2 O 3 , and NiO, the WKM band gaps with fixed/excluded semicore using SG15 are even smaller than that of LDA (see also Supplementary Table 9). Thus, fixing the spin-up and spin-down electron number during WKM SCF calculation does not solve the problem.
The remaining problem can be analyzed by looking at the magnetic moment at each atom during WKM SCF calculations (see Supplementary Tables 10 and 11). Taking MnO as an example, the Mn LDA magnetic moment for SG15 pseudopotential calculation is −4.3 μ B . In Fig. 5a, we show that, in both SG15-WKM and SG15-WKM-fix1 calculations (where the numbers of electron in the spin-  up and spin-down channels are fixed, so there is no spin flow between these two channels), after removing one electron from the valence d x2−y2 Wannier function (shown in Supplementary Fig.  5) in the spin-down channel, the absolute value of magnetic moment for the atom at the center of the Wannier function has changed from −4.3 to −2.5 μ B . This reduction is >1 (the removal of the electron from the down spin). This is caused not by a spin flow between the spin-up and spindown channels, as their numbers are fixed. Rather, it is caused by a spatial charge flowing within the same spin, mostly from the opposite spin (the up spin). More specifically, the charge in the up spin, from the neighboring Mn atoms from the occupied Wannier function, has flowed to the unoccupied up spin Wannier function at the center. In another word, there is an interaction between the occupied and unoccupied d-state Wannier functions between neighboring atoms at the same spin. This is shown in Fig. 5b, where the spin-up electron charge density after the SCF E LDA ({s w }) calculation increases at the center atom and decreases around the neighboring atoms. Such charge flow is driven by the local total charge (spin up plus spin down) neutrality requirement to minimize the total energy and enabled by the interaction between the conduction and valence band d-state Wannier functions. This is not the case for the Wannier functions of other orbitals. For instance, in NiO, after removing one electron from the O-p x orbital Wannier function, the magnetic moment changed from 0 to 0.81 and the change is <1.
The above spatial charge flow can also be analyzed directly using Wannier occupation o w ¼ P i ϕ w h jψ i i j j 2 o i , instead of atomic magnetic moments. Before the WKM SCF calculations (e.g., using LDA ψ i ), o w is 0 for conduction band Wannier function ϕ w , and 1 for valence band Wannier function ϕ w . If we take one electron from one valence spin-down Wannier function ϕ w , then the o w of the neighboring valence spin-up Wannier functions will become less than 1 after SCF calculation, and o w for the center conduction band spin-up Wannier function will have o w > 0 (e.g., 0.36 for NiO). This corroborates well with the atomic magnetic moment analysis.
To avoid such opposite spin spatial charge flow, we carried out a test where the wave functions as well as their occupations of the opposite spin are fixed at their ground-state LDA values when performing the SCF E LDA ({s w }) calculations. We call this the "fix2" procedure in Figs. 4 and 5. The "fix2" procedure has indeed resulted in smaller than one magnetic moment changes when we remove (or add) one electron from one d-state Wannier functions, as shown in Fig. 5a (see also Supplementary Table 12). Fixing the opposite spin wave function (not just its total number of electron) does increase the band gap, but it overestimates the band gap significantly. This is because by fixing the opposite spin wave functions, we have completely removed their dielectric screening effect, which lead to a too big self-interaction energy for the Wannier function, thus too large λ w values and band gaps. So far, we have not found a reliable way to calculate the band gap of these Mott insulators under the current WKM approach. Compared with the d 0 and d 10 systems, one challenge here is the strong interaction and coupling between the two spin channels via their d-state Wannier functions from the opposite sides of the band gap. In a more conventional semiconductor, or in the d 0 and d 10 cases, such strong interaction does not exist.
The difficulty of these Mott insulators might also be related to the DFT ground states themselves. It could be true that these DFT ground states (e.g., the occupied electrons wave function subspace) are not correct. This can be viewed by their atomic magnetic moment (which represents the amount of d-state occupation). Supplementary Table 14 shows that the DFT calculated magnetic moment can be very different from the experimental value. For instance, the magnetic moment in LDA calculation for NiO is 1.1 μ B , whereas the experimental value is from 1.64 to 1.9 μ B . Since our current WKM method does not mix the valence band wave function subspace with conduction band wave function subspace, it is unlikely that it can fix this DFT ground-state problem.
Given the strong coupling between the occupied and unoccupied d-state Wannier functions, and the need to correct the DFT ground-state occupation subspace, it seems like, to find a solution for the partially occupied d-state systems, it might be necessary to use Wannier functions not constructed from occupied and unoccupied orbitals separately, instead d-state Wannier functions constructed from both occupied and unoccupied orbitals might be necessary. Such Wannier functions can be more localized, and the dynamics and coupling between them must be dealt with coherently. The Wannier functions from neighboring atomic sites might also need to be treated together. This is much like the LDA + U approach. For example, Cococcioni et al. 19 have shown that, in an LDA + U method, a linear response approach to calculate U (satisfying the Koopmans requirement) seems to fix the band gap problems in NiO. Such an approach can also change the occupation subspace, thus fixing the DFT groundstate problem.
It is worth to note that, not just WKM has a difficulty in describing such strongly correlated Mott insulator, other meanfield methods, including the GW method, have similar issues. As discussed in the Introduction, the GW results for MnO bandgaps varies from 2.3 13 to 4.4 eV 15 depending on the self-consistent procedure and initial input wave functions. The same is true for Fe 2 O 3 . For example, Liao and Carter 17 investigated the GW band gap for α-Fe 2 O 3 . They found that the calculated GW band gap can range from 1.3 to 4.8 eV depending on the initial input (from PBE, LDA + U to HSE06, PBE0), while the experimental band gap is 2.6 eV.
It is also an interesting general question whether the strongly correlated Mott insulator can be described by a Slater determinant in the wave function expression. Although the current WKM scheme fails to describe such system, this does not exclude the possibility that some variance of it, e.g., the one with Wannier function from both occupied and unoccupied orbitals might be successful. It remains to be seen whether one has to use multiconfiguration or local correlated treatment like in the dynamic mean-field theory to describe such systems.
In conclusion, we tested the WKM for the TM oxide systems. We found that it can predict well the band gaps for the d 0 and d 10 TM oxides, provided that a semicore removal procedure is used in the E hxc calculation. On the other hand, the current WKM procedure fails to predict the band gap for the partially occupied d-state TM oxides, even after tests of several procedures. The problem seems to have originated from the strong interaction between the valence band and conduction band d-state Wannier functions. To solve this problem, one has to go beyond the current WKM procedure. Instead, the d-state Wannier functions constructed from a combination of valence band and conduction band orbitals should be used, and their mutual interactions need to be treated coherently.

METHODS
All our calculations are performed using PWmat 41,42 , which runs on graphics processing unit. Fifty-Ryd plane wave cutoff was used in our calculations for SG15 pseudopotential. For FHI pseudopotential, the pseudopotential recommended cutoff energies for different elements are used accordingly. Monkhorst-Pack 43 method is used for k-points sampling. In order to ensure that the Wannier functions are orthogonal to semicore states, k-points setting for the bulk (for Wannier function generation) was consistent with our supercells (e.g., the bulk k-points are folded from the supercell Γ point). Generally speaking, the number of k-points multiplied by the number of atoms are >500 and <1000. Our supercell size is thus also between 500 and 1000 atoms, and only Γ point is used for supercell calculations. The calculated LDA band gaps between FHI and SG15 are typically very close as shown in Fig. 1. Wannier90 code 37 was used to generate Wannier functions. Initial guess of Wannier functions is needed