Detection and Correction of Delocalization Errors for Electron and Hole Polarons Using Density-Corrected DFT

: Modeling polaron defects is an important aspect of computational materials science, but the description of unpaired spins in density functional theory (DFT) often su ﬀ ers from delocalization error. To diagnose and correct the overdelocalization of spin defects, we report an implementation of density-corrected (DC-)DFT and its analytic energy gradient. In DC-DFT, an exchange-correlation functional is evaluated using a Hartree − Fock density, thus incorporating electron correlation while avoiding self-interaction error. Results for an electron polaron in models of titania and a hole polaron in Al-doped silica demonstrate that geometry optimization with semilocal functionals drives signi ﬁ cant structural distortion, including the elongation of several bonds, such that subsequent single-point calculations with hybrid functionals fail to a ﬀ ord a localized defect even in cases where geometry optimization with the hybrid functional does localize the polaron. This has signi ﬁ cant implications for traditional work ﬂ ows in computational materials science, where semilocal functionals are often used for structure relaxation. DC-DFT calculations provide a mechanism to detect situations where delocalization error is likely to a ﬀ ect the results.

O ne of the most pernicious and long-standing problems in density functional theory (DFT) is that of selfinteraction error (SIE), 1−3 which arises from the incomplete cancellation of Hartree self-repulsion by approximate exchange-correlation (XC) functionals.Among the manifestations of SIE is the overstabilization of fractionally charged systems, 4 leading to an underestimation of reaction barrier heights and motivating the use of hybrid functionals. 5,6−32 Herein, we demonstrate that density-corrected (DC-)DFT can be used both to detect and to correct the overdelocalization of polaron defects.
DC-DFT is a simple procedure whereby a self-consistent Hartree−Fock (HF) density, ρ HF , is used to evaluate the XC energy.This procedure incorporates electron correlation effects but avoids self-consistent iterations at the DFT level, which would introduce SIE into the density, and has been shown to afford reasonable reaction barrier heights even when GGA functionals are used. 33,34−41 In such cases, errors may be ameliorated through the use of an SIE-free density, namely, ρ HF .The usefulness of DC-DFT therefore hinges on identifying cases where the error is densitydriven because otherwise it is probably not advantageous to sacrifice self-consistency. 41,46Delocalization error in systems with one or more unpaired electrons represents a clear case where the error is density-driven, thus we expect such problems to benefit from DC-DFT.
In the present work, we apply this technique to an electron polaron in TiO 2 and a hole polaron in doped SiO 2 .Functionalversus density-driven errors have not previously been considered in these types of systems, yet the open-shell nature of polaron defects suggests that DC-DFT may be advantageous in such cases.(This supposition is supported by results from a density-sensitivity metric, 40 as discussed below.)Within the materials science community there is enormous interest in polaron transport and trapping in transition-metal oxides, 47−51 yet delocalization error in DFT presents a significant conundrum for the computational modeling of these species. 52o introduce the DC-DFT formalism, let us express the energy functional as where the tilde indicates that the functional is approximate in practice because of the approximate nature of E ̃xc [ρ], whereas functional forms for the noninteracting kinetic energy (T s ), the Hartree potential (v H ), and the external potential (v ext ) are known.Error in the approximate energy E ̃may be considered to arise from two contributions: 14 where P μν is the HF density matrix and F μν KS is the Kohn−Sham (KS) Fock matrix constructed from it.Operationally, derivatives ∂P μν /∂x are obtained from molecular orbital (MO) coefficient derivatives, 54 solving the coupled-perturbed equations in Z-vector form for efficiency. 53−56 Although the DC-DFT energy is not variational, its gradient is certainly not "ill-defined" as some have claimed. 57Molecular properties such as multipole moments, polarizabilities, infrared and Raman intensities, and magnetic resonance parameters can each be cast in the form of energy derivatives, 58 and these properties are also well-defined in DC-DFT.
DC-DFT and its analytic gradient have been implemented in the Q-Chem program 59 and will be available in Q-Chem v. 6.0.Validation tests are described in the Supporting Information, including comparisons between analytic and finite-difference gradients (Tables S1−S4) and a demonstration that DC-DFT conserves energy in an ab initio molecular dynamics simulation (Figure S1).We also reproduced previous results 33 demonstrating that DC-BLYP affords improved reaction barriers as compared to either self-consistent BLYP or HF calculations (Table S5).
The remainder of this work describes the application of DC-DFT to polaron defects in oxide materials, where GGAs exhibit significant delocalization error.Periodic DFT is typically the method of choice for modeling solid-state materials but that approach has failed to reproduce correct ground-state electronic and geometric structure in important cases, 25−27,60−63 particularly with respect to the localized nature of O(2p)-derived holes obtained in doped silica and metal oxides.One such example is the polaron in Al-doped silica (SiO 2 ), [24][25][26][27][28]62 where a Si 4+ ion is replaced by Al 3+ .Electron paramagnetic resonance (EPR) spectroscopy suggests that the unpaired spin is localized on a single oxygen atom adjacent to the defect center, 64−66 creating a local structural distortion with a single elongated Al−O bond.
−27 Hyperfine coupling parameters for 17 O, 29 Si, and 27 Al are consistent with experiment when computed at the HF level, 25−27 or alternatively using DFT with a large fraction of exact exchange, 24 and the defect exhibits local symmetry breaking and Jahn−Teller reconstruction. 26Application of semilocal functionals to cluster models, however, affords a hole that is delocalized over all four neighboring oxygen atoms at the Al center, without significant structural distortion. 27,60A localized defect is recovered by applying a self-interaction correction, 28 or by using functionals with a large fraction of exact exchange. 24,26Notably, B3LYP (with 20% exact exchange) fails to localize the hole in Al-doped SiO 2 , 24−26 whereas the functional BB1K (with 42% exact exchange) does localize it, 22,26 as does DFT+U with a large Hubbard U parameter. 62sing the cluster models shown in Figure 1, which are similar to those used in ref 27, we have investigated whether DC-DFT is beneficial in describing the localization of the hole in Al-doped silica.We created defects by replacing the central Si atom with Al and capping the terminal oxygen atoms with hydrogen so that the clusters are charge neutral.All three structures in Figure 1 were initially optimized using secondorder Møller−Plesset perturbation theory (MP2), with the 6-31G* basis set within the resolution-of-identity approximation.
DC-DFT is advantageous only if the density-driven error is large, 38−41 so we examine that question first.A metric −41 In eq 3, this approximate functional is evaluated using the self-consistent density obtained either from a HF calculation or else from the local density approximation (LDA).For small molecules and clusters, values S > 2 kcal/mol have been suggested as cases where DC-DFT may improve the results. 40,41For the structures in Figure 1, values of S computed using either PBE or BLYP are much larger than this threshold, e.g., S > 30 kcal/mol for Al(OH) 4 and S > 120 kcal/mol for Al[OSi(OH) 3 ] 4 (Table S6).This indicates sizable density sensitivity but also suggests that the metric S may increase with system size, at least until the system is large enough for delocalization to reach its asymptotic limit.As such, we also compute S/N atoms as a metric that is more similar The Journal of Physical Chemistry Letters across cluster sizes.This quantity ranges from 2.6 to 7.5 kcal/ mol for the three Al-doped silica models that we consider (Table S6).We suggest that the normalized metric S/N atoms is a more general alternative to S that should be used going forward.
Consistent with experiment, MP2 optimization affords a localized hole (Figure 2a) and a single elongated Al−O bond length (Table 1).Starting from the MP2-optimized structure for each cluster, we next performed geometry optimizations at various DFT levels of theory to investigate the structural distortion and (de)localization of the hole.As indictors of distortion, Al−O bond lengths around the defect center are listed in Table 1 for several different functionals and their DC-DFT analogues.Spin densities (ρ α − ρ β ) for the Al(OSiH 3 ) 4 cluster are shown in Figure 2, and those for the other clusters are shown in Figures S2 and S3.
Cluster structures optimized at the SIE-free HF and MP2 levels of theory exhibit a single elongated Al−O bond, whereas the remaining three Al−O bonds have nearly identical bond lengths.In contrast, geometries obtained with GGA functionals (BLYP and PBE) exhibit four identical Al−O bond lengths and very little structural distortion.This is consistent with a hole that is delocalized over all four bonds, as shown for BLYP in Figure 2b, but is inconsistent with EPR spectroscopy.For hybrid functionals PBE0 and B3LYP, the hole is delocalized over two Al−O bonds rather than four (Figure 2c).This is reflected in two elongated Al−O bonds whose bond lengths have intermediate values between the short and long Al−O bond lengths obtained at the MP2 level (Table 1).
For the hybrid functionals, Al−O distances are omitted from Table 1 for the largest cluster, Al[OSi(OH) 3 ] 4 , because the hole migrates to the surface of the cluster as the geometry is relaxed, ultimately delocalizing over the surface oxygen atoms around the Si centers (Figure S3).−26 In contrast, upon relaxing the structure of each doped silica cluster using DC-DFT, starting from the localized MP2 structure, localization is preserved in all cases: DC-BLYP (Figure 2e), DC-B3LYP (Figure 2f), DC-PBE, and DC-PBE0.(Note that the DC-DFT spin densities in Figure 2 are HF densities because plotting DFT spin densities would be inconsistent with the strategy of DC-DFT, which is to avoid introducing SIE by avoiding the DFT density altogether.)Consistent with the localized spin densities obtained with DC-DFT, the bond lengths in Table 1 show that all four DC-DFT methods afford a single elongated Al−O bond, with bond lengths that are within 0.03 Å (or less) of MP2 results, for both the shorter and longer Al−O bonds.
DC-DFT geometry optimizations afford a properly localized hole in these doped silica clusters precisely because the energy and gradient avoid the use of a DFT density.To emphasize this point, Table 2 reports Mulliken "spin charges" (obtained from ρ α − ρ β ) for the oxygen atoms around the Al defect, computed using either the HF density or else the BLYP density, at the same optimized geometry of the cluster in either case.For the DC-DFT cases in Table 2, this means that we perform the geometry optimization using ρ HF but then evaluate ρ BLYP at the optimized geometry, for demonstrative purposes.(Results based on ρ HF represent the proper DC-DFT procedure.)Despite the limitations of the Mulliken procedure, these spin charges have proven to be an effective way to characterize the (de)localization of a radical center, 16−18 and a similar picture is obtained with the larger def2-TZVPD basis set (Table S8).Spin delocalization moves the defect toward the surface oxygen atoms (Figure S3).

The Journal of Physical Chemistry Letters
Whereas DC-DFT affords a well-localized hole with a single spin charge of ≈1, the BLYP density evaluated at the same geometry reintroduces delocalization of the hole, reducing the spin charge to ≈0.7 in Al(OH) 4 and to ≈0.6 in Al(OSiH 3 ) 4 .The hole does not completely delocalize because the HF-or DC-DFT-optimized geometry elongates one Al−O bond, creating a bond dipole that is evidently sufficient to partially localize the spin density even in BLYP or PBE calculations.When one of those GGA functionals is used to optimize the geometry, however, delocalization becomes complete and results in four equivalent spin charges on the oxygen atoms surrounding the Al defect.For the largest cluster, Al[OSi-(OH) 3 ] 4 , GGA calculations at GGA-optimized geometries afford a hole that is completely delocalized over the entire cluster.Conversely, geometries optimized using either B3LYP or PBE0 result in two elongated Al−O bonds.In these cases, the two exaggerated bond dipoles are large enough that a subsequent HF calculation cannot completely localize the hole, and instead two significant spin charges are obtained, with values of ≈0.6 in Al(OH) 4 and ≈0.5 in Al(OSiH 3 ) 4 .
Whereas calculations reported up to this point use a dopant to create a hole, in transition-metal oxides it is possible to create localized electronic states via structural distortion following excitation at the band gap, 50,51,67,68 with important applications in photocatalysis. 69,70For simplicity, we create electron polarons in the ground state of model systems representing titania by introducing an extra electron into a closed-shell system, and we find that many of the same trends are on display as were observed with dopant-induced holes.For example, consider the small (TiO 2 ) 3 (H 2 O) 3 cluster that is depicted in Figure 3a.By breaking the symmetry of the ring and introducing an unpaired electron, a spin-localized structure for [(TiO 2 ) 3 (H 2 O) 3 ] − can be optimized at the HF level, and this localization is preserved upon structural relaxation using DC-PBE (Figure 3a).Relaxation at the PBE level, however, leads to significant spin delocalization, with Mulliken spin charges s 1 = 0.64 and s 2 = 0.24 = s 3 on the three Ti atoms.PBE0 optimization leads to a much smaller amount of delocalization (s 1 = 0.95 and s 2 = 0.07 = s 3 ).
Figure 3b examines the energetics with various approaches, computed along the DC-PBE geometry optimization pathway starting from the spin-localized HF structure.Nearly identical energetics to those obtained at the DC-PBE level (used to optimize the structure) are obtained when the energy is evaluated at each step along the pathway using a self-consistent PBE calculation.Furthermore, the spin remains largely localized on a single Ti atom despite the fact that a selfconsistent PBE geometry optimization produces a much more delocalized electronic structure.This suggests that SIE drives geometric distortions that, in turn, facilitate additional delocalization, creating a feedback loop when the GGA functional is employed self-consistently.
−72 Understanding the behavior of polaronic states is crucial in understanding charge transport in this and other metal oxide photocatalysts, 47−52 but the degree to which an electron polaron localizes in anatase TiO 2 has been a subject of debate.Some studies predict delocalization of the excess charge over multiple Ti ions, 73−75 while others suggest the formation of a localized or "self-trapped" polaron. 68t has been suggested that localized and delocalized states may be similar in energy, 76 and that the localization properties may vary from one crystal polymorph to another. 74Semilocal functionals overstabilize charge-delocalized structures, 29 and while this propensity is somewhat reduced by employing hybrid functionals, 30,31 or alternatively via DFT+U, 75−80 the use of potentially nontransferrable parameters (such as the fraction of exact exchange or the Hubbard U parameter) means that some experimental data are needed to determine the correct approach. 77However, EPR experiments have produced conflicting evidence regarding the (de)localization of electron polarons in anatase TiO 2 . 73e constructed a cluster model of anatase TiO 2 as follows.Starting from the crystal structure of the bulk material, we extracted a cluster Ti 21 O 70 H 56 that is depicted in Figure 4a, capping the surface oxygen atoms with hydrogen to maintain charge neutrality.This has been suggested to be the smallest cluster model where the central Ti atom retains bulk TiO 2 character. 81The geometry of this cluster model was then relaxed using HF theory but constraining the surface Ti−O and O−H bond lengths to provide some mechanical stability at the surface.Following ref 81, these optimizations used the  S3).

The Journal of Physical Chemistry Letters
LANL2DZ basis set along with the eponymous effective core potential for Ti. 82 model for an electron polaron is obtained by adding an extra electron to the Ti 21 O 70 H 56 cluster model in a manner such that the polaron is associated with the central (bulk-like) Ti atom.Naive approaches to adding this electron cause it to delocalize at the surface of the cluster, which is an edge effect rather than an SIE artifact, as evidenced by problematic delocalization even at the HF level.Various strategies have been devised to obtain a localized polaron in a finite cluster model, 83 typically involving either bond distortion or else changes to the nuclear charges.Following other theoretical studies of TiO 2 , in which Ti−O bond lengths were elongated by 2−8% in order to localize the defect, 31,84 we opted to elongate all six Ti−O bonds around the central Ti atom in the initial geometry of [Ti 21 O 70 H 56 ] − .The elongated Ti−O bond lengths at the initial geometry can be seen in Figure S4c, and the structure was subsequently relaxed at the HF level using the Mulliken spin charge s(Ti*) on the central atom (Ti*) to monitor (de)localization as the optimization proceeds.This analysis demonstrates that the localized polaron that is initially induced by bond elongation remains localized when the structure is relaxed at the HF level (Figure 4b).The density sensitivity at the HF-optimized polaron geometry is S/N atoms ≈ 8 kcal/mol using either BLYP or PBE, which is even larger than was found for the Al-doped silica clusters.
Spin charges computed at the PBE0 level along the HF optimization trajectory (Figure S4b) also correspond to a localized defect.In contrast, the PBE functional partially delocalizes the polaron even when the HF structure is used, with s(Ti*) < 0.6 at the HF-optimized geometry as compared to values s(Ti*) ≈ 1 that are obtained using either HF or PBE0.This is the same behavior that was observed for Aldoped silica.In the HF-optimized structure of [Ti 21 O 70 H 56 ] − , all of the Ti*−O bonds are longer than their values in the neutral structure, but two of them exhibit significantly greater elongation, suggesting the participation of two Ti−O bonds (Figure S4c).
We next investigate whether this level of localization persists with DFT.Starting from the HF-optimized structure of [Ti 21 O 70 H 56 ] − , we reoptimized the geometry using either PBE or PBE0 and monitored s(Ti*) as a proxy for polaron localization.The results are plotted in Figure 5. Bond-length analysis along the PBE0 optimization pathway (Figure S5a) suggests a localized polaron, and a plot of s(Ti*) along the same pathway (Figure 5a) demonstrates that the extra electron indeed remains localized on the central atom, Ti*.If s(Ti*) is computed instead using the PBE functional, along the same PBE0 optimization pathway, then the polaron seems stable but a significant fraction of the spin density is transferred onto other atoms, resulting in s(Ti*) < 0.6.
When the optimization is performed instead with PBE, then not only is the PBE polaron signif icantly delocalized, with s(Ti*) < 0.2, but also the spin charge computed using PBE0 indicates that PBE-driven changes in bond lengths cause delocalization that persists at the PBE0 level.The value of s(Ti*) starts at ≈1 when computed with PBE0 but falls to <0.5 along the PBE optimization pathway, with significant fluctuations in between (Figure 5b).When the HF density is evaluated along this PBE optimization pathway, Mulliken charges indicate hopping of the polaron between Ti atoms.This behavior likely arises from a combination of HF theory's strong tendency to localize charge, combined with a GGA optimization that is driving the elongation of multiple Ti−O bonds.At any given step along the pathway, the HF calculation localizes the polaron on whichever Ti−O bond is longest, and that location switches several times over the course of the optimization (Figure S5b).
These observations are significant in view of a standard operating procedure in computational materials science, which is to perform geometry optimizations using GGA functionals  The Journal of Physical Chemistry Letters and then to use hybrid functionals for single-point calculations at GGA-optimized geometries.(This procedure mitigates the high cost of hybrid DFT under periodic boundary conditions.)Our results demonstrate that PBE optimization drives the system into structures that are severely affected by SIE, with elongated bonds (consistent with a fully delocalized spin defect) from which a PBE0 single-point calculation cannot recover.A PBE0 calculation at the PBEoptimized geometry thus affords a delocalized polaron, which is a qualitatively different result as compared to what is obtained when PBE0 is used for the optimization.
We next evaluate the performance of DC-DFT along the PBE and PBE0 optimization pathways.The PBE pathway is more interesting insofar as a major change in polaron localization is observed in this case, as evident from the spin charges in Figure 5b.We seek to assess the extent to which this density-driven error is mitigated by DC-DFT.This assessment requires a benchmark that is SIE-free, for which we choose MP2/def2-SVP.The treatment of dynamical correlation may help to assuage a tendency toward overlocalization that is exhibited by the HF method, which is evident from its behavior for fractional electron systems. 85−90 Figure 6 reports the results of single-point energy calculations along the PBE and PBE0 optimization pathways.The starting point is a HF-optimized structure, hence the HF energy immediately rises as either a PBE-or PBE0-driven optimization moves the system away from the minimumenergy HF geometry.Notable discontinuities in the HF energy along the PBE pathway (Figure 6a) arise when the localized HF polaron jumps between Ti centers, as a result of HF theory's strong tendency to localize unpaired spins, combined with the simultaneous distortion of Ti−O bonds due to PBE delocalization error.In Figure S8, one can observe that the HF polaron is localized on a different Ti atom in the starting and ending structures optimized at the PBE level.Jumps in the HF energy in Figure 6a correlate with jumps in the gap between the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO), computed at the HF level (Figure S9a).These gaps are absent at the PBE level, consistent with the wellknown band gap problem associated with GGA functionals, 91−93 which reflects the absence of a proper derivative discontinuity and is thus intimately connected to SIE and delocalization error. 92A Born−Oppenheimer molecular dynamics simulation using a GGA functional might therefore be anticipated to afford incorrect polaron dynamics  The Journal of Physical Chemistry Letters changes in the electronic configuration are rendered thermally accessible by artificially small energy gaps.
Even a modest fraction of exact exchange introduces something of a derivative discontinuity, and the energetics at the PBE0 level of theory closely track the MP2 results.The MP2, hybrid DFT, and DC-DFT methods do jumps in energy when the polaron localizes on a different site, although MP2 substantially mutes this effect as compared to HF theory.These jumps are essentially artifacts of the GGA-driven optimization and are absent when either PBE0 or HSE06 is used to optimize the geometry (Figures 6b and S10a, respectively).Note also that in the latter two cases the PBE functional lowers the energy even faster than the hybrid functional that is being used to optimize the geometry!This is an SIE artifact (overstabilization of fractional charges), and PBE predicts a delocalized polaron even at the PBE0optimized geometry, where PBE0 itself localizes the polaron (Figure 5a).Overall, the close match between PBE0 and MP2 suggests that the former is a good level of theory for this system, although that fact is evident only in light of the MP2 results.
When DC-DFT is applied to this system, energetics for the PBE optimization pathway are shifted upward by about 0.1 hartree, placing DC-PBE in reasonable agreement with MP2.The same is true for DC-PBE along the PBE0 optimization pathway, only in that case the DC-PBE energy profile is free of the jumps that stem from hopping of the localized HF spin density between Ti sites.At least in this system, DC-PBE appears to offer an alternative means to achieve MP2-quality results at HF cost.PBE0 already affords accurate energetics (comparable to MP2) along the PBE, PBE0, and HSE06 optimization pathways, but DC-PBE0 shifts these energetics upward by 0.02−0.03hartree.The use of DC-PBE0 therefore does not offer any advantage over PBE0 itself because both methods require the computational equivalent of a HF calculation.
Of the ∼0.2 hartree of energy lowering that is obtained at the PBE level when the HF-optimized starting structure is relaxed (Figure 6a), comparison to the DC-PBE result suggests that about half of that amount is density-driven error.In the PBE0 optimization, the difference between PBE and DC-PBE is much smaller, consistent with the idea that the densitydriven error is reduced when a hybrid functional is used for the optimization.Shifts between DFT and DC-DFT are plotted in Figure S11 for two functionals and two optimization pathways, demonstrating that the shift for PBE0 is consistently smaller than that for PBE, as a result of the smaller density-driven error in the hybrid functional.
The optimization pathway obtained using the HSE06 functional 94,95 affords results that are very similar to PBE0driven optimizations.This is not entirely surprising, given that HSE06 is a short-range hybrid designed to mimic PBE0 while reducing cost under periodic boundary conditions.The SCAN functional 96 affords energetics that are close to MP2 along various optimization pathways (Figure S12), but that means that DC-SCAN shifts the energies upward relative to MP2.The same is true for the hybrid SCAN0 functional. 97The DC-SCAN procedure has been recommended in some recent work, 43−45 but in this particular application the quality of SCAN energetics is somewhat degraded by density correction.
In summary, we have implemented DC-DFT along with its analytic gradient in the Q-Chem program and have assessed this technique as a means to describe the electronic and geometrical structure of model systems containing polaron defects.In Al-doped silica, where EPR measurements suggest that the dopant-induced hole is localized on a single oxygen atom, GGA functionals and hybrids such as B3LYP and PBE0 instead delocalize the hole over more than one oxygen atom.DC-DFT versions of these same functionals afford a localized hole and structural distortion that is consistent with experiment, i.e., a single elongated Al−O bond.This suggests that DC-DFT might be an alternative to hybrid functionals in these systems.
As a second application, we studied localization of an extra electron in a sizable cluster model of anatase TiO 2 .EPR spectroscopy of electron polarons in this system has certain features that are thought to be consistent with a localized Ti 3+ defect but has other features consistent with delocalization and a "large polaron". 73Starting from a HF-optimized structure with a localized polaron, we find that optimization with PBE0 or HSE06 preserves the defect localization whereas optimization with PBE drives delocalization onto multiple Ti atoms, leading to the elongation of multiple Ti−O bonds.
In both titania 2 and Al-doped silica, GGA-induced structural distortion is severe enough that a subsequent PBE0 singlepoint calculation cannot recover a localized polaron because exaggerated bond dipoles over multiple stretched bonds provide a driving force for delocalization and cause fractionation of the spin density.This is a significant observation insofar as a standard procedure in computational materials science is to perform structure relaxations using a GGA (often PBE) followed by single-point energy calculations with a hybrid functional (often PBE0 or HSE06).That procedure affords qualitatively incorrect electronic structure in several of the examples considered here.
For relaxation of the polaron in TiO 2 , DC-PBE affords energetics in good agreement with MP2.For the PBE0 and SCAN functionals, however, DC-DFT worsens what had been favorable agreement with the MP2 results.Because the cost of DC-DFT in periodic calculations will be similar to the cost of plain PBE0, we suggest that the role of DC-DFT should be a diagnostic one: where SIE problems are anticipated, DC-DFT provides an indication of whether defect delocalization is likely to persist as the fraction of exact exchange is adjusted.It does so without adjustable parameters and without the need to abandon electron correlation.The DC-DFT analytic gradient, missing in widely used quantum chemistry programs until now, 41,42 affords the ability to describe structural relaxation with this method and thus to interrogate how the elimination of SIE impacts the geometric structure of molecules and materials.
Computational details and additional analysis (PDF) Transparent Peer Review report available (PDF) ■ AUTHOR INFORMATION a functional-driven error E ̃xc [ρ] − E xc [ρ] that originates in the approximate nature of the XC functional, where ρ is the exact density obtained by minimizing the exact functional E[ρ], and a density-driven error E ̃[ρ] − E ̃[ρ], where ρ̃indicates the approximate density obtained by self-consistent variational minimization of eq 1. DC-DFT reduces the density-driven error by eliminating ρ̃in favor of ρ HF .The DC-DFT energy expression, E DC−DFT = E ̃[ρ HF ], is not variational and therefore its analytic gradient requires the solution of coupled-perturbed equations.53,54Following the formalism in ref 33, the DC-DFT analytic energy gradient is

Figure 1 .
Figure 1.Cluster models of Al-doped silica, where the central Si atom of the SiO 2 crystal structure is replaced by Al, which is the central atom in each cluster.

Figure 2 .
Figure 2. Spin densities for Al(OSiH 3 ) 4 using geometries that have been optimized at the indicated levels of theory.For DC-BLYP and DC-B3LYP, the HF spin density is plotted, and for MP2 we plot the "relaxed" (correlated) one-electron spin density.Spin densities computed at PBE-and PBE0-optimized structures are not shown because they are essentially indistinguishable from the BLYP and B3LYP results, respectively.The 6-31G* basis is used for MP2 and 6-31++G* is used for the other methods.

Figure 3 .
Figure 3. (a) Spin densities [(TiO ) 3 (H 2 ] − , with positive and negative values of ρ α − ρ β indicated in purple and green, respectively.The central structure emphasizes the Ti atoms as larger gray spheres, and around the periphery are structures optimized at the indicated levels of theory.(The LANL2DZ basis set and pseudopotential are used in each case.)The quantity s 1 is the Mulliken spin charge on the Ti atom at the apex of the ring.(b) Relative energies ΔE (open circles, to be read from the scale on the left) and spin charges s 1 (filled squares, scale on the right) along the DC-PBE relaxation pathway starting from the HF-optimized structure.

Figure 4 .
Figure 4. Cluster model of an electron polaron in anatase TiO 2 .(a) Capped cluster Ti 21 O 70 H 56 extracted from the bulk TiO 2 crystal structure.(b) HF spin density for a localized polaron in [Ti 21 O 70 H 56 ] − .

Figure 5 .
Figure 5. Variation of the spin charge s on the central Ti atom (Ti*) in [Ti 21 O 70 H 56 ] − , computed at different levels of theory as the structure is optimized using (a) PBE0 versus (b) PBE.The starting structure in either case is optimized at the HF level and exhibits a localized polaron.Insets on the right show the spin densities at various ending points.In (b), the HF result for s(Ti*) is not shown because the spin density hops between being localized on different Ti atoms as the optimization proceeds.

Figure 6 .
Figure 6.Relative along (a) PBE and (b) PBE0 optimization pathways for an electron polaron in anatase TiO 2 , starting from a HFoptimized geometry that localizes the polaron on the central Ti atom.Single-point energy calculations (def2-SVP basis) are performed using each of the methods indicated.All energies are shifted to a common value at the initial structure.

Table 1 .
Optimized Al−O Bond Distances (in Å) in Al-Doped Silica Models