Inapplicability of exact constraints and a minimal two-parameter generalization to the DFT+ U based correction of self-interaction error
IInapplicability of exact constraints and a minimal two-parameter generalization to theDFT+ U based correction of self-interaction error Glenn Moynihan, ∗ Gilberto Teobaldi,
2, 3 and David D. O’Regan School of Physics, CRANN and AMBER, Trinity College Dublin, Dublin 2, Ireland Stephenson Institute for Renewable Energy and Department of Chemistry,The University of Liverpool, L69 3BX Liverpool, United Kingdom Beijing Computational Science Research Center, Beijing 100094, China (Dated: December 15, 2016)In approximate density functional theory (DFT), the self-interaction error is an electron delo-calization anomaly associated with underestimated insulating gaps. It exhibits a predominantlyquadratic energy-density curve that is amenable to correction using efficient, constraint-resemblingmethods such as DFT + Hubbard U (DFT+ U ). Constrained DFT (cDFT) enforces conditions onDFT exactly, by means of self-consistently optimized Lagrange multipliers, and while its use toautomate error corrections is a compelling possibility, we show that it is limited by a fundamentalincompatibility with constraints beyond linear order. We circumvent this problem by utilizing sep-arate linear and quadratic correction terms, which may be interpreted either as distinct constraints,each with its own Hubbard U type Lagrange multiplier, or as the components of a generalizedDFT+ U functional. The latter approach prevails in our tests on a model one-electron system, H +2 ,in that it readily recovers the exact total-energy while symmetry-preserving pure constraints failto do so. The generalized DFT+ U functional moreover enables the simultaneous correction of thetotal-energy and ionization potential, or the correction of either together with the enforcement ofKoopmans’ condition. For the latter case, we outline a practical, approximate scheme by which therequired pair of Hubbard parameters, denoted as U and U , may be calculated from first-principles. PACS numbers: 71.15.-m, 31.15.E-, 71.15.Qe, 71.15.Dx
Approximate density-functional theory (DFT) un-derlies much of contemporary quantum-mechanicalatomistic simulation, providing a widespread and valu-able complement to experiment . The predictive capac-ity of DFT is severely limited, however, by systematic er-rors exhibited by tractable exchange-correlation func-tionals such as the local-density approximation (LDA) or generalized-gradient approximations . Perhaps themost prominent of these pathologies is the delocalizationor many-electron self-interaction error (SIE) , which isdue to spuriously curved rather than correctly piecewise-linear total-energy profiles with respect to the total elec-tron number . This gives rise to the underestimatedfundamental gaps , charge-transfer energies , and re-action barriers characteristic of practical DFT. Whilethe construction of viable, explicit density functionalsfree of pathologies such as SIE is extremely challeng-ing , significant progress has been made in the devel-opment of implicit functionals in the guise of correctiveapproaches. Examples include methods that operate bycorrecting SIE on a one-electron basis according to vari-ationally optimized definitions, such as generalized Perdew-Zunger approaches, in which much progress hasrecently been made by generalizing to complex-valued or-bitals , and those which address many-electron SIEdirectly, such as Koopman’s compliant functionals .An established, computationally very efficient DFTcorrection scheme is DFT+ U , originally developedto restore the Mott-Hubbard effects absent in the LDAdescription of transition-metal oxides. A simplified for-mulation , in which the required Hubbard U param- eter is a linear-response property of the system underscrutiny , is now routinely and diversely applied .Beginning with Ref. 11, Marzari, Kulik, and co-workershave suggested and extensively developed the inter-pretation of DFT+ U as a correction for SIE, for sys-tems in which it may be primarily attributed to dis-tinct subspaces (otherwise, the related Koopman’s com-pliant functionals are available ). The SIE correctingDFT+ U functional is given, where ˆ n Iσ = ˆ P I ˆ ρ σ ˆ P I , by E U = (cid:88) Iσ U I (cid:2) ˆ n Iσ − ˆ n Iσ ˆ n Iσ (cid:3) . (1)Here, ˆ ρ σ is the Kohn-Sham density-matrix for spin σ and ˆ P I is a projection operator for the subspace I .DFT+ U attains the status of an automatable, first-principles method when it is provided with calculatedHubbard U parameters (particularly at theirself-consistency ), which may be thought of assubspace-averaged SIEs quantified in situ. The subspacesare usually pre-defined for corrective treatment, havingbeen deemed responsible for the dominant SIEs on thebasis of physical intuition and experience, although a fur-ther level of self-consistency over subspaces is possibleusing Wannier functions . DFT+ U effectively adds aset of penalty functionals promoting integer eigenvaluesin ˆ n Iσ , and it replicates the effect of a derivative discon-tinuity in the energy, for each subspace I , by adding anoccupancy-dependent potential ˆ v Iσ = U I ( ˆ P I / − ˆ n Iσ ).While DFT+ U is effective and computationally effi-cient, even linear-scaling , a considerable degree of careis needed to calculate the required U parameters, which a r X i v : . [ phy s i c s . c h e m - ph ] D ec sometimes pose numerical challenges . Fully self-contained calculations of the Hubbard U by means ofautomated variational extremization would be extremelyuseful for many practitioners, and expedient in high-throughput materials search contexts . The constraint-like functional form of DFT+ U , where the U I resem-ble the Lagrange multipliers of penalty functionals onthe eigenvalues of ˆ n Iσ , suggests the possible viability ofsuch a method. Constrained density-functional theory(cDFT) formalizes and automates the use of self-consistent penalty functionals in DFT, enforcing themas exact constraints by locating the lowest-energy com-patible excited state of the underlying functional. It is ef-fective for treating SIE, in its own right, for systems com-prising well-separated fragments, where it may be usedto break physical symmetries and to explore the integer-occupancy states at which SIE is typically reduced .However, as we now demonstrate, cDFT is fundamentallyincompatible with constraints beyond linear order, andtherefore exact constraints cannot be used to correct SIEin an automated fashion. As a result, it seems that wecannot excite a SIE affected system to a state that willreliably exhibit less SIE, without breaking a symmetry.The simplest conceivable SIE-targeting constraintfunctional is the quadratic form C = (cid:80) I ( N I − N I c ) ,where N I = Tr[ˆ n I ] is the total occupancy of a particu-larly error-prone subspace I and N I c is its targeted value,neglecting the spin index for concision. This constraintis a functional of subspace total occupancies, rather thanthe occupancy eigenvalues as in DFT+ U , which is animportant distinction for all but single-orbital sites. For N sites symmetry-equivalent subspaces with N I c = N c for all I , the total-energy of the system is given by W = E DFT + V c C , where V c is a common cDFT La-grange multiplier. This gives rise to a constraining po-tential of the form ˆ v c = 2 V c (cid:80) I ( N I − N c ) ˆ P I , making ex-plicit its dependence on the constraint non-satisfaction.This, in turn, implies an externally imposed interactioncorrection given by ˆ f c = 2 V c (cid:80) I ˆ P I ˆ P I , which acts tomodify the energy-density profile, and which is identicalto that generated by DFT+ U (Eq. 1) when V c = − U I / dW/dV c = C n , so that the total-energy W ( V c ) alwaysattains a stationary point upon constraint satisfaction,in this case when C = 0. Fig. 1 depicts this functionfor an ideal system for the study of one-electron SIE , H +2 , simulated using the PBE functional . At the con-sidered, intermediate bond-length of 4 a , the overlap ofthe two atom-centered PBE 1 s orbital subspaces yields atotal occupancy double-counting of 24%, accounting forspillage. The observed asymptotic behavior of W ( V c )demonstrates that the C constraint is unenforceable.Here, a target occupancy of N c = 0 . V c , butthe same qualitative outcome arises for any N c (cid:54) = N DFT .The key to the failure of the C constraint is the ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● - - - - - - V c ( eV ) T o t a l E n e r gy ( e V ) Energy WC ( V c ) dC / dV c d C / dV c2 - - - - - - | d n C / d V c n ( e V - n ) | FIG. 1. (Color online) The constrained total-energy of thestretched H +2 system, with a target occupancy of N c =0 . s -orbital subspace, against thecDFT Lagrange multiplier V c . Also shown is the constraintfunctional C = ( N − N c ) , averaged over the two atoms, andits first and second derivatives, which all fall off rapidly with V c . fall-off of all self-consistent cDFT response functions d n N/dV n c = d n +1 W/dV n +1 c , as depicted in Fig. 2. Thisresults in a diminishing returns as V c is increased, i.e.,as the constraint asymptotically approaches satisfaction.Motivated by this, we investigate whether non-linear con-straints of form C n = ( N − N c ) n are unsatisfiable for anyorder n (cid:54) = 1 and for any target choice N c (cid:54) = N DFT . The n = 0 case is trivial, and the constraint is ill-definedfor n < C n becomes imaginary for non-integer n with negative ( N − N c ), so that we may limitour discussion to integers n ≥
2. We begin by analyzingthe derivatives of the total-energy W ( V c ). The secondderivative follows directly from the above discussion, as d WdV c = dC n dV c = n ( N − N c ) n − dNdV c . (2)The energy derivative of order m generally involves cDFTresponse functions up to order m −
1, and positive inte-ger powers of ( N − N c ) which may vanish, dependingon m and n , but not diverge. The cDFT response func-tion dN/dV c may be gainfully expanded, if ˆ v ext is theexternal potential, in terms of the intrinsic subspace-projected interacting response function defined by χ =Tr[( dN/d ˆ v ext ) ˆ P ], since this object is independent of theform of the constraint. The first-order cDFT response dN/dV c is thus expressed, by means of the chain rule inˆ v ext = ˆ v c = V c ( δC n /δ ˆ ρ ) = nV c ( N − N c ) n − ˆ P , as dNdV c = Tr (cid:20) dNd ˆ v ext d ˆ v ext dV c (cid:21) = nχ ddV c [ V c C n − ] ⇒ dNdV c = nχC n − (cid:16) − n ( n − χV c C n − (cid:17) − , (3)an expression which we have verified numerically. At anyvalid stationary point, C n = 0 and each of V c , χ , and itsderivatives must remain finite. Thus, for n ≥
2, theresponse dN/dV c and energy curvature d W/dV c boththen vanish. The latter is therefore not a stationary pointdiscriminant, and we move to higher derivatives, such as d WdV c = n ( n − C n − (cid:18) dNdV c (cid:19) + nC n − d NdV c , (4)where the required second-order response is given by d NdV c = n (cid:20) dχdV c C n − + ( n − dNdV c × (cid:18) dχdV c V c C n − + 2 χC n − + ( n − χV c C max( n − , dNdV c (cid:19)(cid:21) × (cid:16) − n ( n − χV c C n − (cid:17) − . (5)This object, and thus d W/dV c both vanish at stationarypoints for all n ≥
2, due to the vanishing C n − in thefirst term of Eq. 5, and due to the vanishing first-orderresponse (for which, see Eq. 3) in all remaining terms.In general, the cDFT response function d m N/dV m c comprises terms proportional to response functions of thesame type but of lower order, plus a single term whichis proportional to a potentially non-vanishing mixed re-sponse function d m − χ/dV m − c multiplied by the nec-essarily vanishing C n − . This serves as an inductiveproof that response functions at all orders, beginningwith dN/dV c , vanish as we approach a vanishing C n , asillustrated in Fig. 2. Then, since each term in the m th energy derivative is always proportional to non-divergentpowers of ( N − N c ) and response functions of at most or-der m −
1, all energy derivatives tend to zero, as depictedin Fig. 1, proving the conjecture. Thus, non-linear con-straints of SIE-targeting C n form cannot be enforced.In order to cast the SIE-targeting C n functionalinto a viable form, one possible option remains. Wemay expand the single-site C , for example, as C = − N c V c ( N − N c ) − V c ( N c − N ), and afford an additionaldegree of freedom to the system by decoupling these twoterms. Writing the result in the notation of DFT+ U , bychange of variables, we arrive at the constraint energy (cid:88) I U (cid:0) N I − N c (cid:1) + (cid:88) I U (cid:0) N c − N I (cid:1) . (6)The vanishing response problem is now circumvented,by interpreting the Hubbard U parameters for linearand quadratic terms as separate Lagrange multipliers.Adapting Eq. 6 to multiple, multi-orbital sites and ne-glecting inter-eigenvalue terms, in the spirit of DFT+ U ,while retaining only the free-energy (setting N I c = 0),we arrive at the generalized DFT+ U correction given by E U U = (cid:88) Iσ U I (cid:2) ˆ n Iσ (cid:3) − (cid:88) Iσ U I (cid:2) ˆ n Iσ (cid:3) . (7) χ ( V c ) dN / dV c d N / dV c2 d N / dV c3 - - - - - Lagrange Multiplier V c ( eV ) | R e s pon s e F un c ti on | ( e V - n ) FIG. 2. (Color online) The magnitudes of the interactingdensity response χ and cDFT response functions d m N/dV m c ,calculated from a polynomial fit to the average subspace oc-cupancy for the same system as in Fig. 1. The d m N/dV m c fall off as we asymptotically approach constraint satisfaction,while the occupancy, and hence χ , tends to a constant value. Here, the DFT+ U functional of Eq. 1 is recovered bysetting U I = U I . Otherwise, the corrective potential ismodified to ˆ v IU U = U I ˆ P / − U I ˆ n I , so that the charac-teristic occupancy eigenvalue dividing an attractive froma repulsive potential is changed from 1 / U I / U I .Self-consistency effects aside, the U parameters are re-sponsible for correcting the interaction and for any gapmodification, while the U parameters may be used toadjust the linear dependence of the energy on the sub-space occupancies, and thereby to refine eigenvalue de-rived properties such as the ionization potential. We notea resemblance between Eq. 7 and the three-parameterDFT+ U αβ functional proposed in Ref. 45, where here athird degree of freedom may be retained by using N I c (cid:54) = 0.Fig. 3 shows the total-energy W of the H +2 system asbefore, but now against the U and U defined in Eq. 6,with a subspace target occupancy of N c = N exact =0 .
602 e, (the population of each of the two PBE 1 s or-bital subspaces, calculated using the exact functional).The total-energy is non-uniquely maximized along theheavy white line where the constraint is satisfied, at ∼ .
92 eV below the exact energy. To understand whythe total-energy is always degenerate, and hence whythe occupancy condition under-defines the pair ( U , U ),it suffices to show that the Hessian of the constraintfunctional, H ij = d W/dU i dU j , is everywhere sin-gular. The determinant of H ij is conveniently calcu-lated in terms of the response functions, i.e., by us-ing the ground-state expressions dW/dU = N/ dW/dU = − N /
2, as | H | = 12 (cid:12)(cid:12)(cid:12)(cid:12) dN/dU − dN /dU dN/dU − dN /dU (cid:12)(cid:12)(cid:12)(cid:12) = 0 , (8)as required, for all U and U . This implies a vanish- U ( eV ) Q u a d r a ti c U ( e V ) - - - - FIG. 3. (Color online) The constrained total-energy ofstretched H +2 against the Lagrange multipliers U and U defined in Eq. 6. The subspace target occupancy is set to N c = N exact , and the zero is set to the exact total-energy. Theconstraint is satisfied at the total-energy maximum along thesolid white line. The ionization potential is exact along thethick dashed line. The linear-response Hubbard U , togetherwith the U value needed to correspondingly recover the exactsubspace density, are shown using thin dashed lines. ing energy curvature along the lines on which the cor-rective potential is constant. The linear-response Hub-bard U at this bond length, calculated using a methodadapted from Ref. 25, is 4 .
84 eV. If we intuitively set U = U , then a corresponding U = 3 .
16 eV is re-quired to recover the exact subspace density. The lineon which the ionization potential is exact (for the specialcase of H +2 , this means that the occupied Kohn-Shamstate eigenvalue and the ion-ion energy sum to the ex-act total-energy, written ε DFT + E ion-ion = E exact ), inter-cepts U ≈ U = U , echoingthe ‘SIC’ double-counting correction proposed in Ref. 47.Finally, for the constrained total-energy, we note thatwhile it can be tuned to reach a maximum at the exactenergy for a plausible target, N c = 0 .
511 e, an unrea-sonable U = U = − . U term of Eq. 7, the total-energy generated bywhich is shown in Fig. 4. The zero of energy and the U ( eV ) Q u a d r a ti c U ( e V ) - - - FIG. 4. (Color online) As per Fig. 3 but showing the free-energy obtained by setting N c = 0 e in Eq. 6, i.e., using thegeneralized DFT+ U of Eq. 7. The thick dashed line is theexact energy intercept, and the thin dashed lines show thelinear-response U together with the corresponding U neededto recover the exact energy. The solid white line, as in Fig. 3,indicates where the exact subspace occupancy is recovered. heavy dashed line show E exact , and the thin dashed linesindicate the Hubbard U = 4 .
84 eV and corresponding U = 4 .
44 eV required to recover it. E exact is attained bya traditional DFT+ U calculation at U = U = 3 .
85 eV,and the intersection of the heavy solid and dashed linesyields the pair, U = 5 .
73 eV and U = 6 .
98 eV, at which E exact and N exact are located. At the same point, theKohn-Sham eigenvalue ε DFT lies at ∼ . ε exact ,reflecting that an accurate total-energy at a particularoccupancy may coincide with an inaccurate ionizationenergy, and vice versa, as was recently shown in detailedanalyses of the residual SIE in hybrid functionals andin DFT+ U itself, at fractional total occupancies .The generalized DFT+ U functional enables simulta-neous correction of the ionization potential and total-energy, or the correction of either together with Koop-mans’ condition . In the one-electron case, as in H +2 ,Koopmans’ condition may be enforced at ε exact . Theapproximately parallel solid curves of Fig. 5 illustratethe large U and U values required to do so, as a func-tion of internuclear distance. To estimate these, wehave used the convenient feature of H +2 that the PBE1 s orbital subspace projectors closely match Kohn-Shamorbitals in spatial profile, almost exactly so at disso-ciation. As demonstrated by the occupancy curves inFig. 5, this implies a negligible charge and kinetic self- U K / UN ( U = ) N ( U = ) U K / U a nd A v e r a g e O cc up a n c y U U ( a ) H ubb a r d P a r a m e t e r ( e V ) FIG. 5. (Color online) Generalized Hubbard U parametersestimated for the PBE H +2 molecule at varying bond lengths.Approximately parallel solid curves (blue, left) depict the U and U values required to recover the exact total-energy andKoopmans’ condition, assuming constant PBE subspace oc-cupancies. Open circles show the corresponding quantitiescalculated using a practical scheme based on the conventionalHubbard U and the PBE occupied eigenvalue (see text). Thetopmost curves (sand, right) show the fraction of the latter U (dashed) and U (solid) due to the Koopmans term, U K .On the same axis (dark red, right), we show that the averageDFT+ U subspace occupancy is insensitive to the U value. consistency effect, vanishing entirely for U = 2 N DFT U ,in subspace-uniform corrections such as those in ques-tion. For N sites equivalent one-orbital subspaces span-ning the energy window responsible for ε DFT , the densitynon-self-consistent U and U are derived from E exact ≈ E DFT + N sites (cid:0) U N DFT − U N (cid:1) /
2, where N DFT =Tr [ˆ n DFT ], and ε exact ≈ ε DFT + ( U − U N DFT ) /
2, inwhich the subspace overlap and spillage are also ne-glected.Since subspace response stiffening very typically re-sults from the application of a conventional DFT+ U cor-rection , it is promising to construct a charge non-self-consistent first-principles calculation scheme for U and U , e.g., for use in refining DFT+ U calculations in or-der to approximately enforce Koopmans’ condition. Forthis, let us suppose we have calculated a conventional U which reconciles the total energy reasonably, so that E exact ≈ E DFT + U N sites (cid:0) N DFT − N (cid:1) /
2. We maycombine this with the previous two equations and a fur-ther requirement for Koopmans’ condition at an accurateenergy, i.e., ε exact = E DFT [ N ] − E DFT [ N − U ≈ U (1 − N DFT ) (2 − N sites N DFT ) + U K , and N DFT U ≈ U (1 − N DFT ) (1 − N sites N DFT ) + U K , (9)where, with U K = 2 ( E DFT [ N − − E DFT [ N ] + ε DFT ),we define the ‘Koopmans U ’. We emphasize that onlyconvenient, approximate DFT quantities are used in these formulae. The interdependence U − U N DFT ≈ U (1 − N DFT ), for any value of N sites , reveals the roleof U in splitting U and U . In H +2 , Koopmans’ condi-tion pushes both up to considerably higher values thanare commonplace for the conventional U of DFT+ U ,which, following Ref. 32 and given N PBE , lies close toits regime of minimal efficacy for eigenvalue correction.The Koopmans fraction of each parameter, U K /U or U K / ( N DFT U ), generically denoted by ‘ U K /U ’ in Fig. 5,lies close to unity for U at all H +2 bond lengths, and itexceeds unity slightly when the U K and U -related con-tributions tend to cancel. U is also U K -dominated atshort bond lengths, at which U ≈ N DFT U , before ul-timately falling off to the average of U and U in thefully dissociated limit. If the outlined proposed schemeis applied to finesse an existing DFT+ U calculation thatis already accurate for recovering the total-energy, usinga linear-response or otherwise calculated U , call it U , then U = 0 eV is the appropriate value to use in ourapproximate formulae of Eq. 9. An approximately Koop-mans’ compliant DFT+ U calculation then results fromthe use of the parameters U + U K and U + U K /N DFT+ U ,in place of U and U . The proposed scheme may begeneralized to multi-orbital subspaces straightforwardly,in terms of the eigenvalues of ˆ n I DFT instead of N DFT .The constant- N DFT approximation may be replaced bya linear-response approximation, in terms of χ , or liftedentirely by means of a parametrization of the occupanciesand a numerical solution of the resulting equations.To conclude, we have proven analytically, with strin-gent numerical tests, that non-linear constraints are in-compatible with cDFT. It is not possible, therefore, toautomate systematic SIE corrections of DFT+ U type bymeans of cDFT, notwithstanding the great utility of thelatter, e.g., for correcting SIE by promoting broken sym-metry, integer-occupancy configurations well describedby approximate functionals . Nonetheless, we havefound that the cDFT free-energy functionals, dubbed‘generalized DFT+ U ’ functionals, offer the intriguing ca-pability of simultaneously correcting two central quanti-ties in DFT, the total-energy and the highest occupied or-bital energy. Our approximate formulae for the requiredparameters, which may differ greatly from the familiarHubbard U , offer a framework within which to furtherdevelop double-counting techniques and first-principlesschemes for the promising class of SIE correcting meth-ods based on DFT+ U , as well as opening uppossibilities for their diverse application. We envisagethat SIE correction schemes of two or more parametersmay also be useful for generalizing the exchange frac-tion of hybrid functionals , and for DFT+ U type cor-rections of perturbative many-body approximations suchas GW , the deviation from linearity of which is some-what analogous to that of approximate DFT . Forthe analysis and correction of spuriously self-interactingmulti-reference systems, we may learn much from the ex-act solution of minimal models .This work was enabled by the Royal Irish Academy –Royal Society International Exchange Cost Share Pro-gramme (IE131505). GT acknowledges support fromEPSRC UK (EP/I004483/1 and EP/K013610/1). GMand DDO’R wish to acknowledge support from the Sci-ence Foundation Ireland (SFI) funded centre AMBER (SFI/12/RC/2278). All calculations were performed onthe Lonsdale cluster maintained by the Trinity Centre forHigh Performance Computing. This cluster was fundedthrough grants from Science Foundation Ireland. ∗ [email protected] P. Hohenberg and W. Kohn, Phys. Rev. , B864 (1964). W. Kohn and L. J. Sham, Phys. Rev. , A1133 (1965). R. O. Jones, Rev. Mod. Phys. , 897 (2015). A. Jain, Y. Shin, and K. A. Persson, Nat. Rev. Mater. ,15004 (2016). J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz, Phys.Rev. Lett. , 1691 (1982). A. J. Cohen, P. Mori-S´anchez, and W. Yang, Science ,792 (2008). I. Dabo, A. Ferretti, N. Poilvert, Y. Li, N. Marzari, andM. Cococcioni, Phys. Rev. B , 115121 (2010). J. P. Perdew and A. Zunger, Phys. Rev. B , 5048 (1981). J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.Lett. , 3865 (1996). B. Kaduk, T. Kowalczyk, and T. Van Voorhis, Chem. Rev. , 321 (2012). H. J. Kulik, M. Cococcioni, D. A. Scherlis, and N. Marzari,Phys. Rev. Lett. , 103001 (2006). A. J. Cohen, P. Mori-S´anchez, and W. Yang, Chem. Rev. , 289 (2012). A. Svane and O. Gunnarsson, Phys. Rev. Lett. , 1148(1990). M. R. Pederson, A. Ruzsinszky, and J. P. Perdew, J.Chem. Phys. , 121103 (2014). T. Schmidt and S. K¨ummel, Phys. Rev. B , 165120(2016). D. Hofmann, S. Kl¨upfel, P. Kl¨upfel, and S. K¨ummel, Phys.Rev. A , 062514 (2012). D. Hofmann, T. K¨orzd¨orfer, and S. K¨ummel, Phys. Rev.Lett. , 146401 (2012). S. Kl¨upfel, P. Kl¨upfel, and H. J´onsson, J. Chem. Phys. , 124102 (2012). G. Borghi, A. Ferretti, N. L. Nguyen, I. Dabo, andN. Marzari, Phys. Rev. B , 075135 (2014). V. I. Anisimov and O. Gunnarsson, Phys. Rev. B , 7570(1991). V. I. Anisimov, J. Zaanen, and O. K. Andersen, Phys.Rev. B , 943 (1991). V. I. Anisimov, I. V. Solovyev, M. A. Korotin, M. T.Czy˙zyk, and G. A. Sawatzky, Phys. Rev. B , 16929(1993). S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J.Humphreys, and A. P. Sutton, Phys. Rev. B , 1505(1998). W. E. Pickett, S. C. Erwin, and E. C. Ethridge, Phys.Rev. B , 1201 (1998). M. Cococcioni and S. de Gironcoli, Phys. Rev. B ,035105 (2005). D. J. Cole, D. D. O’Regan, and M. C. Payne, J. Phys.Chem. Lett. , 1448 (2012). J. M. Garcia-Lastra, J. S. G. Myrdal, R. Christensen, K. S.Thygesen, and T. Vegge, J. Phys. Chem. C , 5568(2013). B. Himmetoglu, A. Floris, S. de Gironcoli, and M. Cococ-cioni, Int. J. Quantum Chem. , 14 (2014). M. Setvin, C. Franchini, X. Hao, M. Schmid, A. Janotti,M. Kaltak, C. G. Van de Walle, G. Kresse, and U. Diebold,Phys. Rev. Lett. , 086402 (2014). M. Shishkin and H. Sato, Phys. Rev. B , 085135 (2016). H. J. Kulik and N. Marzari, J. Chem. Phys. , 114103(2010). Q. Zhao, E. I. Ioannidis, and H. J. Kulik, J. Chem. Phys. , 054109 (2016). V. I. Anisimov, F. Aryasetiawan, and A. I. Lichtenstein,J. Phys. Condens. Matter , 767 (1997). K. Nakamura, R. Arita, Y. Yoshimoto, and S. Tsuneyuki,Phys. Rev. B , 235113 (2006). D. A. Scherlis, M. Cococcioni, P. Sit, and N. Marzari, J.Phys. Chem. B , 7384 (2007). D. D. O’Regan, N. D. M. Hine, M. C. Payne, and A. A.Mostofi, Phys. Rev. B , 081102 (2010). D. D. O’Regan, N. D. M. Hine, M. C. Payne, and A. A.Mostofi, Phys. Rev. B , 085107 (2012). S. Curtarolo, G. L. W. Hart, M. B. Nardelli, N. Mingo,S. Sanvito, and O. Levy, Nat. Mater. , 191 (2013). P. H. Dederichs, S. Bl¨ugel, R. Zeller, and H. Akai, Phys.Rev. Lett. , 2512 (1984). P.-L. Sit, M. Cococcioni, and N. Marzari, J. Electroanal.Chem. , 107 (2007). D. D. O’Regan and G. Teobaldi, Phys. Rev. B , 035159(2016). P. H.-L. Sit, M. Cococcioni, and N. Marzari, Phys. Rev.Lett. , 028303 (2006). The DFT+ U functionality available in the ONETEP linear-scaling DFT package was used with a hard (0 .
65 a cutoff) norm-conserving pseudopotential , 10 a Wannierfunction cutoff radii, and open boundary conditions . Q. Wu and T. Van Voorhis, Phys. Rev. A , 024502(2005). I. Dabo,
Towards first-principles electrochemistry , Ph.D.thesis, Massachusetts Institute of Technology (2008). Total-derivatives are used here to indicate that self-consistent density response effects are included. The Hub-bard parameters remain independent variables. D.-K. Seo, Phys. Rev. B , 033102 (2007). V. Atalla, I. Y. Zhang, O. T. Hofmann, X. Ren, P. Rinke,and M. Scheffler, Phys. Rev. B , 035140 (2016). H. Jiang, R. I. Gomez-Abal, P. Rinke, and M. Scheffler,Phys. Rev. B , 045108 (2010). W. Nelson, P. Bokes, P. Rinke, and R. W. Godby, Phys.Rev. A , 032505 (2007). M. Dauth, F. Caruso, S. K¨ummel, and P. Rinke, Phys.Rev. B , 121115 (2016). T. Olsen and K. S. Thygesen, J. Chem. Phys. , 164116(2014). C.-K. Skylaris, P. D. Haynes, A. A. Mostofi, and M. C.Payne, J. Chem. Phys. , 084119 (2005). A. M. Rappe, K. M. Rabe, E. Kaxiras, and J. D.Joannopoulos, Phys. Rev. B , 1227 (1990). G. J. Martyna and M. E. Tuckerman, J. Chem. Phys.110