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 TiO_{2} 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.

## I. INTRODUCTION

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 *F*^{0} in the calculation of electron–electron interactions.^{4} In this simplification, only *U*_{eff} ≡ *U* − *J* 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 MnO_{2} 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 TiO_{2} rutile and anatase phases using all-electron WIEN2k. TiO_{2} 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 TiO_{2} 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 TiO_{2} 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, TiO_{2} 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 (*P*4_{2}/*mnm*) whereas anatase to the body-centered tetragonal system (*I*4_{1}/*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.

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 *R*_{MT}(Ti) = 1.78 *a*_{B} (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 *R*_{MT}(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.

## II. METHOD AND COMPUTATIONAL DETAILS

^{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,

*u*

_{l}(

*r*,

*E*

_{l}) is a solution of the radial Schrödinger equation for energy

*E*

_{l}, and $u\u0307l(r,El)$ is the energy derivative of

*u*

_{l}at the energy

*E*

_{l}. 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,

*lo*, the standard APW basis is used with the usual fixed energy

*E*

_{l}along with a new local orbital (

*lo*) for variational flexibility,

*A*

_{lm}and

*B*

_{lm}are now independent of

**k**

_{n}but determined by the requirement that $\varphi 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 *a*_{B}, respectively, to represent a reference system with a small correlation space. The accuracy of the calculation is controlled by the dimensionless cutoff parameter *R*_{MT}*K*_{max}, which determines the size of the basis through the product of the smallest *R*_{MT} value and the largest *K* vector in the plane wave expansion. The value of *R*_{MT}*K*_{max} 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 TiO_{2} 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 *R*_{MT}(Ti) values were systematically increased from 1.78 to 2.30 *a*_{B}. Because the L/APW method requires non-overlapping spheres for atomic regions, *R*_{MT} values for oxygen were reduced accordingly from 1.61 to 1.30 *a*_{B}. For larger *R*_{MT} (Ti) values, the pronounced dissimilarity in muffin-tin radii between Ti and O could be problematic, leading to very different effective *R*_{MT}*K*_{max} values. To avoid the linearization of the basis, *R*_{MT}*K*_{max} was thus reduced to 7.5, and a high energy local orbital was included for the largest *R*_{MT}(*Ti*) = 2.3 *a*_{B}. 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 *R*_{MT}-independence of the GGA calculations.

*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}

^{,}

*l*and spin index

*σ*, and

*V*

_{ee}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

*E*

_{dc}. The most widely used method of eliminating double counting for strongly correlated insulating systems was applied

^{39}

*n*=

*n*

^{↑}+

*n*

^{↓}are the occupation numbers of localized orbitals.

*F*

^{0}[4,2] in evaluating the interaction matrix elements ⟨

*m*,

*m*″|

*V*

_{ee}|

*m*′,

*m*‴⟩. Under this isotropic scheme of the DFT+U, the Hubbard correction with the double counting term included is written as

*U*

_{eff}regarded as an effective parameter, including the exchange

*J*such that

*U*

_{eff}=

*U*−

*J*. Consequently, by differentiating with respect to each orbital occupancy, one obtains the additional orbital-dependent potential energy

## III. RESULTS AND DISCUSSION

### A. The effect of Hubbard 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.

. | . | . | Lattice constants . | . | |
---|---|---|---|---|---|

Phase . | Method . | U (eV) . | a (Å) . | c (Å) . | c/a . |

Rutile | GGA | 0 | 4.647 | 2.966 | 0.638 |

GGA+U | 3 | 4.650 (+0.06) | 2.985 (+0.64) | 0.642 | |

GGA+U | 5 | 4.649 (+0.04) | 2.995 (+0.98) | 0.644 | |

GGA+U | 8 | 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 | 0 | 3.801 | 9.705 | 2.553 |

GGA+U | 3 | 3.819 (+0.47) | 9.693 (−0.12) | 2.538 | |

GGA+U | 5 | 3.828 (+0.71) | 9.695 (−0.10) | 2.532 | |

GGA+U | 8 | 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 . | . | |
---|---|---|---|---|---|

Phase . | Method . | U (eV) . | a (Å) . | c (Å) . | c/a . |

Rutile | GGA | 0 | 4.647 | 2.966 | 0.638 |

GGA+U | 3 | 4.650 (+0.06) | 2.985 (+0.64) | 0.642 | |

GGA+U | 5 | 4.649 (+0.04) | 2.995 (+0.98) | 0.644 | |

GGA+U | 8 | 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 | 0 | 3.801 | 9.705 | 2.553 |

GGA+U | 3 | 3.819 (+0.47) | 9.693 (−0.12) | 2.538 | |

GGA+U | 5 | 3.828 (+0.71) | 9.695 (−0.10) | 2.532 | |

GGA+U | 8 | 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.

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.

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 studies^{34,42–44} but well short of the experimental values of 3.03 eV^{45} 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}

The valence band is mainly made of O 2*p* states. The empty conduction band consists of predominantly Ti 3*d* 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 *t*_{2g} band and the upper one as the Ti *e*_{g} band in approximating the site symmetry of Ti as the octahedral symmetry *O*_{h}.^{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 *d*_{xy}, *d*_{yz}, and *d*_{zx} characters. Likewise, the upper conduction band is predominantly composed of $dz2$ and $dx2\u2212y2$ characters, indeed justifying the assignment of the *O*_{h} symmetry. The calculated crystal field splitting (CFS) of 2.3 eV between the *t*_{2g} and *e*_{g} bands is in fair agreement with the experimental value of 2.6 eV.^{51}

It is notable that the Ti 3*d* states and the O 2*p* states hybridize substantially in the middle and lower valence band regions.^{50,54} The overlap between the Ti *e*_{g} states and O 2*p* states around −4 eV mainly represents the *σ* bonding in the molecular orbital (MO) scheme. In the higher energy region of the valence band, Ti *d*_{yz} and *d*_{zx} form weak *π*-bonding with O 2*p* states, and at the top of the valence band, non-bonding Ti *d*_{xy} 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 *e*_{g} states and O 2*p* 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 *t*_{2g} and *e*_{g} 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 *t*_{2g} and *e*_{g} sub-bands are affected disproportionately by the Hubbard U interaction term. For example, the anti-bonding *π** *t*_{2g} band is pushed up to a greater extent relative to the *σ** *e*_{g} 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 method^{34} 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 Kulik^{40} reported that the CFS between the two *d* sub-bands decreases even more rapidly as the U value increases. The *t*_{2g} and *e*_{g} 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 *t*_{2g}-type orbital symmetry.^{50} For oxygen, three distinct regions of electron surplus along the three nearby Ti atoms represent the *sp*^{2}-type bonding for which oxygen forms a planar O–Ti_{3} cluster. One also observes excess electronic charges at Ti directed along the octahedrally coordinated O atoms, corresponding to the *e*_{g} orbitals observed in PDOS. The significant charge accumulation between Ti and O along the bond directions clearly supports substantial covalency through *e*_{g} 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}

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 (*e*_{g} 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 *a*_{B}. The occupancy in the *e*_{g} orbitals is much higher than that in the *t*_{2g} 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.12*e* per spin or 0.24*e* 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).

U(R_{MT}(Ti) = 1.78 a_{B})
. | Total d
. | $dz2$ . | $dx2\u2212y2$ . | d_{xy}
. | d_{yz}
. | d_{zx}
. |
---|---|---|---|---|---|---|

0 | 0.516 | 0.139 | 0.148 | 0.070 | 0.078 | 0.079 |

3 | 0.481 | 0.134 | 0.145 | 0.062 | 0.070 | 0.071 |

5 | 0.459 | 0.130 | 0.140 | 0.057 | 0.065 | 0.066 |

8 | 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(R_{MT}(Ti) = 2.3 a_{B}) | Total d | $dz2$ | $dx2\u2212y2$ | d_{xy} | d_{yz} | d_{zx} |

0 | 0.759 | 0.207 | 0.223 | 0.104 | 0.112 | 0.112 |

3 | 0.709 | 0.205 | 0.217 | 0.089 | 0.098 | 0.099 |

5 | 0.680 | 0.202 | 0.219 | 0.080 | 0.089 | 0.089 |

8 | 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(R_{MT}(Ti) = 1.78 a_{B})
. | Total d
. | $dz2$ . | $dx2\u2212y2$ . | d_{xy}
. | d_{yz}
. | d_{zx}
. |
---|---|---|---|---|---|---|

0 | 0.516 | 0.139 | 0.148 | 0.070 | 0.078 | 0.079 |

3 | 0.481 | 0.134 | 0.145 | 0.062 | 0.070 | 0.071 |

5 | 0.459 | 0.130 | 0.140 | 0.057 | 0.065 | 0.066 |

8 | 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(R_{MT}(Ti) = 2.3 a_{B}) | Total d | $dz2$ | $dx2\u2212y2$ | d_{xy} | d_{yz} | d_{zx} |

0 | 0.759 | 0.207 | 0.223 | 0.104 | 0.112 | 0.112 |

3 | 0.709 | 0.205 | 0.217 | 0.089 | 0.098 | 0.099 |

5 | 0.680 | 0.202 | 0.219 | 0.080 | 0.089 | 0.089 |

8 | 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.28*e* and −1.14*e*, respectively, using the PBE functional. These values are in good agreement with previously calculated values, for instance, +2.24*e* for Ti.^{40} As U increases to 10 eV, the Bader charges increase to +2.58*e* and −1.29*e*, respectively, indicating that Ti loses an additional 0.3*e* to two O atoms. This value is also comparable to about 0.25*e* calculated using the pseudopotential method.^{40}

Rutile . | . | Bader charge (e)
. | Critical points (3, −1) from (in Å) . | ||
---|---|---|---|---|---|

R_{MT}(Ti) (a_{B})
. | U (eV) . | Ti . | O . | Ti; O(1) . | Ti; O(2) . |

1.78 | 0 | +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 | 0 | +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 (e) | Critical points (3, −1) from (in Å) | |||

1.78 | 0 | +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 | 0 | +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 |

Rutile . | . | Bader charge (e)
. | Critical points (3, −1) from (in Å) . | ||
---|---|---|---|---|---|

R_{MT}(Ti) (a_{B})
. | U (eV) . | Ti . | O . | Ti; O(1) . | Ti; O(2) . |

1.78 | 0 | +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 | 0 | +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 (e) | Critical points (3, −1) from (in Å) | |||

1.78 | 0 | +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 | 0 | +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.24*e* 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 *R*_{MT}(Ti) = 1.78 *a*_{B} (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 TiO_{2} 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 $dx2\u2212y2$ 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 *d*_{xy}, *d*_{yz}, and *d*_{zx} 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 *e*_{g} and *t*_{2g} orbitals is reflected in the density difference map (Fig. 5), which shows that the depleted region of *t*_{2g} symmetry (blue) has enlarged to a greater extent than the region of *e*_{g} 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.

. | . | Rutile . | Anatase . | ||
---|---|---|---|---|---|

Method . | U (eV) . | Ti–O(1) (Å) . | Ti–O(2) (Å) . | Ti–O(1) (Å) . | Ti–O(2) (Å) . |

GGA | 0 | 2.006 | 1.958 | 1.995 | 1.948 |

GGA+U | 3 | 2.007 | 1.966 | 2.002 | 1.955 |

GGA+U | 5 | 2.006 | 1.970 | 2.003 | 1.959 |

GGA+U | 8 | 2.007 | 1.977 | 2.006 | 1.966 |

GGA+U | 10 | 2.008 | 1.981 | 2.008 | 1.970 |

. | . | Rutile . | Anatase . | ||
---|---|---|---|---|---|

Method . | U (eV) . | Ti–O(1) (Å) . | Ti–O(2) (Å) . | Ti–O(1) (Å) . | Ti–O(2) (Å) . |

GGA | 0 | 2.006 | 1.958 | 1.995 | 1.948 |

GGA+U | 3 | 2.007 | 1.966 | 2.002 | 1.955 |

GGA+U | 5 | 2.006 | 1.970 | 2.003 | 1.959 |

GGA+U | 8 | 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 TiO_{2}, Ti_{2}O_{3}, 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 Ti_{2}O_{3}, 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 TiO_{2}, 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 *e*_{g} as well as non-bonding *t*_{2g} 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 $dx2\u2212y2$ predicted by GGA (Table II).

### B. The effect of the local-projection size

The comparison of the structural parameters and the electronic structures of TiO_{2} 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 *a*_{B}: *R*_{MT}(Ti) = 1.78, *R*_{MT}(O) = 1.61; *R*_{MT}(Ti) = 1.91, *R*_{MT}(O) = 1.73; *R*_{MT}(Ti) = 2.1, *R*_{MT}(O) = 1.5; *R*_{MT}(Ti) = 2.3, *R*_{MT}(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 *R*_{MT}(Ti) = 1.78 (Fig. 2). It increases from 4.655 Å with *R*_{MT}(Ti) = 1.78 *a*_{B} to 4.696 Å with *R*_{MT}(Ti) = 2.3 *a*_{B}. The lattice constant *c* also shows its increasing value with a larger muffin-tin radius of Ti, from 3.021 Å with *R*_{MT}(Ti) = 1.78 *a*_{B} to 3.077 Å with *R*_{MT}(Ti) = 2.3 *a*_{B} [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 *R*_{MT}(Ti) = 2.3 *a*_{B} 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 *R*_{MT}(Ti) = 1.78 *a*_{B} to *R*_{MT}(Ti) = 2.3 *a*_{B} [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}

The lattice parameters obtained using the largest muffin-tin radius of 2.3 *a*_{B} 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 *R*_{PAW} = 2.5 *a*_{B}.^{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 *R*_{MT}(Ti) = 2.5 *a*_{B}. However, a bigger Ti muffin-tin sphere would necessitate further reduction of *R*_{MT}(O) in the APW scheme, risking approximate linear dependencies of the basis due to very different effective *R*_{MT}*K*_{max}. We believe that *R*_{MT}(Ti) = 2.3 *a*_{B} is about the maximum size for trustworthy results. Nevertheless, the unmistakable trend of the lattice constants with an increasing *R*_{MT}(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 *R*_{MT}(Ti) = 1.78 *a*_{B} and *R*_{MT}(Ti) = 2.3 *a*_{B} (Fig. 7). With the PBE functional and the small *R*_{MT} 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: *E*_{rutile} < *E*_{anatase} < *E*_{columbite}. 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: *E*_{rutile} < *E*_{anatase} < *E*_{brookite} < *E*_{columbite}.^{32}

Using *R*_{MT}(Ti) = 2.3 *a*_{B} (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 *R*_{MT} 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-workers^{61} 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 TiO_{2}, the polaronic gap state energy, the formation energy, and the polaronic hopping barrier all showed a strong dependence, as *R*_{PAW} for Ti varied from 1.9 to 2.8 *a*_{B}. 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 Amadon^{66} compared the truncated/renormalized atomic orbitals with the projected localized Wannier functions as the local orbitals in the study of UO_{2} and TiO_{2}. Using the U values obtained from the constrained random phase approximation (cRPA), they reported that the equilibrium volume can either expand (as previously reported^{67}) 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 *R*_{MT}(Ti) on our PBE+U results, the difference electron density map calculated with U = 8 eV using *R*_{MT}(Ti) = 2.30 *a*_{B} is stacked on top of that using *R*_{MT}(Ti) = 1.78 *a*_{B} in Fig. 8. The result using the large *R*_{MT}(Ti) value exhibits a similar pattern of overall charge transfer from Ti to O as that using the small *R*_{MT}(Ti). However it is clearly noticeable that the calculation with large *R*_{MT}(Ti) predicts significantly more charges remaining near the Ti atom in the *e*_{g} 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.

The calculation using *R*_{MT}(Ti) = 2.30 *a*_{B} shows that the occupancy in all *d* orbitals decreases with increasing *U*, as observed with *R*_{MT}(Ti) = 1.78 *a*_{B}. However, a detailed examination reveals that the occupancy in the *σ*-bonding *e*_{g} orbitals ($dz2$ and $dx2\u2212y2$) 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 *R*_{MT} value. On the other hand, the occupancy in the *π*-type *t*_{2g} orbitals (*d*_{xy}, *d*_{yz}, *d*_{zx}) decreases to 58%–61% of the initial values, more than 64%–67% of the initial values using the small *R*_{MT} = 1.78 *a*_{B}. Therefore, the calculation with the large *R*_{MT}(Ti) produces a greater variation in occupancy between the *e*_{g} and *t*_{2g} orbitals upon the U term. The large variation is also seen in the charge density difference map, where the electron depletion in *t*_{2g} orbitals reaches further into the interstitial regions between two oxygens in the calculation with the large *R*_{MT}(Ti) (Fig. 8).

The larger disparity in occupancy between the *e*_{g} and *t*_{2g} orbitals can also explain a bigger reduction in the CFS predicted by the calculation using the large *R*_{MT}(Ti) value. It is evident from Eq. (8) that the larger the occupancies in *e*_{g} orbitals, the smaller the Hubbard potential shifts $Veg$ are. On the contrary, the potential shifts $Vt2g$ of the *t*_{2g} 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 *R*_{MT}(Ti) = 2.3 *a*_{B}, compared to the calculation using *R*_{MT}(Ti) = 1.78 *a*_{B}. 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 *t*_{2g} and the O 2*p* orbitals. Therefore, the bandgap is less sensitive to the size of the correlation sphere for Ti than the crystal field splitting.

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 *R*_{MT}(Ti) values (Table III). Using *R*_{MT}(Ti) = 2.30 *a*_{B} and PBE, the Bader analysis indicates that the charges on Ti and O are exactly the same as in the calculation with *R*_{MT}(Ti) = 1.78 *a*_{B}: +2.28*e* and −1.14*e* for Ti and O, respectively. After U is increased to 10 eV, the Bader charges on Ti and O are increased to +2.60*e* and −1.30*e*, respectively, indicating the transfer of 0.32*e* to both oxygen atoms. The amount is slightly larger than the charge transfer predicted using the small *R*_{MT}(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 *R*_{MT}(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 *R*_{MT}(Ti), the Bader charges on Ti and O are increased to +2.56*e* and −1.28*e*, respectively, indicating the transfer of an additional 0.3*e* to both oxygen atoms. The calculation using the large *R*_{MT}(Ti) similarly predicts a slightly larger charge transfer of 0.32*e* 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 *R*_{MT}(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 TiO_{2} 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 $EMA\u2212EMR$ increases from 1.33 eV for U = 0 to 1.59 eV for U = 10 eV, as calculated using *R*_{MT}(Ti) = 2.30 *a*_{B}. The increased lattice energy difference points toward the reduced stability of anatase at higher U.

On the other hand, the different rates of anatase destabilization remain unaccounted for in the present study, as both the large and small *R*_{MT}(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., *D*_{2h} for Ti (rutile) and *C*_{2v} 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 *R*_{MT}(Ti) and smaller *R*_{MT}(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 *R*_{MT} dependence, for example, by dividing *U*_{eff} with the total occupancy *n*_{d}. However, Nawa *et al.*^{70} pointed out that the on-site Coulomb interaction parameter *U*_{eff} itself can vary as much as 2–3 eV as the occupancy *n*_{d} (or *R*_{MT}) increases in their study of transition metal monoxides using linear response theory. In order to minimize the dependence of the electronic properties on *R*_{MT}, they proposed to use suitably scaled values of *U*_{eff} for different *R*_{MT} values. On the other hand, Wang and Jiang^{71} 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 *R*_{MT} 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 *R*_{MT} 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 *R*_{MT}(Ti) = 1.9 *a*_{B} is included in Fig. 7. Our results show that the stability gap decreases more rapidly than that calculated with *R*_{MT}(Ti) = 1.78 *a*_{B} but does not get reversed as predicted by the calculation with *R*_{MT}(Ti)= 2.3 *a*_{B}.

## IV. CONCLUSIONS

We have performed a systematic investigation of TiO_{2} 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 *e*_{g} and the nearly empty *t*_{2g} 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 *R*_{MT}(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 *e*_{g} orbitals shows higher sensitivity than that of *t*_{2g} orbitals to varying the muffin-tin radius. On the other hand, a larger muffin-tin radius produces a depleted region of charges in the *t*_{2g} 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.

## SUPPLEMENTARY MATERIAL

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.

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

## REFERENCES

**LDA**+

*U*method

*J*dependence in the LSDA+

*U*treatment of noncollinear magnets

*U*approach: A critical assessment

*ab initio*total-energy calculations using a plane-wave basis set

_{2}photocatalysis: Concepts, mechanisms, and challenges

_{2}films

_{2}nanomaterials

_{2}reduction selectivity

_{2}for stable and effective microwave absorption

_{2}(110) driven by suboxides

_{2}(110) via Ti interstitials and line defects

_{2}(110): Adsorption of trimethyl acetic acid

_{2}(1 1 0) probed by scanning tunneling microscopy

_{2}stability landscape: Polymorphism, surface energy, and bound water energetics

_{2}

_{2}phase stability: Reliability of the density functional theory predictions and persistence of the self-interaction error

_{2}polymorphs using DFT+

*U*and hybrid functionals

*U*approach

*U*calculations of crystal lattice, electronic structure, and phase stability under pressure of TiO

_{2}polymorphs

*VESTA*3 for three-dimensional visualization of crystal, volumetric and morphology data

*U*study on structural, electronic, vibrational and thermodynamic properties of TiO

_{2}polymorphs and hydrogen titanate: Tuning the Hubbard ‘

*U*-term’

_{2}rutile and anatase polytypes: Performances of different exchange-correlation functionals

_{2}: Effect of Hartree–Fock exchange

*U*for DFT+

*U*calculations for titanium oxides

_{2}anatase single crystals

*s*X-ray-absorption edges of transition-metal oxides

_{3,2}and O K electron energy-loss spectra of Ti

_{x}O

_{y}phases

_{2}

_{2}

_{2}O

_{3}: Structural changes accompanying the semiconductor–metal transition

_{2}O

_{3}by materials project.

*The Oxide Handbook*

*U*method on polaron properties in energy materials

*U*method for a local-orbital basis

*U*method

*U*study of self-trapping, trapping, and mobility of oxygen-type hole polarons in barium stannate

*U*

_{pp}on oxygen

*p*orbitals in oxides: Role of correlated orbitals on the example of UO

_{2}and TiO

_{2}

_{2}within the LDA+

*U*method

*U*method for correlated materials

*d*-electron systems

*U*energy and its application to actinide compounds

*Pr*

^{3+},

*U*

^{4+}, and

*UPt*

_{3}