The method of density functional theory (DFT) combined with Hubbard U correction has been employed in wide-ranging applications in correlated systems. Although the DFT+U method can dramatically improve the treatment of many strongly correlated systems, careful attention must be paid to those parameters that significantly influence results. By varying the local-projection size with the muffin-tin radius, we examine the effect of the Hubbard U correction on lattice parameters, electronic structure, and the relative stability of TiO2 between rutile and anatase. Our results show that different ranges of the local projection can produce strikingly different quantitative results, leading to distinct conclusions. The quantitative discrepancies are attributed to the sensitivity of the occupancy numbers for the correlated orbitals toward the size of the projection subspace.

The method of density functional theory (DFT) combined with Hubbard U correction has been widely employed in many applications studying correlated systems.1,2 The DFT+U method typically combines the standard local density approximations (LDA) or generalized gradient approximations (GGA) with an explicit treatment of the correlation effect of the electrons in specific local orbitals only. Therefore, its popularity derives in part from the resulting high computational efficiency. The DFT+U method can dramatically improve the treatment of many strongly correlated systems. However, careful attention must be paid to the interpretation of the results, as a number of parameters can significantly influence the results.

Among the many input parameters used in the DFT+U, it has been well documented that the results of the DFT+U approach crucially depend on the values of the on-site Coulomb interaction represented by U and the exchange interaction J. While substantial efforts have been made to develop first-principle ways to compute the parameters for specific systems, it has become common practice to treat the values of U and J as adjustable parameters in search of agreement with particular experimental properties. Second, the complete formulation of DFT+U in rotationally invariant form describes fully orbital-dependent electronic interactions.3 In practice, many programs utilize a rather simplified functional that retains only the lowest order Slater integrals F0 in the calculation of electron–electron interactions.4 In this simplification, only UeffUJ parameter for a given l is needed. Although the simplified functional has yielded results similar to those of the fully rotationally invariant version for many materials, some studies have reported that the explicit inclusion of the Hund J was necessary to describe the ground state of strongly correlated systems properly.5,6

Another critical parameter lies within the definition of the local occupation matrix. In order to calculate the occupation matrix, a localized basis set has to be chosen with the defined correlated space, in which the Hubbard potential is applied. It is believed that the DFT+U calculations should not depend on the choice of the basis set as long as the localized basis sets are properly chosen and carefully tested. However, Wang et al.7 reported significant differences in the properties of MnO and MnO2 calculated from two popular pseudopotential-based codes, Quantum ESPRESSO (QE)8 and Vienna ab initio Simulation Package (VASP),9 as well as all-electron WIEN2k.10,11 The authors showed that the relative stability between ferromagnetic and antiferromagnetic phases, the bandgap, and the oxidation reaction energy crucially depend on the specific implementation. The observed discrepancies were attributed to the choice of local orbitals and their effective spatial ranges among the codes.

The aforementioned study has motivated us to examine the effect of spatial ranges of the local projection subspace in the DFT+U method systematically and to investigate the cause of dependence on the spatial extent. In this paper, we report our investigation of TiO2 rutile and anatase phases using all-electron WIEN2k. TiO2 is a quintessential transition metal oxide. It has attracted a long and sustained interest in the past due to its rich practical and potential applications in many areas, including heterogeneous and photo-catalysis.12–28 In addition, the relative stability of TiO2 among various phases has been of significant theoretical interest in the past.29 Although rutile is regarded as the most stable phase of bulk, as evidenced by experimental data,30 it has been well known that standard LDA and GGA, as well as hybrid functionals, could not reproduce the relative stability of TiO2 among various phases, particularly between anatase and rutile. More recently, a number of research groups have reported that the DFT+U approach allows the correct prediction of phase stability with a specific range of U.31–34 With a large number of theoretical and experimental studies carried out on these phases and readily available for comparison, TiO2 makes an ideal case to investigate the effect of the local-projection size in the DFT+U approach.

The structural models of rutile and anatase are briefly reviewed in Fig. 1. Rutile belongs to the tetragonal system (P42/mnm) whereas anatase to the body-centered tetragonal system (I41/amd).35 In both structures, Ti atom is surrounded by six oxygen atoms. In rutile, the approximate oxygen octahedra share edges that form a ribbon running along c direction. In anatase, edge shared oxygen octahedra run in zigzag along a and b directions. In both rutile and anatase, the oxygen octahedron is distorted with two distinct Ti–O bonds: Ti–O(1) from Ti to apical oxygen and Ti–O(2) from Ti to equatorial oxygen.14 It is worth noting that Ti–O(1) points along [110] or equivalent directions in rutile but along [001] in anatase.

FIG. 1.

Structures of TiO2 (a) rutile and (b) anatase are shown using VESTA.36 Two inequivalent Ti–O bonds and the coordinate axes for Ti d states are also shown.

FIG. 1.

Structures of TiO2 (a) rutile and (b) anatase are shown using VESTA.36 Two inequivalent Ti–O bonds and the coordinate axes for Ti d states are also shown.

Close modal

The paper is organized as follows. Section II describes the augmented plane waves plus the local orbital approach implemented in Wien2k and the simplified isotropic formulation of DFT+U. In Sec. III A, the results using a muffin-tin (MT) radius for Ti RMT(Ti) = 1.78 aB (Bohr radius) are compared to those from pseudopotential methods to demonstrate systematic, quantitative discrepancies. In Sec. III B, the size of the Hubbard U correlation sphere is varied via RMT(Ti) values. The results from different ranges of the Hubbard U interaction are presented and discussed. The main findings of the analysis are summarized in Sec. IV.

The spin-polarized calculations have been performed using the all-electron, augmented plane wave, and local orbital programs implemented in Wien2k.10,11 In the augmented plane wave scheme, the unit cell is divided into non-overlapping MT spheres, which are centered at atoms and interstitial regions. The linearized augmented plane wave (LAPW) approach utilizes a linear combination of radial functions multiplied by spherical harmonics within the atomic region,
ϕkn=lmAlm,knul(r,El)+Blm,knu̇l(r,El)Ylm(r̂),
(1)
where ul(r, El) is a solution of the radial Schrödinger equation for energy El, and u̇l(r,El) is the energy derivative of ul at the energy El. The coefficients Alm,kn and Blm,kn are then determined by matching the basis function to the plane wave of the interstitial region in value and slope at the muffin-tin radius,
ϕkn=n1Ωeiknr.
(2)
In a more efficient method of APW+lo, the standard APW basis is used with the usual fixed energy El along with a new local orbital (lo) for variational flexibility,
ϕkn=lmAlm,knul(r,El)Ylm(r̂),
(3)
ϕlmlo=lmAlmul(r,El)+Blmu̇l(r,El)Ylm(r̂),
(4)
where the coefficients Alm and Blm are now independent of kn but determined by the requirement that ϕlmlo vanishes at the sphere boundary and is normalized within the sphere. In all our calculations, APW+lo is used.

The core states are distinguished from valence states with an energy separation of −6.0 Ry, resulting in Ti 3s and 3p as semi-core states, whereas Ti 3d, 4s, O 2s, and 2p are valence states. The core electrons are treated in a fully relativistic manner using the Dirac–Fock method. The valence electrons are treated with a scalar relativistic approach. In treating the exchange-correlation effect at the GGA-level of approximation, the Perdew–Burke–Ernzerhof (PBE) functional was used.37,38

The muffin-tin radii for Ti and O were initially chosen as 1.78 and 1.61 aB, respectively, to represent a reference system with a small correlation space. The accuracy of the calculation is controlled by the dimensionless cutoff parameter RMTKmax, which determines the size of the basis through the product of the smallest RMT value and the largest K vector in the plane wave expansion. The value of RMTKmax was increased and tested up to 10. During self-consistency cycles, energy converged to 0.1 mRy. The internal coordinates were optimized by minimizing the forces to less than 1 mRy/bohr. Because the energy difference between the TiO2 phases can be small, the convergence of the calculated energy with respect to various parameters was carefully checked. The Brillouin Zone was sampled from 700 to 2000 points for convergence.

In order to investigate the effect of varying the spatial range of the Hubbard interaction on the DFT+U calculations, the RMT(Ti) values were systematically increased from 1.78 to 2.30 aB. Because the L/APW method requires non-overlapping spheres for atomic regions, RMT values for oxygen were reduced accordingly from 1.61 to 1.30 aB. For larger RMT (Ti) values, the pronounced dissimilarity in muffin-tin radii between Ti and O could be problematic, leading to very different effective RMTKmax values. To avoid the linearization of the basis, RMTKmax was thus reduced to 7.5, and a high energy local orbital was included for the largest RMT(Ti) = 2.3 aB. In addition, the increased values for the spherical expansion in the MT spheres, the potential, and the charge density expansion were tested and used to achieve the good convergence and RMT-independence of the GGA calculations.

The electron correlation effect of d electrons was taken into account using the orbital dependent potential formulated with the Hubbard on-site interaction term.39 The interaction is based on the rotationally or unitary-transformation invariant formulation,3,
EHub{nmm}=12{m},σm,m|Vee|m,mnmmσnmmσ+m,m|Vee|m,mm,m|Vee|m,mnmmσnmmσ,
(5)
where nmmσ is the density matrix for correlated electrons with the magnetic quantum number of the orbital angular momentum l and spin index σ, and Vee are the screened Coulomb interactions. Because the standard density functional approximations such as LDA and GGA already include the part of the interaction energy that the Hubbard model aims to better describe, it must be subtracted via the double-counting term Edc. The most widely used method of eliminating double counting for strongly correlated insulating systems was applied39 
Edc{nmm}=U2n(n1)J2n(n1)n(n1),
(6)
where nσ=Tr(nmmσ) and n = n + n are the occupation numbers of localized orbitals.
The interaction energy functional can be simplified by keeping only the lowest order Slater integral F0 [4,2] in evaluating the interaction matrix elements ⟨m, m″|Vee|m′, m‴⟩. Under this isotropic scheme of the DFT+U, the Hubbard correction with the double counting term included is written as
EHub=σUeff2Trnσ(1nσ),
(7)
with Ueff regarded as an effective parameter, including the exchange J such that Ueff = UJ. Consequently, by differentiating with respect to each orbital occupancy, one obtains the additional orbital-dependent potential energy
Vmσ=Ueff12nmσ.
(8)

In this study, the orbital dependent correction by Eq. (8) was applied to the localized orbitals described by the basis function [Eqs. (3) and (4)] within the augmentation sphere defined by RMT. Throughout the paper, Ueff is referred to simply as U.

Structural optimization of rutile and anatase phases has been performed with fully relaxing internal parameters using the PBE functional, with the Hubbard U value increasing from 0 up to 10 eV (Table I). The calculated lattice constants for rutile are in good agreement with the experimental values of a = 4.593 Å and c = 2.959 Å,35 only slightly overestimated by a maximum of 1.2%. The lattice constants for anatase are also overestimated by about 2%, compared to the experimental values of a = 3.785 Å and c = 9.512 Å.35 As the U value increases, lattice constants for both rutile and anatase gradually expand, making the theoretical values more unsatisfactory in comparison with the GGA as well as the experimental ones. More specifically, the lattice constant a for rutile increases to only 0.2% with U = 10 eV from the GGA lattice value. On the other hand, the lattice constant c increases to 1.9% with U = 10 eV, compared to that with GGA only. For anatase, the lattice constant a exhibits a larger expansion of about 1.3% with U varying from 0 to 10 eV, whereas the lattice constant c shows little or no increase.

TABLE I.

Calculated lattice parameters of rutile and anatase as a function of U. The number in parenthesis is the percentage change from the GGA value (U = 0 eV). Experimental values are denoted by italics.

Lattice constants
PhaseMethodU (eV)a (Å)c (Å)c/a
Rutile GGA 4.647 2.966 0.638 
 GGA+U 4.650 (+0.06) 2.985 (+0.64) 0.642 
 GGA+U 4.649 (+0.04) 2.995 (+0.98) 0.644 
 GGA+U 4.652 (+0.11) 3.011 (+1.52) 0.647 
 GGA+U 10 4.655 (+0.17) 3.021 (+1.85) 0.649 
 Expt.35   4.593 2.959 0.644 
Anatase GGA 3.801 9.705 2.553 
 GGA+U 3.819 (+0.47) 9.693 (−0.12) 2.538 
 GGA+U 3.828 (+0.71) 9.695 (−0.10) 2.532 
 GGA+U 3.843 (+1.10) 9.705 (−0.004) 2.526 
 GGA+U 10 3.851 (+1.31) 9.715 (+0.10) 2.522 
 Expt.35   3.785 9.512 2.513 
Lattice constants
PhaseMethodU (eV)a (Å)c (Å)c/a
Rutile GGA 4.647 2.966 0.638 
 GGA+U 4.650 (+0.06) 2.985 (+0.64) 0.642 
 GGA+U 4.649 (+0.04) 2.995 (+0.98) 0.644 
 GGA+U 4.652 (+0.11) 3.011 (+1.52) 0.647 
 GGA+U 10 4.655 (+0.17) 3.021 (+1.85) 0.649 
 Expt.35   4.593 2.959 0.644 
Anatase GGA 3.801 9.705 2.553 
 GGA+U 3.819 (+0.47) 9.693 (−0.12) 2.538 
 GGA+U 3.828 (+0.71) 9.695 (−0.10) 2.532 
 GGA+U 3.843 (+1.10) 9.705 (−0.004) 2.526 
 GGA+U 10 3.851 (+1.31) 9.715 (+0.10) 2.522 
 Expt.35   3.785 9.512 2.513 

The calculated lattice parameters for rutile are further compared to those reported in the literature using pseudopotential and projector augmented wave (PAW) methods [Fig. 2(a)]. The lattice values among those selected references spread more or less by about 1%. Our values (a = 4.647 Å and c = 2.966 Å) are well within the range of reported values using the same PBE functional. However, with the Hubbard U interaction included, the pseudopotential and PAW calculations appear to exhibit a far greater dependence of the lattice constants a and c on the Hubbard U values, in sharp contrast to the current study (blue circles). For example, using the PAW method in VASP, Arroyo-de Dompablo et al.34 reported that the lattice constants a and c expand from 4.650 to 4.725 Å and from 2.968 to 3.108 Å, respectively (orange squares). These increases represent a 1.6% and 4.7% expansion, respectively, as U increases from 0 to 10 eV.

FIG. 2.

(a) Comparison of the calculated lattice constants for rutile with those using pseudopotential and PAW methods as a function of U and (b) the percentage changes in the lattice constants from their respective GGA values. The dashed lines are only to guide the eyes. References a: 34, b: 40, c: 31, d: 33, e: 41.

FIG. 2.

(a) Comparison of the calculated lattice constants for rutile with those using pseudopotential and PAW methods as a function of U and (b) the percentage changes in the lattice constants from their respective GGA values. The dashed lines are only to guide the eyes. References a: 34, b: 40, c: 31, d: 33, e: 41.

Close modal

In order to make the comparison among various studies easier on the effect of U, the reported lattice constants are normalized to their GGA values. The percentage changes in the lattice constants are plotted in [Fig. 2(b)]. It is clear that the lattice constants reported in this study are much less sensitive to the U values compared to those from other studies using pseudopotential and PAW methods. With U increasing up to 10 eV, the other studies show that the lattice constant a expands as much as about 1.6% and the lattice constant c expands to a far greater extent to 6%. The rather large discrepancy in the predicted values for the lattice constants is disturbing, but it is also noteworthy that all the studies paint the same qualitative trend. The lattice constants generally increase with increasing U, and the lattice constant c is apparently more sensitive to the U values than a.

The comparison of the calculated lattice constants for anatase shows a similar quantitative discrepancy between the present and previously reported results (Fig. 3). The pseudopotential and PAW methods apparently predict a much greater dependence of the lattice constants on the Hubbard U values than this study predicts. The lattice constants a and c increase by up to 3.7% and 1.8%, respectively, for U = 10 eV,34 in contrast to about 1.2% for a and little or no change in c predicted in the current study. Again, despite the clear deviation in the quantitative prediction, both this work and previously reported studies show the same qualitative trend of lattice constants expanding with U. For anatase, it is the lattice constant a expanding by a greater amount with the Hubbard U value.

FIG. 3.

(a) Comparison of the calculated lattice constants for anatase with those using pseudopotential and PAW methods as a function of U and (b) the percentage changes in the lattice constants from their respective GGA values. The dashed lines are only to guide the eyes. References a: 34, c: 31, d: 33, e: 41.

FIG. 3.

(a) Comparison of the calculated lattice constants for anatase with those using pseudopotential and PAW methods as a function of U and (b) the percentage changes in the lattice constants from their respective GGA values. The dashed lines are only to guide the eyes. References a: 34, c: 31, d: 33, e: 41.

Close modal

The electronic structures of rutile and anatase are also examined and compared to the results from pseudopotential and PAW methods. The density of states (DOS) of rutile and anatase are presented in Fig. 4. The calculated bandgaps with the PBE functional are 1.93 and 2.17 eV for rutile and anatase, respectively. They are similar to the previously reported values of the DFT studies34,42–44 but well short of the experimental values of 3.03 eV45 and 3.2 eV,46 respectively. The severe underestimation of bandgaps is attributed to self-interaction error (SIE), known as the failure of standard density functional approximation.47 

FIG. 4.

The projected density of states for rutile (left) and anatase (right) with the selected U values: Ti total d (blue), dz2 (black), dx2y2 (gray), dxy (purple), dyz (yellow), dzx (green), and O p (red). The Fermi level is arbitrary and set to the top of the valence band maximum. In addition, the bandgap energy and the crystal field splitting between t2g and eg are indicated.

FIG. 4.

The projected density of states for rutile (left) and anatase (right) with the selected U values: Ti total d (blue), dz2 (black), dx2y2 (gray), dxy (purple), dyz (yellow), dzx (green), and O p (red). The Fermi level is arbitrary and set to the top of the valence band maximum. In addition, the bandgap energy and the crystal field splitting between t2g and eg are indicated.

Close modal

The valence band is mainly made of O 2p states. The empty conduction band consists of predominantly Ti 3d states, thus indicating a strong charge-transfer insulator.48,49 For rutile, the conduction band is split into two sub-bands: the lower one between 1.9 and 4.4 eV and the upper one between 4.4 and higher. The lower subband is identified as the Ti t2g band and the upper one as the Ti eg band in approximating the site symmetry of Ti as the octahedral symmetry Oh.50–53 Further decomposition of the Ti d projected DOS (PDOS) using the local coordinate axes (Fig. 1) reveals that the lower conduction band is mainly composed of dxy, dyz, and dzx characters. Likewise, the upper conduction band is predominantly composed of dz2 and dx2y2 characters, indeed justifying the assignment of the Oh symmetry. The calculated crystal field splitting (CFS) of 2.3 eV between the t2g and eg bands is in fair agreement with the experimental value of 2.6 eV.51 

It is notable that the Ti 3d states and the O 2p states hybridize substantially in the middle and lower valence band regions.50,54 The overlap between the Ti eg states and O 2p states around −4 eV mainly represents the σ bonding in the molecular orbital (MO) scheme. In the higher energy region of the valence band, Ti dyz and dzx form weak π-bonding with O 2p states, and at the top of the valence band, non-bonding Ti dxy is observed in Figs. S1 and S2 in the supplementary material.50 The anatase PDOS shows a qualitatively similar picture of a strong charge-transfer insulator with a significant σ-type hybridization between Ti eg states and O 2p states (Fig. 4 right), which is consistent with early studies.54 

As the Hubbard U interaction is included and eventually increases to 10 eV, there are a couple of prominent changes taking place in the electronic structure. The bandgaps for both rutile and anatase widen gradually and reach 2.98 and 3.18 eV at U = 10 eV, closer to the experimental values (Fig. S3 in the supplementary material). The opening or widening of a bandgap has been a hallmark of the Hubbard U correction applied to the conventional DFT due to alleviated SIE. Our results are generally in good agreement with the theoretical values previously reported in the literature.34,42,43 Another interesting change is the decrease in the CFS between the t2g and eg states in the conduction band. The CFS values decrease from 2.3 to 2.1 eV for rutile and anatase, respectively, to 1.8 and 1.5 eV at U = 10 eV (Fig. 4). The narrowing of the CFS with increasing U suggests that the t2g and eg sub-bands are affected disproportionately by the Hubbard U interaction term. For example, the anti-bonding π* t2g band is pushed up to a greater extent relative to the σ* eg band.

The previous calculations show again a qualitatively similar trend of the CFS but quantitatively distinct results upon increasing the Hubbard U interaction. The calculation using the PAW method34 started out with a similar CFS of 2.4 eV for rutile using the PBE functional but reported its decrease to 0.9 eV at U = 10 eV, in sharp contrast to the 1.8 eV shown in Fig. 4. Using the pseudopotential method in QE, Zhao and Kulik40 reported that the CFS between the two d sub-bands decreases even more rapidly as the U value increases. The t2g and eg peaks are completely overlapped by the time U reaches 8 eV, whereas their splitting in our study decreases much more gradually and still remains well-resolved to be 1.9 eV (Fig. S4 in the supplementary material). Qualitatively similar but quantitatively diverging results in both DOS and the optimized structural parameters hint at the fact that systematic inconsistencies may play a role in the Hubbard U implementation among the methods.

The effect of the Hubbard U interaction on DFT calculations can be further appreciated by studying the difference electron density plot, which shows the difference in electron density between the crystal and the superposed atoms. In Fig. 5, the difference electron density plots are compared before and after applying the Hubbard interaction with U = 10 eV. Without the Hubbard U interaction, the plot of difference electron density shows the large deficiency and excess of the electronic charges around Ti and O, respectively, consistent with a strong charge transfer nature. The angular profile of charge deficiency around Ti reflects the t2g-type orbital symmetry.50 For oxygen, three distinct regions of electron surplus along the three nearby Ti atoms represent the sp2-type bonding for which oxygen forms a planar O–Ti3 cluster. One also observes excess electronic charges at Ti directed along the octahedrally coordinated O atoms, corresponding to the eg orbitals observed in PDOS. The significant charge accumulation between Ti and O along the bond directions clearly supports substantial covalency through eg and sp hybridization. Furthermore, a closer examination reveals that the accumulated charges between Ti and the equatorial O are greater than those between Ti and the apical O. The observed charge density difference is in good agreement with previously reported experimental and theoretical charge-density maps.50,53

FIG. 5.

The plots of the difference electron density for rutile in the (110) plane using (a) PBE and (b) PBE+U (10 eV). The blue indicates the electronic charge deficiency, and the red indicates the surplus.

FIG. 5.

The plots of the difference electron density for rutile in the (110) plane using (a) PBE and (b) PBE+U (10 eV). The blue indicates the electronic charge deficiency, and the red indicates the surplus.

Close modal

With the inclusion of U = 10 eV, both areas showing the charge deficiency around Ti (blue) and the excess around O (red) are significantly enlarged [Fig. 5(b)]. The enhancement of charge depletion around Ti and charge excess around O indicates an increased charge transfer between Ti and O brought about by the Hubbard U interaction. At the same time, the small charge surplus seen at Ti along the Ti–O bond directions (eg symmetry) is strongly suppressed, suggesting the curtailment of the d-sp hybridization.

The partial charges of Ti are calculated by integrating the PDOS (Fig. 4) for a quantitative assessment of the charges in the octahedral symmetry. Table II (top) lists the occupancy of each d sub-orbital per spin, projected within the muffin-tin sphere with a radius of 1.78 aB. The occupancy in the eg orbitals is much higher than that in the t2g orbitals, consistent with the dominant σ bonding observed in Fig. 5. As U increases eventually up to 10 eV, the occupancy of each d sub-orbital decreases as expected from the effect of the Hubbard potential on the less than half-filled orbital. The total d occupancy decreases by 0.12e per spin or 0.24e for both spins, confirming enhanced charge depletion at Ti. A very similar effect of Hubbard U on anatase is also observed (Table S1 in the supplementary material).

TABLE II.

The occupancies of d orbitals (per spin) of rutile as a function of U at two different RMT (Ti) values.

U(RMT(Ti) = 1.78 aB)Total ddz2dx2y2dxydyzdzx
0.516 0.139 0.148 0.070 0.078 0.079 
0.481 0.134 0.145 0.062 0.070 0.071 
0.459 0.130 0.140 0.057 0.065 0.066 
0.424 0.123 0.135 0.050 0.057 0.058 
10 0.394 0.118 0.125 0.045 0.052 0.053 
U(RMT(Ti) = 2.3 aBTotal d dz2 dx2y2 dxy dyz dzx 
0.759 0.207 0.223 0.104 0.112 0.112 
0.709 0.205 0.217 0.089 0.098 0.099 
0.680 0.202 0.219 0.080 0.089 0.089 
0.617 0.194 0.204 0.067 0.076 0.076 
10 0.580 0.190 0.196 0.060 0.068 0.067 
U(RMT(Ti) = 1.78 aB)Total ddz2dx2y2dxydyzdzx
0.516 0.139 0.148 0.070 0.078 0.079 
0.481 0.134 0.145 0.062 0.070 0.071 
0.459 0.130 0.140 0.057 0.065 0.066 
0.424 0.123 0.135 0.050 0.057 0.058 
10 0.394 0.118 0.125 0.045 0.052 0.053 
U(RMT(Ti) = 2.3 aBTotal d dz2 dx2y2 dxy dyz dzx 
0.759 0.207 0.223 0.104 0.112 0.112 
0.709 0.205 0.217 0.089 0.098 0.099 
0.680 0.202 0.219 0.080 0.089 0.089 
0.617 0.194 0.204 0.067 0.076 0.076 
10 0.580 0.190 0.196 0.060 0.068 0.067 

Because the calculated partial charges depend on the muffin-tin radius, the Bader charge analysis has also been performed.55,56 Table III shows that the Bader charges for Ti and O are +2.28e and −1.14e, respectively, using the PBE functional. These values are in good agreement with previously calculated values, for instance, +2.24e for Ti.40 As U increases to 10 eV, the Bader charges increase to +2.58e and −1.29e, respectively, indicating that Ti loses an additional 0.3e to two O atoms. This value is also comparable to about 0.25e calculated using the pseudopotential method.40 

TABLE III.

Bader charges and bond critical points are provided for U = 0 and 10 eV, calculated with two muffin-tin radii of Ti. For each BCP along Ti–O(1) (apical O) and Ti–O(2) (equatorial O), the two distances from Ti and O are listed.

RutileBader charge (e)Critical points (3, −1) from (in Å)
RMT(Ti) (aB)U (eV)TiOTi; O(1)Ti; O(2)
1.78 +2.28 −1.14 1.002; 1.004 0.982; 0.977 
1.78 10 +2.58 −1.29 0.992; 1.016 0.981; 1.001 
2.3 +2.28 −1.14 1.008; 0.995 0.991; 0.966 
2.3 10 +2.60 −1.30 1.012; 1.021 1.005; 1.008 
Anatase  Bader charge (eCritical points (3, −1) from (in Å) 
1.78 +2.26 −1.13 0.997; 0.998 0.975; 0.974 
1.78 10 +2.56 −1.28 0.992; 1.017 0.974; 0.997 
2.3 +2.26 −1.13 1.004; 0.993 0.983; 0.962 
2.3 10 +2.58 −1.29 1.012; 1.017 0.997; 1.002 
RutileBader charge (e)Critical points (3, −1) from (in Å)
RMT(Ti) (aB)U (eV)TiOTi; O(1)Ti; O(2)
1.78 +2.28 −1.14 1.002; 1.004 0.982; 0.977 
1.78 10 +2.58 −1.29 0.992; 1.016 0.981; 1.001 
2.3 +2.28 −1.14 1.008; 0.995 0.991; 0.966 
2.3 10 +2.60 −1.30 1.012; 1.021 1.005; 1.008 
Anatase  Bader charge (eCritical points (3, −1) from (in Å) 
1.78 +2.26 −1.13 0.997; 0.998 0.975; 0.974 
1.78 10 +2.56 −1.28 0.992; 1.017 0.974; 0.997 
2.3 +2.26 −1.13 1.004; 0.993 0.983; 0.962 
2.3 10 +2.58 −1.29 1.012; 1.017 0.997; 1.002 

The amount of charge transfer from the Bader analysis is greater than 0.24e from the partial charge analysis (Table II). The difference is likely due to the fact that the Ti atom estimated in the Bader analysis is bigger than the Ti atomic sphere in the APW method using RMT(Ti) = 1.78 aB (0.94 Å). To get an effective size of Ti in the Bader analysis, a bond critical point (BCP) (3, −1), where the zero flux surface intercepts, is located. Along the bond between Ti and apical O(1), a BCP is located at 1.002 Å from Ti and 1.004 Å from O(1) before the Hubbard U interaction turned on (Table III). Another BCP is found at 0.982 Å from Ti and 0.977 Å from O(2) along the bond between Ti and equatorial O. These BCPs are well outside the Ti augmentation sphere used in the APW. With the inclusion of U = 10 eV, the BCP moves farther from O at 1.016 Å along the bond direction Ti–O(1) and slightly closer to 0.992 Å from Ti. Similarly, along the bond direction Ti–O(2), the BCP is found substantially farther at 1.001 Å from O(2) but approximately at the same distance of 0.981 Å from Ti. The substantially bigger oxygen with a slightly smaller Ti is congruous with a more ionic TiO2 resulting from the Hubbard U interaction.

Besides the overall effect of enhanced charge transfer, a closer inspection of the sub-orbitals in the partial charge analysis (Table II) reveals a significant variation in the occupancies caused by the Hubbard U. By the time U increases to 10 eV, the occupancies in the dz2 and dx2y2 orbitals decrease to 0.118 and 0.125, respectively, representing 84%–85% of the U = 0 eV values. On the other hand, the occupancies in the dxy, dyz, and dzx orbitals decline much sharply to 0.045, 0.052, and 0.053, respectively, only 64%–67% of the initial occupancies. The uneven reduction of charge population between the eg and t2g orbitals is reflected in the density difference map (Fig. 5), which shows that the depleted region of t2g symmetry (blue) has enlarged to a greater extent than the region of eg symmetry (red) has depleted. A similar effect of U on the density difference map for anatase is also observed (Fig. S5 in the supplementary material).

The charge redistribution and change in the bonding nature between Ti and O can be used to understand the accompanying structural changes that the Hubbard U brings about. Table IV lists the two distinctive Ti–O bond lengths: Ti–O(1) and Ti–O(2) between Ti and apical and equatorial oxygen atoms, respectively (Fig. 1), as a function of U. For both rutile and anatase, the long Ti–O(1) bond shows little or no changes, whereas the short Ti–O(2) bond in the equatorial plane exhibits much greater sensitivity to the U values. The lack of significant changes in the Ti–O(1) bond mirrors more or less constant values of a for rutile and c for anatase observed as a function of U in Figs. 2 and 3. On the other hand, a higher sensitivity of the Ti–O(2) bond on U mimics more remarkable changes in c for rutile and a for anatase, as the short Ti–O(2) bonds in edge-shared octahedral ribbons run along c for rutile and a and b axes for anatase, respectively.

TABLE IV.

The two inequivalent Ti–O bond lengths are listed as a function of U (eV) calculated using RMT = 1.78 aB for rutile and anatase.

RutileAnatase
MethodU (eV)Ti–O(1) (Å)Ti–O(2) (Å)Ti–O(1) (Å)Ti–O(2) (Å)
GGA 2.006 1.958 1.995 1.948 
GGA+U 2.007 1.966 2.002 1.955 
GGA+U 2.006 1.970 2.003 1.959 
GGA+U 2.007 1.977 2.006 1.966 
GGA+U 10 2.008 1.981 2.008 1.970 
RutileAnatase
MethodU (eV)Ti–O(1) (Å)Ti–O(2) (Å)Ti–O(1) (Å)Ti–O(2) (Å)
GGA 2.006 1.958 1.995 1.948 
GGA+U 2.007 1.966 2.002 1.955 
GGA+U 2.006 1.970 2.003 1.959 
GGA+U 2.007 1.977 2.006 1.966 
GGA+U 10 2.008 1.981 2.008 1.970 

Yoshiya et al.57 noted in their investigation of TiO2, Ti2O3, and TiO that the covalency decreases as the formal charge decreases from +4 to +2. The reduction in covalency was attributed to more d electrons occupying the anti-bonding σ∗ orbitals in the conduction band. In corundum Ti2O3, the two distinctive Ti–O bonds are 2.07 and 2.02 Å.58,59 For the cubic TiO, the Ti–O bond is 2.09 Å long.60 In the present study of rutile and anatase TiO2, the shorter Ti–O bond lengths correspond to substantial covalency, as noted earlier. Therefore, the elongation of the Ti–O bond lengths can be explained in terms of the increased ionicity with the diminished covalent characteristics upon the application of the Hubbard interaction. The change in bonding nature is facilitated by on-site Coulomb repulsion to depopulate the bonding eg as well as non-bonding t2g orbitals. A more pronounced expansion in the equatorial Ti–O(2) results from its much shorter initial bond-length in comparison with Ti–O(1). The shorter Ti–O(2) bond length may be viewed as a stronger covalent character through dx2y2 predicted by GGA (Table II).

The comparison of the structural parameters and the electronic structures of TiO2 between the present and previously reported studies shows a similar qualitative trend for the effect of the Hubbard U. Yet, the quantitative discrepancy can become disturbingly large. In order to directly investigate the effect of the spatial range of the local projection in the DFT+U calculations, the augmentation sphere of Ti is systematically varied using the following set of muffin-tin radii in units of aB: RMT(Ti) = 1.78, RMT(O) = 1.61; RMT(Ti) = 1.91, RMT(O) = 1.73; RMT(Ti) = 2.1, RMT(O) = 1.5; RMT(Ti) = 2.3, RMT(O) = 1.3.

Figure 6 exhibits the changes in the lattice constants of rutile and anatase with an increasing muffin-tin radius of Ti using U= 10 eV. The lattice constant a for rutile displays a dramatic dependence on the muffin-tin radius of Ti [Fig. 6(a)], in sharp contrast to its little or no dependence on the U value using RMT(Ti) = 1.78 (Fig. 2). It increases from 4.655 Å with RMT(Ti) = 1.78 aB to 4.696 Å with RMT(Ti) = 2.3 aB. The lattice constant c also shows its increasing value with a larger muffin-tin radius of Ti, from 3.021 Å with RMT(Ti) = 1.78 aB to 3.077 Å with RMT(Ti) = 2.3 aB [Fig. 6(b)]. The increased lattice constants become more comparable to the values reported with PAW and pseudopotential methods in Fig. 2. For example, the lattice constants of 4.696 and 3.077 Å for a and c, respectively, using RMT(Ti) = 2.3 aB are much closer to 4.725 and 3.108 Å for a and c, respectively, for U = 10 eV reported using the PAW method in Ref. 34. Likewise, the lattice constants of anatase exhibit their pronounced dependence on the muffin-tin radius of Ti. For U = 10 eV, the values of a and c increase from 3.849 and 9.719 Å, respectively, to 3.908 and 9.814 Å as the muffin-tin radius of Ti is increased from RMT(Ti) = 1.78 aB to RMT(Ti) = 2.3 aB [Figs. 6(c) and 6(d)]. These values are again much better matched with the values obtained from PAW and pseudopotential methods (Fig. 3), for example, 3.947 and 9.871 Å with U = 10 eV.34 

FIG. 6.

Lattice constants (a) a and (b) c for rutile and (c) and (d) similarly for anatase plotted against various muffin-tin radii of Ti using U = 10 eV. The blue dashed line is to only guide eyes.

FIG. 6.

Lattice constants (a) a and (b) c for rutile and (c) and (d) similarly for anatase plotted against various muffin-tin radii of Ti using U = 10 eV. The blue dashed line is to only guide eyes.

Close modal

The lattice parameters obtained using the largest muffin-tin radius of 2.3 aB for Ti are still a bit smaller compared to those from the PAW method.34 In the PAW formalism, the Hubbard U interaction is implemented using the projection function and the pseudo-wavefunctions inside the augmentation sphere. While the PAW augmentation radius for Ti can vary depending on the inclusion of semi-core electrons, among other factors, the widely used potential Ti pv has RPAW = 2.5 aB.61 Therefore, we ascribe the remaining discrepancy in the lattice parameters to the fact that the muffin-tin sphere used in this study is still likely smaller than the PAW augmentation sphere used in the previous studies. For a better agreement, one could calculate RMT(Ti) = 2.5 aB. However, a bigger Ti muffin-tin sphere would necessitate further reduction of RMT(O) in the APW scheme, risking approximate linear dependencies of the basis due to very different effective RMTKmax. We believe that RMT(Ti) = 2.3 aB is about the maximum size for trustworthy results. Nevertheless, the unmistakable trend of the lattice constants with an increasing RMT(Ti) value in Fig. 6 demonstrates the crucial role that the size of the Hubbard U correlation space can play in quantitative predictions using DFT+U approaches.

The effect of the size of the correlation space on the relative energetics between rutile and anatase is also investigated by comparing the calculations using RMT(Ti) = 1.78 aB and RMT(Ti) = 2.3 aB (Fig. 7). With the PBE functional and the small RMT value for Ti, the anatase is more stable than rutile by −0.098 eV (blue bar), which is comparable to those reported in the literature.31,34 With an increasing U value, the energy difference between rutile and anatase gradually diminishes to −0.031 eV at U = 10 eV. However, anatase still remains more stable than rutile, in clear contrast to the results from the pseudopotential and PAW methods. For instance, Arroyo-de Dompablo et al.34 reported that the total energy of rutile becomes lower than that of anatase with any U value larger than 5 eV in the GGA+U approach. The authors further suggested that with the range of U values between 5 and 8 eV, their study predicts the correct ordering in energetics among the phases, including columbite: Erutile < Eanatase < Ecolumbite. Similarly, Curnan and Kitchin reported that the U values between 4.7 and 7.0 eV with the PBE functional reproduce the ordering of the stability among rutile, anatase, brookite, and columbite phases consistent with experiments: Erutile < Eanatase < Ebrookite < Ecolumbite.32 

FIG. 7.

The energy of anatase relative to that of rutile (per formula unit) as a function of U values, calculated with three sets of muffin-tin radii for Ti: 1.78, O: 1.61 (blue), Ti: 1.91, O: 1.73 (black), and Ti: 2.30, O: 1.30 (red) aB. The experimental value is from Ref. 30.

FIG. 7.

The energy of anatase relative to that of rutile (per formula unit) as a function of U values, calculated with three sets of muffin-tin radii for Ti: 1.78, O: 1.61 (blue), Ti: 1.91, O: 1.73 (black), and Ti: 2.30, O: 1.30 (red) aB. The experimental value is from Ref. 30.

Close modal

Using RMT(Ti) = 2.3 aB (red bars in Fig. 7), the calculation exhibits that the energy difference between rutile and anatase decreases more rapidly with increasing U values. The relative stability is now reversed between U = 5 and 8 eV, making the rutile more stable, in better agreement with the aforementioned results. At U = 10 eV, the energy difference between rutile and anatase becomes about 53 meV. The dependency of the relative energetics between rutile and anatase, together with their lattice parameters, on RMT is unsettling. It underscores the imperative role of the correlation-space size, which can lead to a different conclusion using the DFT+U approaches.

The sensitivity of the DFT+U results on the definition of the local projection space was hinted at early.62–64 However, a growing number of studies in recent years have carried out more focused studies on the different local projection methods and their spatial extents affecting the DFT+U results in various materials and implementations. Wang and co-workers61 carried out a PBE+U study on the polaronic properties of energy materials using various PAW augmentation radii and found that the size of the project space can influence many calculated properties substantially. For instance, for rutile TiO2, the polaronic gap state energy, the formation energy, and the polaronic hopping barrier all showed a strong dependence, as RPAW for Ti varied from 1.9 to 2.8 aB. The authors attributed the observed connection to the Hubbard U term, which depends quadratically on the local project density matrix. G. Geneste et al.65 applied the PBE+U method to the oxygen p orbitals to investigate the properties of the oxygen-type hole polarons in doped barium stannate. They reported that the Hubbard U value required to localize the hole as a small polaron depends significantly on whether the density matrix of the truncated atomic orbital is renormalized or not. Recently, Outerovitch and Amadon66 compared the truncated/renormalized atomic orbitals with the projected localized Wannier functions as the local orbitals in the study of UO2 and TiO2. Using the U values obtained from the constrained random phase approximation (cRPA), they reported that the equilibrium volume can either expand (as previously reported67) or contract upon the addition of the oxygen U potential, whether the correlated orbitals are atomic or Wannier. Similarly, Park et al.68 used three different correlated subspaces constructed by un-normalized, orthonormalized, and maximally localized Wannier function (MLWF) orbitals to study the properties of rare earth nickelates. They found that the use of the un-normalized PAW method leads to especially different results compared to the results from the orthonormalized and MLWF projectors.

As for the first step to gain insight into the effect of RMT(Ti) on our PBE+U results, the difference electron density map calculated with U = 8 eV using RMT(Ti) = 2.30 aB is stacked on top of that using RMT(Ti) = 1.78 aB in Fig. 8. The result using the large RMT(Ti) value exhibits a similar pattern of overall charge transfer from Ti to O as that using the small RMT(Ti). However it is clearly noticeable that the calculation with large RMT(Ti) predicts significantly more charges remaining near the Ti atom in the eg orbitals along the bonding directions to O after the Hubbard U application. For a quantitative analysis, various d orbital populations calculated using the two different muffin-tin radii in Table II are also compared.

FIG. 8.

The plot of the difference electron density for rutile in the (110) plane with U = 8 eV using (a) RMT(Ti) = 2.3 aB and (b) RMT(Ti) = 1.78 aB. The blue indicates the electronic charge deficiency and the red surplus, compared to the superposed atomic densities.

FIG. 8.

The plot of the difference electron density for rutile in the (110) plane with U = 8 eV using (a) RMT(Ti) = 2.3 aB and (b) RMT(Ti) = 1.78 aB. The blue indicates the electronic charge deficiency and the red surplus, compared to the superposed atomic densities.

Close modal

The calculation using RMT(Ti) = 2.30 aB shows that the occupancy in all d orbitals decreases with increasing U, as observed with RMT(Ti) = 1.78 aB. However, a detailed examination reveals that the occupancy in the σ-bonding eg orbitals (dz2 and dx2y2) reduces to only 88%–92% of the initial values at U = 10 eV as compared to 84%–85% of the initial values using the small RMT value. On the other hand, the occupancy in the π-type t2g orbitals (dxy, dyz, dzx) decreases to 58%–61% of the initial values, more than 64%–67% of the initial values using the small RMT = 1.78 aB. Therefore, the calculation with the large RMT(Ti) produces a greater variation in occupancy between the eg and t2g orbitals upon the U term. The large variation is also seen in the charge density difference map, where the electron depletion in t2g orbitals reaches further into the interstitial regions between two oxygens in the calculation with the large RMT(Ti) (Fig. 8).

The larger disparity in occupancy between the eg and t2g orbitals can also explain a bigger reduction in the CFS predicted by the calculation using the large RMT(Ti) value. It is evident from Eq. (8) that the larger the occupancies in eg orbitals, the smaller the Hubbard potential shifts Veg are. On the contrary, the potential shifts Vt2g of the t2g orbitals are affected to a lesser extent as they are maxed out for being nearly empty. Consequently, the CFS between the two bands should decrease more rapidly for RMT(Ti) = 2.3 aB, compared to the calculation using RMT(Ti) = 1.78 aB. Indeed, this is qualitatively what is observed in the PDOS of rutile with U = 8 eV when compared using the two different muffin-tin radii (Fig. 9). It is also noteworthy to observe that the bandgap is largely determined by the separation between the t2g and the O 2p orbitals. Therefore, the bandgap is less sensitive to the size of the correlation sphere for Ti than the crystal field splitting.

FIG. 9.

The projected density of states for rutile with U = 8 eV using RMT(Ti) = 2.3 aB (top) and RMT(Ti) = 1.78 aB (bottom): Ti d (blue), dz2 (black), dx2y2 (gray), dxy (purple), dyz (yellow), dzx (green), and O p (red).

FIG. 9.

The projected density of states for rutile with U = 8 eV using RMT(Ti) = 2.3 aB (top) and RMT(Ti) = 1.78 aB (bottom): Ti d (blue), dz2 (black), dx2y2 (gray), dxy (purple), dyz (yellow), dzx (green), and O p (red).

Close modal

The effect of the muffin-tin radius on the calculated charge transfer between Ti and O is studied by comparing the Bader analyses using the two different RMT(Ti) values (Table III). Using RMT(Ti) = 2.30 aB and PBE, the Bader analysis indicates that the charges on Ti and O are exactly the same as in the calculation with RMT(Ti) = 1.78 aB: +2.28e and −1.14e for Ti and O, respectively. After U is increased to 10 eV, the Bader charges on Ti and O are increased to +2.60e and −1.30e, respectively, indicating the transfer of 0.32e to both oxygen atoms. The amount is slightly larger than the charge transfer predicted using the small RMT(Ti) calculation. For anatase, the Bader charges on Ti and O are predicted to be slightly smaller than those of rutile using the PBE functional, regardless of the RMT(Ti) value. The smaller ionic charges are consistent with the shorter bond lengths observed in anatase. After U is increased to 10 eV using the small RMT(Ti), the Bader charges on Ti and O are increased to +2.56e and −1.28e, respectively, indicating the transfer of an additional 0.3e to both oxygen atoms. The calculation using the large RMT(Ti) similarly predicts a slightly larger charge transfer of 0.32e as in rutile.

The inclusion of the on-site Hubbard U interaction raises the total energies for rutile and anatase. However, it pushes up the energy of anatase by a greater amount than that of rutile for a given increment of U, resulting in a decrease in the stability gap (Fig. 7). Because both the large and small RMT(Ti) calculations predict the stability gap to reduce as U increases, we attribute the destabilization of anatase relative to rutile to the increased ionic bonding character in TiO2 induced by the Hubbard U. Using the Bader charges in Table III, we have calculated the Madelung potentials on the Ti and O sites (Table S2 in the supplementary material) as well as the total electrostatic lattice energy (Fig. 10).36 As the Hubbard U value increases, the contribution from the Madelung energy indeed increases for both rutile and anatase. In a more detailed examination, the difference in the lattice energy between rutile and anatase EMAEMR increases from 1.33 eV for U = 0 to 1.59 eV for U = 10 eV, as calculated using RMT(Ti) = 2.30 aB. The increased lattice energy difference points toward the reduced stability of anatase at higher U.

FIG. 10.

The Madelung lattice energy of rutile (black square) and anatase (blue circle) as a function of U, calculated using RMT(Ti) = 1.78 aB (open) and RMT(Ti) = 2.30 aB (filled).

FIG. 10.

The Madelung lattice energy of rutile (black square) and anatase (blue circle) as a function of U, calculated using RMT(Ti) = 1.78 aB (open) and RMT(Ti) = 2.30 aB (filled).

Close modal

On the other hand, the different rates of anatase destabilization remain unaccounted for in the present study, as both the large and small RMT(Ti) calculations yield practically the same result (Fig. 10). One possible explanation is that the Madelung calculation is the simplest one for evaluating the lattice electrostatic energy. It is based on the model of point charges or charges distributed with spherical symmetry, corresponding to the electric monopole term in the multipole expansion of electrostatic interaction. For a more rigorous assessment of the lattice energy, one can include higher-order terms allowed by the site symmetry, e.g., D2h for Ti (rutile) and C2v for O.53,69 It remains for a future study whether the inclusion of higher-order multipole terms in the lattice energy could explain more rapid destabilization of anatase relative to rutile calculated with larger RMT(Ti) and smaller RMT(O) (Fig. 7).

In order to alleviate the disconcerting sensitivity of the DFT+U calculations on the muffin-tin radius, one may consider a simple way to reduce the RMT dependence, for example, by dividing Ueff with the total occupancy nd. However, Nawa et al.70 pointed out that the on-site Coulomb interaction parameter Ueff itself can vary as much as 2–3 eV as the occupancy nd (or RMT) increases in their study of transition metal monoxides using linear response theory. In order to minimize the dependence of the electronic properties on RMT, they proposed to use suitably scaled values of Ueff for different RMT values. On the other hand, Wang and Jiang71 adopted the reformulated DFT+U approach and determined on-site Coulomb interaction parameters based on the Thomas–Fermi screening model, without U and J.72–74 They reported that the method produces the energy difference between the ferromagnetic and antiferromagnetic phases much less dependent on the size of RMT in a number of oxide systems. Further studies are clearly necessary to produce the DFT+U and like results independent of the size of the correlation space.

The dependence on RMT may be viewed alternatively with a related question: what should be the appropriate size of the correlation space in the DFT+U approach. The concept of the correlation space is not well-defined just like the atom in molecules or solids because both concepts have to deal with spatially extended charge densities. The comparison between the partial charge analysis inside the correlation sphere and the Bader charge analysis presented in this study suggests that one can apply a consistent treatment by choosing the two spaces as similar as possible. After all, there is no justification why one should choose to apply the Hubbard U correction in a space much smaller or larger than the atomic sphere for the rest of the local orbitals. In this perspective, the relative stability between rutile and anatase calculated using RMT(Ti) = 1.9 aB is included in Fig. 7. Our results show that the stability gap decreases more rapidly than that calculated with RMT(Ti) = 1.78 aB but does not get reversed as predicted by the calculation with RMT(Ti)= 2.3 aB.

We have performed a systematic investigation of TiO2 rutile and anatase phases using the PBE+U method. The facilitation of on-site Coulomb repulsion leads both rutile and anatase to be more ionic by transferring additional electrons from Ti to O. The enhanced ionic character with diminished covalency increases the Ti–O bond lengths. The accompanied expansion of lattice parameters and the reduction of the stability gap between rutile and anatase derive from the transition to higher ionicity in Ti and O. In addition, both rutile and anatase show substantial unevenness in occupancy among d sub-orbitals due to considerable covalency. Because the Coulomb repulsive potentials are dependent on the orbital occupation, the corresponding potential shifts are also unequal between the bonding eg and the nearly empty t2g orbitals.

By varying the spatial range of the local orbital projection via different muffin-tin radii of Ti, our study further elucidates qualitatively similar but quantitatively distinct effects of the RMT(Ti) values on the structural parameters, the electronic structures, and the relative stability. Since the occupancy of the local orbitals directly depends on the muffin-tin radius, the potential energy shift of the eg orbitals shows higher sensitivity than that of t2g orbitals to varying the muffin-tin radius. On the other hand, a larger muffin-tin radius produces a depleted region of charges in the t2g orbitals further into the interstitial region. Consequently, the calculations using different sizes for the correlation sphere can give rise to distinct maps of the charge distribution, which are accountable for the quantitative discrepancies in the structural parameters, the electronic structures, and the relative stability.

A proper local projection is indeed of critical importance in applying the Hubbard U interaction in order to obtain coherent and ascribable results in the DFT+U approach. As shown in this work, one possibility is to choose the range of local projection consistent with the electronic charge distribution for a given system, e.g., the Bader bond critical point. Another possibility is to utilize the MLWF, which is designed to capture the localization of the electrons for specific electronic characteristics that are sought after, such as the correlated subspace. Whether the correlated orbitals are treated using DFT+U methods or beyond, such as the combination of DFT and dynamics mean field theory, the results from this study should be valuable for the further development of first-principle modeling of strongly correlated systems.

See the supplementary material for additional data from rutile band structure calculations and k-point specific electron density plots, as well as the effect of U and RMT(Ti) on the density difference plots and the occupancies of the d orbitals for anatase.

K.P. is grateful for the helpful and stimulating discussions with P. Blaha and F. Tran. M.R. is thankful for the travel support from the Department of Physics and the School of Graduate Studies at Baylor University. The authors acknowledge C. Bell at the High Performance Computing Center of Baylor University for technical support.

The authors have no conflicts to disclose.

Kenneth Park: Conceptualization (lead); Data curation (lead); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Project administration (lead); Resources (equal); Supervision (lead); Validation (equal); Writing – original draft (lead); Writing – review & editing (lead). Manjula Raman: Conceptualization (lead); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (lead); Validation (equal); Visualization (lead); Writing – original draft (equal); Writing – review & editing (equal). Anjy-Joe Olatunbosun: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Software (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). Jared Pohlmann: Formal analysis (equal); Investigation (equal); Methodology (equal); Resources (equal); Software (equal); Writing – review & editing (equal).

The data that support the findings of this study are available within the article and its supplementary material.

1.
V. I.
Anisimov
,
F.
Aryasetiawan
, and
A. I.
Lichtenstein
, “
First-principles calculations of the electronic structure and spectra of strongly correlated systems: The LDA + U method
,”
J. Phys.: Condens. Matter
9
(
4
),
767
(
1997
).
2.
B.
Himmetoglu
,
A.
Floris
,
S.
de Gironcoli
, and
M.
Cococcioni
, “
Hubbard-corrected DFT energy functionals: The LDA+U description of correlated systems
,”
Int. J. Quantum Chem.
114
(
1
),
14
49
(
2014
).
3.
A. I.
Liechtenstein
,
V. I.
Anisimov
, and
J.
Zaanen
, “
Density-functional theory and strong interactions: Orbital ordering in mott-hubbard insulators
,”
Phys. Rev. B
52
(
8
),
R5467
(
1995
).
4.
S. L.
Dudarev
,
G. A.
Botton
,
S. Y.
Savrasov
,
C. J.
Humphreys
, and
A. P.
Sutton
, “
Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+U study
,”
Phys. Rev. B
57
(
3
),
1505
(
1998
).
5.
E.
Bousquet
and
N.
Spaldin
, “
J dependence in the LSDA+U treatment of noncollinear magnets
,”
Phys. Rev. B
82
(
22
),
220402
(
2010
).
6.
L.
de’Medici
, “
Hund’s coupling and its key role in tuning multiorbital correlations
,”
Phys. Rev. B
83
(
20
),
205112
(
2011
).
7.
Y.-C.
Wang
,
Z.-H.
Chen
, and
H.
Jiang
, “
The local projection in the density functional theory plus U approach: A critical assessment
,”
J. Chem. Phys.
144
(
14
),
144106
(
2016
).
8.
P.
Giannozzi
,
S.
Baroni
,
N.
Bonini
,
M.
Calandra
,
R.
Car
,
C.
Cavazzoni
,
D.
Ceresoli
,
G. L.
Chiarotti
,
M.
Cococcioni
,
I.
Dabo
et al, “
QUANTUM ESPRESSO: A modular and open-source software project for quantum simulations of materials
,”
J. Phys.: Condens. Matter
21
(
39
),
395502
(
2009
).
9.
G.
Kresse
and
J.
Furthmüller
, “
Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set
,”
Phys. Rev. B
54
(
16
),
11169
(
1996
).
10.
P.
Blaha
,
K.
Schwarz
,
G. K. H.
Madsen
,
D.
Kvasnicka
,
J.
Luitz
et al wien2k,
2001
.
11.
P.
Blaha
,
K.
Schwarz
,
F.
Tran
,
R.
Laskowski
,
G. K. H.
Madsen
, and
L. D.
Marks
, “
WIEN2k: An APW+lo program for calculating the properties of solids
,”
J. Chem. Phys.
152
(
7
),
074101
(
2020
).
12.
A.
Fujishima
and
K.
Honda
, “
Electrochemical photolysis of water at a semiconductor electrode
,”
Nature
238
(
5358
),
37
38
(
1972
).
13.
M.
Grätzel
, “
Dye-sensitized solar cells
,”
J. Photochem. Photobiol., C
4
(
2
),
145
153
(
2003
).
14.
U.
Diebold
, “
The surface science of titanium dioxide
,”
Surf. Sci. Rep.
48
,
53
(
2003
).
15.
X.
Chen
,
L.
Liu
,
P. Y.
Yu
, and
S. S.
Mao
, “
Increasing solar absorption for photocatalysis with black hydrogenated titanium dioxide nanocrystals
,”
Science
331
(
6018
),
746
750
(
2011
).
16.
J. C.
Matsubu
,
S.
Zhang
,
L.
DeRita
,
N. S.
Marinkovic
,
J. G.
Chen
,
G. W.
Graham
,
X.
Pan
, and
P.
Christopher
, “
Adsorbate-mediated strong metal–support interactions in oxide-supported Rh catalysts
,”
Nat. Chem.
9
(
2
),
120
127
(
2017
).
17.
M. N.
Chong
,
B.
Jin
,
C. W. K.
Chow
, and
C.
Saint
, “
Recent developments in photocatalytic water treatment technology: A review
,”
Water Res.
44
(
10
),
2997
3027
(
2010
).
18.
Q.
Guo
,
C.
Zhou
,
Z.
Ma
, and
X.
Yang
, “
Fundamentals of TiO2 photocatalysis: Concepts, mechanisms, and challenges
,”
Adv. Mater.
31
(
50
),
1901997
(
2019
).
19.
R.
Lang
,
X.
Du
,
Y.
Huang
,
X.
Jiang
,
Q.
Zhang
,
Y.
Guo
,
K.
Liu
,
B.
Qiao
,
A.
Wang
, and
T.
Zhang
, “
Single-atom catalysts based on the metal–oxide interaction
,”
Chem. Rev.
120
(
21
),
11986
12043
(
2020
).
20.
B.
O’regan
and
M.
Grätzel
, “
A low-cost, high-efficiency solar cell based on dye-sensitized colloidal TiO2 films
,”
Nature
353
(
6346
),
737
740
(
1991
).
21.
K.
Lee
,
A.
Mazare
, and
P.
Schmuki
, “
One-dimensional titanium dioxide nanomaterials: Nanotubes
,”
Chem. Rev.
114
(
19
),
9385
9454
(
2014
).
22.
J.
Wen
,
X.
Li
,
W.
Liu
,
Y.
Fang
,
J.
Xie
, and
Y.
Xu
, “
Photocatalysis fundamentals and surface modification of TiO2 nanomaterials
,”
Chin. J. Catal.
36
(
12
),
2049
2070
(
2015
).
23.
J. C.
Matsubu
,
V. N.
Yang
, and
P.
Christopher
, “
Isolated metal active site concentration and stability control catalytic CO2 reduction selectivity
,”
J. Am. Chem. Soc.
137
(
8
),
3076
3084
(
2015
).
24.
J.
Xu
,
L.
Sun
,
X.
Qi
,
Z.
Wang
,
Q.
Fu
, and
C.
Pan
, “
A novel strategy to enhance the multiple interface effect using amorphous carbon packaged hydrogenated TiO2 for stable and effective microwave absorption
,”
J. Mater. Chem. C
7
(
20
),
6152
6160
(
2019
).
25.
K. T.
Park
,
M. H.
Pan
,
V.
Meunier
, and
E. W.
Plummer
, “
Surface reconstructions of TiO2(110) driven by suboxides
,”
Phys. Rev. Lett.
96
(
22
),
226105
(
2006
).
26.
K. T.
Park
,
M.
Pan
,
V.
Meunier
, and
E. W.
Plummer
, “
Reoxidation of TiO2(110) via Ti interstitials and line defects
,”
Phys. Rev. B
75
(
24
),
245415
(
2007
).
27.
K.
Zhu
,
Y.
Xia
,
Z.
Zhang
, and
K.
Park
, “
Probing structure of cross-linked (1 × 2) rutile TiO2(110): Adsorption of trimethyl acetic acid
,”
J. Phys. Chem. C
120
(
28
),
15257
15264
(
2016
).
28.
Y.
Xia
,
K.
Zhu
,
Z.
Zhang
, and
K.
Park
, “
Photo-stimulated desorption of trimethyl acetic acid on cross-linked (1 × 2) TiO2(1 1 0) probed by scanning tunneling microscopy
,”
Appl. Surf. Sci.
511
,
145553
(
2020
).
29.
A. A.
Levchenko
,
G.
Li
,
J.
Boerio-Goates
,
B. F.
Woodfield
, and
A.
Navrotsky
, “
TiO2 stability landscape: Polymorphism, surface energy, and bound water energetics
,”
Chem. Mater.
18
(
26
),
6324
6332
(
2006
).
30.
M. R.
Ranade
,
A.
Navrotsky
,
H. Z.
Zhang
,
J. F.
Banfield
,
S. H.
Elder
,
A.
Zaban
,
P. H.
Borse
,
S. K.
Kulkarni
,
G. S.
Doran
, and
H. J.
Whitfield
, “
Energetics of nanocrystalline TiO2
,”
Proc. Natl. Acad. Sci. U. S. A.
99
(
suppl_2
),
6476
6481
(
2002
).
31.
Y.
Zhang
,
J. W.
Furness
,
B.
Xiao
, and
J.
Sun
, “
Subtlety of TiO2 phase stability: Reliability of the density functional theory predictions and persistence of the self-interaction error
,”
J. Chem. Phys.
150
(
1
),
014105
(
2019
).
32.
M. T.
Curnan
and
J. R.
Kitchin
, “
Investigating the energetic ordering of stable and metastable TiO2 polymorphs using DFT+U and hybrid functionals
,”
J. Phys. Chem. C
119
(
36
),
21060
21071
(
2015
).
33.
N. H.
Vu
,
H. V.
Le
,
T. M.
Cao
,
V. V.
Pham
,
H. M.
Le
, and
D.
Nguyen-Manh
, “
Anatase–rutile phase transformation of titanium dioxide bulk material: A DFT+U approach
,”
J. Phys.: Condens. Matter
24
(
40
),
405501
(
2012
).
34.
M. E.
Arroyo-de Dompablo
,
A.
Morales-García
, and
M.
Taravillo
, “
DFT+U calculations of crystal lattice, electronic structure, and phase stability under pressure of TiO2 polymorphs
,”
J. Chem. Phys.
135
(
5
),
054503
(
2011
).
35.
J. K.
Burdett
,
T.
Hughbanks
,
G. J.
Miller
,
J. W.
Richardson
, Jr.
, and
J. V.
Smith
, “
Structural-electronic relationships in inorganic solids: Powder neutron diffraction studies of the rutile and anatase polymorphs of titanium dioxide at 15 and 295 K
,”
J. Am. Chem. Soc.
109
(
12
),
3639
3646
(
1987
).
36.
K.
Momma
and
F.
Izumi
, “
VESTA 3 for three-dimensional visualization of crystal, volumetric and morphology data
,”
J. Appl. Crystallogr.
44
(
6
),
1272
1276
(
2011
).
37.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
, “
Generalized gradient approximation made simple
,”
Phys. Rev. Lett.
77
(
18
),
3865
(
1996
).
38.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
, “
Perdew, Burke, and Ernzerh of Reply:
,”
Phys. Rev. Lett.
80
(
4
),
891
(
1998
).
39.
V. I.
Anisimov
,
I. V.
Solovyev
,
M. A.
Korotin
,
M. T.
Czyżyk
, and
G. A.
Sawatzky
, “
Density-functional theory and NiO photoemission spectra
,”
Phys. Rev. B
48
(
23
),
16929
(
1993
).
40.
Q.
Zhao
and
H. J.
Kulik
, “
Where does the density localize in the solid state? Divergent behavior for hybrids and DFT+U
,”
J. Chem. Theory Comput.
14
(
2
),
670
683
(
2018
).
41.
E.
German
,
R.
Faccio
, and
A. W.
Mombrú
, “
A DFT+U study on structural, electronic, vibrational and thermodynamic properties of TiO2 polymorphs and hydrogen titanate: Tuning the Hubbard ‘U-term’
,”
J. Phys. Commun.
1
(
5
),
055006
(
2017
).
42.
F.
Labat
,
P.
Baranek
,
C.
Domain
,
C.
Minot
, and
C.
Adamo
, “
Density functional theory analysis of the structural and electronic properties of TiO2 rutile and anatase polytypes: Performances of different exchange-correlation functionals
,”
J. Chem. Phys.
126
(
15
),
154703
(
2007
).
43.
Y.-f.
Zhang
,
W.
Lin
,
Y.
Li
,
K.-n.
Ding
, and
J.-q.
Li
, “
A theoretical study on the electronic structures of TiO2: Effect of Hartree–Fock exchange
,”
J. Phys. Chem. B
109
(
41
),
19270
19277
(
2005
).
44.
Z.
Hu
and
H.
Metiu
, “
Choice of U for DFT+U calculations for titanium oxides
,”
J. Phys. Chem. C
115
(
13
),
5841
5845
(
2011
).
45.
D. C.
Cronemeyer
, “
Electrical and optical properties of rutile single crystals
,”
Phys. Rev.
87
(
5
),
876
(
1952
).
46.
H.
Tang
,
H.
Berger
,
P. E.
Schmid
,
F.
Levy
, and
G.
Burri
, “
Photoluminescence in TiO2 anatase single crystals
,”
Solid State Commun.
87
(
9
),
847
850
(
1993
).
47.
J. P.
Perdew
and
M.
Levy
, “
Physical content of the exact Kohn-Sham orbital energies: Band gaps and derivative discontinuities
,”
Phys. Rev. Lett.
51
(
20
),
1884
(
1983
).
48.
J. B.
Goodenough
, “
Metallic oxides
,”
Prog. Solid State Chem.
5
,
145
399
(
1971
).
49.
P.
Cox
,
Transition Metal Oxides
(
Oxford Science Publications
,
1992
).
50.
P. I.
Sorantin
and
K.
Schwarz
, “
Chemical bonding in rutile-type compounds
,”
Inorg. Chem.
31
(
4
),
567
576
(
1992
).
51.
F. M. F.
De Groot
,
M.
Grioni
,
J. C.
Fuggle
,
J.
Ghijsen
,
G. A.
Sawatzky
, and
H.
Petersen
, “
Oxygen 1s X-ray-absorption edges of transition-metal oxides
,”
Phys. Rev. B
40
(
8
),
5715
(
1989
).
52.
E.
Stoyanov
,
F.
Langenhorst
, and
G.
Steinle-Neumann
, “
The effect of valence state and site geometry on Ti L3,2 and O K electron energy-loss spectra of TixOy phases
,”
Am. Mineral.
92
(
4
),
577
586
(
2007
).
53.
B.
Jiang
,
J. M.
Zuo
,
N.
Jiang
,
M.
O’keeffe
, and
J. C. H.
Spence
, “
Charge density and chemical bonding in rutile, TiO2
,”
Acta Crystallogr., Sect. A: Found. Crystallogr.
59
(
4
),
341
350
(
2003
).
54.
R.
Asahi
,
Y.
Taga
,
W.
Mannstadt
, and
A. J.
Freeman
, “
Electronic and optical properties of anatase TiO2
,”
Phys. Rev. B
61
(
11
),
7459
(
2000
).
55.
R. F. W.
Bader
, “
Atoms in molecules
,”
Acc. Chem. Res.
18
(
1
),
9
15
(
1985
).
56.
J.
Sofo
and
G.
Garcia
, “
Implementation of Bader theory in Wien package
,” http://www.wien2k.at/lapw/; accessed 15 12 2023.
57.
M.
Yoshiya
,
I.
Tanaka
,
K.
Kaneko
, and
H.
Adachi
, “
First principles calculation of chemical shifts in ELNES/NEXAFS of titanium oxides
,”
J. Phys.: Condens. Matter
11
(
16
),
3217
(
1999
).
58.
C. E.
Rice
and
W. R.
Robinson
, “
High-temperature crystal chemistry of Ti2O3: Structural changes accompanying the semiconductor–metal transition
,”
Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem.
33
(
5
),
1342
1348
(
1977
).
59.
The Materials Project, Materials data on Ti2O3 by materials project.
60.
G. V.
Samsonov
,
The Oxide Handbook
, 2nd ed.
(IFI/Plenum, 1982), translated from Russian by
R. K.
Johnston
.
61.
Z.
Wang
,
C.
Brock
,
A.
Matt
, and
K. H.
Bevan
, “
Implications of the DFT+U method on polaron properties in energy materials
,”
Phys. Rev. B
96
(
12
),
125150
(
2017
).
62.
W. E.
Pickett
,
S. C.
Erwin
, and
E. C.
Ethridge
, “
Reformulation of the LDA+U method for a local-orbital basis
,”
Phys. Rev. B
58
(
3
),
1201
(
1998
).
63.
S.
Fabris
,
G.
Vicario
,
G.
Balducci
,
S.
de Gironcoli
, and
S.
Baroni
, “
Electronic and atomistic structures of clean and reduced ceria surfaces
,”
J. Phys. Chem. B
109
(
48
),
22860
22867
(
2005
).
64.
C.
Tablero
, “
Representations of the occupation number matrix on the LDA/GGA+U method
,”
J. Phys.: Condens. Matter
20
(
32
),
325205
(
2008
).
65.
G.
Geneste
,
B.
Amadon
,
M.
Torrent
, and
G.
Dezanneau
, “
DFT+U study of self-trapping, trapping, and mobility of oxygen-type hole polarons in barium stannate
,”
Phys. Rev. B
96
(
13
),
134123
(
2017
).
66.
R.
Outerovitch
and
B.
Amadon
, “
Electronic interaction Upp on oxygen p orbitals in oxides: Role of correlated orbitals on the example of UO2 and TiO2
,”
Phys. Rev. B
107
(
23
),
235126
(
2023
).
67.
S.-G.
Park
,
B.
Magyari-Köpe
, and
Y.
Nishi
, “
Electronic correlation effects in reduced rutile TiO2 within the LDA+U method
,”
Phys. Rev. B
82
(
11
),
115109
(
2010
).
68.
H.
Park
,
A. J.
Millis
, and
C. A.
Marianetti
, “
Density functional versus spin-density functional and the choice of correlated subspace in multivariable effective action theories of electronic structure
,”
Phys. Rev. B
92
(
3
),
035146
(
2015
).
69.
M.
Kara
and
K.
Kurki-Suonio
, “
Symmetrized multipole analysis of orientational distributions
,”
Acta Crystallogr., Sect. A: Found. Adv.
37
(
2
),
201
210
(
1981
).
70.
K.
Nawa
,
T.
Akiyama
,
T.
Ito
,
K.
Nakamura
,
T.
Oguchi
, and
M.
Weinert
, “
Scaled effective on-site Coulomb interaction in the DFT+U method for correlated materials
,”
Phys. Rev. B
97
(
3
),
035117
(
2018
).
71.
Y.-C.
Wang
and
H.
Jiang
, “
Local screened Coulomb correction approach to strongly correlated d-electron systems
,”
J. Chem. Phys.
150
(
15
),
154116
(
2019
).
72.
F.
Bultmark
,
F.
Cricchio
,
O.
Grånäs
, and
L.
Nordström
, “
Multipole decomposition of LDA+U energy and its application to actinide compounds
,”
Phys. Rev. B
80
(
3
),
035121
(
2009
).
73.
M. S. S.
Brooks
, “
Thomas-Fermi screening of exchange interactions
,”
J. Phys.: Condens. Matter
13
(
22
),
L469
(
2001
).
74.
M. R.
Norman
, “
Calculation of effective Coulomb interaction for Pr3+, U4+, and UPt3
,”
Phys. Rev. B
52
(
3
),
1421
(
1995
).