Many people in the materials science and solid-state community are familiar with the acronym “DFT+U.” For those less familiar, this technique uses ideas from model Hamiltonians that permit the description of both metals and insulators to address problems of electron over-delocalization in practical implementations of density functional theory (DFT). Exchange-correlation functionals in DFT are often described as belonging to a hierarchical “Jacob’s ladder” of increasing accuracy in moving from local to non-local descriptions of exchange and correlation. DFT+U is not on this “ladder” but rather acts as an “elevator” because it systematically tunes relative energetics, typically on a localized subshell (e.g., *d* or *f* electrons), regardless of the underlying functional employed. However, this tuning is based on a metric of the local electron density of the subshells being addressed, thus necessitating physical or chemical or intuition about the system of interest. I will provide a brief overview of the history of how DFT+U came to be starting from the origin of the Hubbard and Anderson model Hamiltonians. This history lesson is necessary because it permits us to make the connections between the “Hubbard U” and fundamental outstanding challenges in electronic structure theory, and it helps to explain why this method is so widely applied to transition-metal oxides and organometallic complexes alike.

## I. INTRODUCTION

I will start by clarifying terminology and concepts for scientists like myself who are increasingly interested in research on transition-metal containing systems, be they molecular or solid-state (Fig. 1) in nature, at the intersection of the physics and chemistry communities. The physicist’s terminology is that a large number of these systems are characterized by *strong correlation*, where this phrase refers to the need *both* to describe atomic-scale Coulomb repulsion and exchange interactions rigorously as well as to describe the correlated probability densities of individual electrons. The local density approximation (LDA) description of the exchange-correlation (XC)^{1,2} in density functional theory (DFT) (see Becke’s recent perspective^{3}) is well-known to work best in systems with delocalized, slowly varying densities, such as metals. The paramount importance of understanding metals in the solid-state materials community led to early adoption of DFT within the LDA, but DFT was largely ignored by the theoretical chemistry community until the advent of the generalized gradient approximation (GGA) and hybrid functionals. Even to this day, the LDA still sees wide use in the physics community, and GGAs are unlikely to improve upon the LDA for description of metallic systems.

The chemist’s perspective differs significantly: here, we start primarily from a Hartree-Fock (HF) point of view, where HF provides the exact time-independent Schrodinger equation solution in the Born-Oppenheimer approximation for a system that can be described by a single Slater determinant. Coulomb repulsion and quantum-mechanical exchange, which arises directly from the use of the Slater determinant form, are exactly described. The exact electron repulsion integrals of Hartree-Fock are the missing ingredients of LDAs/GGAs that would make HF well-suited to describe much of the exchange and repulsion components of strong correlation that physicists refer to in these Mott insulator and other materials that are poorly described by a LDA or band theory^{4} description. In chemistry, the word correlation more commonly refers to two alternative and distinct phenomena beyond the Coulomb exchange and repulsion already incorporated into Hartree-Fock theory. The first is *dynamic* (short-range) correlation, wherein the probability density of one electron reduces the probability of finding another electron in the same space above and beyond the reduction in probability due to exchange. This feature is often added back into a predominantly Hartree-Fock wavefunction through expanding the total wavefunction in a basis of other Slater determinants that preferentially occupy increasingly diffuse molecular orbitals. The second is *static* correlation, wherein no single Slater determinant adequately describes the wavefunction, and a superposition of nearly equally dominant determinants is needed, which may be either purely long-range, as in H_{2}^{+} dissociation,^{5} or may also present at short-range, as in cases with degenerate frontier orbitals (e.g., the singlet diradical trimethylenethane^{6}).

From both the chemistry and the physics perspectives, there is a clear need for methods that correctly describe electron-repulsion of well-localized electrons, and Hartree-Fock succeeds quite well at this. However, it is also necessary to correctly describe subtle dynamic correlation effects that govern bonding and permit some measure of delocalization, either in a bulk metal or in a molecule, but are absent from Hartree-Fock. The Hubbard model first arose out of the recognition that rather than describing all electrons as a sea in a weak field of point nuclei, appropriately modeled via the homogeneous electron gas through the LDA in DFT; many electrons retain their well-localized, atomic or molecular character in certain classes of solid-state materials. For this reason, I will make the argument here that the language that both a physicist and a chemist should be able to understand is that the degree of localization or delocalization is a problem in both approximate density functional theory and in Hartree-Fock. Hartree-Fock or even a hybrid functional may struggle to describe a metal, i.e., these approaches will not properly delocalize the electrons, while modern implementations of density functional theory with common XC approximations will often over-delocalize electrons in cases where the electrons should be well-localized. Here, we restrict ourselves to cases where the problem is again local correlation and exchange, not long-range delocalization problems as in H_{2}^{+}. This localization/delocalization problem is pernicious: two primary examples are that over-delocalization turns solid state insulators into metals in condensed-matter physics and the preferentially stabilizes of low-spin transition-metal complexes in inorganic chemistry. Conversely, over-localization in Hartree-Fock can overestimate barrier heights and is generally not useful for understanding trends in chemical bonding, where dynamical correlation is crucial.

Chemists and physicists are in full agreement that an electronic structure method should have a good balance in describing these competing effects, but the approaches to solve this problem may not be universally agreed upon. One of the most popular solutions to this problem, which originates from the chemistry community, has been to find a marriage between the two using a weighted average of DFT and HF exchange. Examples of GGA hybrid XC functionals are the nearly exclusively employed B3LYP in chemistry^{7–9} or the equivalent PBE0 hybrid^{10,11} in physics. More recently, meta-GGAs^{12} that incorporate higher order dependence on the density (e.g., the popular M06 family of “Minnesota functionals”^{13,14}), or range-corrected functionals that use an error function to switch between a short-range and a long-range definition of the functional,^{15} have also been developed both as pure and hybrid XCs. Chemists may be curious why physicists have lagged behind their advances in employing these hybrid functionals. In fact, the PBE^{16} generalized-gradient approximation remains the most widely used functional in the solid-state community,^{17} despite its well-known shortcomings for describing many key systems of interest. While plane wave periodic boundary condition codes exhibit challenges in describing core electrons^{18} and require careful use of pseudopotentials,^{19–21} these approaches remain a necessity for describing extended bulk systems and have performance advantages^{22} for describing well-delocalized electrons versus using a localized-basis set approach with large basis sets.

Within the plane wave framework, calculating *exact exchange* from Hartree-Fock on the Kohn-Sham single-particle orbitals is very costly^{23} and intractable for most systems of interest. Considerably more recently, physicists do now also employ hybrid exchange on smaller systems or put in the Herculean effort to do the same for larger systems, most commonly by using range-separated hybrids^{24,25} that only incorporate HF exchange in the short range. These short-range hybrids work well in the solid because effects sensitive to descriptions of dynamic correlation and exchange dominate over the accuracy recovered through correct description of the computationally costly long-range, static correlation (see a more thorough discussion in Ref. 26). Hybrids largely remain an enigma because the parameter describing the “correct” amount of Hartree-Fock exchange is highly system dependent,^{27} and range-separation adds additional complexity.^{28} Only recent work^{29,30} has started to point to ways these parameters can be tuned in response to system properties. In this perspective, I will describe DFT+U, an alternative method that originates from successful model Hamiltonians in the solid-state physics community but has broad applicability to systems of relevance in both physics and chemistry.

Like hybrid functionals, DFT+U has its own advantages and disadvantages, but they are less well understood. I will provide a brief overview of the history of DFT+U: first, its origins from the Hubbard and Anderson models to approximate electron repulsion integrals locally for well-localized electrons in a solid; next, its migration to a simplified model Hamiltonian with parameters extracted from the LDA; and finally, its evolution to a simple, practical tool. I will also provide an outlook on the strength, weaknesses, and practical challenges as well as some discussion of relative benefit of key extensions to DFT+U theory.

## II. THE HUBBARD MODEL IN SOLID-STATE PHYSICS

Those with a strong solid-state physics background are likely familiar with both the Hubbard model^{31} and its predecessor, the Anderson impurity model.^{32} These models carry gravitas in the physics community due to their successes in describing complex materials properties, such as correctly describing insulators normally predicted by band theory to be metals (e.g., Mott insulators^{33}), describing the Kondo effect,^{34} or even demonstrating some ability to characterize high-Tc superconductors. In all of these cases, the model is heuristically describing how a magnetic moment will localize on an atom rather than delocalizing across the solid. For everyone feeling a little lost, I will reiterate the physical and chemical motivation of these models and provide a brief explanation of the model Hamiltonian terminology that subsequently became associated with the DFT+U electronic structure method.

^{35}Hubbard model is able to explain the transition between the limits of the band model in a conductor and the localized-limit in an insulator. These two limits are achieved by tuning competing interactions: the kinetic energy operator and weak electron-nuclear attraction favor delocalization, but there is also a Coulomb repulsion term associated with pairing electrons on an atomic “site.” In commonly employed second quantization notation, the full Hubbard Hamiltonian may be expressed

^{36}terminology where the electrons are described as being able to move only between orbitals μ and ν on nearest neighbor sites labeled

*i*and

*j*(Fig. 2). The explicit form of the hopping integral is

*t*. The second part of the expression in the Hubbard Hamiltonian corresponds to electron-electron repulsion integrals of orbitals located on up to four distinct atomic sites

*i*,

*j*,

*k*,

*l*with orbitals μ, ν, σ, τ.

^{31}and was therefore the only interaction that required treatment. The average value of the Coulomb energy required to pair two electrons on a single atom is usually also treated as a single parameter,

*U*, although Hubbard referred

^{31}to this parameter only as

*I*. While in the original work, a single pair of electrons was considered, the model is traditionally employed on an entire band of localized

*d*or

*f*electrons. The Hubbard model with modern notation is greatly simplified

*i*is denoted

*N*

^{I↑}or

*N*

^{I↓}, respectively.

In 1961, the related Anderson impurity model was first introduced,^{32} two years prior to the 1963 work of Hubbard on localized electrons in solids.^{31} Anderson also first proposed describing the on-site Coulomb integrals with the notation “*U*,” but Anderson’s model is not just about *U* terms between electrons on the same site. The Anderson impurity model considers interactions absent from the Hubbard model, which are relevant to recent extensions^{37} to DFT+U that have begun to consider again the importance of balancing on-site, or intra-subshell, interactions with inter-site, or inter-subshell, hybridization. Specifically, the Anderson impurity model incorporated the same elements as those in Hubbard’s model Hamiltonian but also considered the hybridization of localized *d* electrons with more delocalized *s* electrons. Therefore, the Anderson model provides a more nuanced understanding of chemical bonding, which we found to be necessary for describing transition-metal hydrides (a discussion of employing *U* on both 3*d* and 4*s* electrons simultaneously is covered in our Ref. 38 and also in Ref. 39). Additionally, Anderson first posed that it was more important to describe repulsion of opposite-spin electrons (Coulomb integrals or *J*), rather than attraction of parallel ones (exchange integrals or *K*) to describe localized magnetic moments correctly, a key approximation often used in modern DFT+U practice, denoted as *U*_{eff} = *U* − *J*.

While there have been occasional suggestions that for specific materials, an explicit exchange term, *J* in Hubbard model notation, is critical for describing key interactions, I have only limited experience with the method. In the few cases, I have had a chance to test, the two effects counterbalance each other. That is, at its heart, as first proposed by Hubbard and Anderson, the *U* discourages pairing of electrons, while the *J* encourages pairing of electrons. I have anecdotally observed preference of low-spin states in GGAs for molecules. In a solid, over-pairing of electrons corresponds to an antiferromagnetic insulator or to a metal. A strong *U* term will therefore prefer ferromagnetic ordering. If using tuned and not calculated values, the only discrepancy that arises between using an effective *U _{eff}* vs.

*U*alongside a

*J*term would be for cases where at least three distinct spin states are considered, which is quite uncommon in most applications of DFT+U.

## III. MODEL PARAMETERS FROM ELECTRONIC STRUCTURE THEORY

Both Hubbard and Anderson proposed to approximate terms in their Hamiltonians from a single-particle, Hartree-Fock approximation. Careful readers may have already noted that Kohn-Sham^{40,41} DFT was introduced after these model Hamiltonians. As both DFT and the Hubbard model matured in parallel in the 1970s and into the 1980s, a merger between the two theories was imminent. The local density approximation had become increasingly popular in physics for describing the band structure of simple metals. Concurrently, the Hubbard and Anderson models were being increasingly applied to insulating^{33} and superconducting^{42} materials that captured much interest in the solid-state physics community. In the early 1970s,^{43} a number of scientists began^{43–50} in earnest to compute model Hamiltonian parameters directly from electronic structure calculations. Much of what is known today as DFT+U was first developed in this context, but few people are aware of this work. So, I will provide a review of the major developments.

The earliest parameter calculations involved the use of Hartree-Fock on atoms in modified environments to estimate values of *U* on 3*d*^{43} and 4*f*^{44,45} electrons that were found to be much lower than Hubbard’s early estimates of *U* ≈ 20 eV,^{31} a technique that has recently been revisited in modern DFT+U.^{51,52} The lower *U* value takes into account how interactions of atomic electrons will be reduced due to *screening* or *renormalization* from the response of all other delocalized electrons. Specifically, an increase in the number of electrons on the Hubbard atom site will shift the electrostatic potential upward but that shift is partially compensated by the other electrons moving away from the atom center. The importance of screening is also sometimes emphasized in modern linear-response *U* calculations.^{53–55}

^{46}within the linear muffin tin orbital (LMTO)

^{56}approximation were used early on for the purpose of estimating

*U*in the Anderson impurity model. There, constrained DFT was used to obtain energies at integer values of charges

^{46}to obtain

*U*as

*U*is the difference between the ionization potential (

*IP*) and electron affinity (

*EA*) for a particular atom (

*I*) and subshell (

*nl*) with respect to the rest of the molecule with respect to a neutral atom with $ N nl I $ electrons, including other subshells on the atom

*I*. That is,

*U*here is the energy required to pair two electrons on a site by removing one electron from an equivalent site. This definition has several other representations that may look more familiar to the reader. For instance, one may recognize that Eq. (4) is a finite difference representation of the second derivative of the total energy with respect to the number of electrons in the

*nl*subshell,

*U*in terms of orbital energies are also possible via Janak’s theorem

^{57}in finite difference form

*I*with subshell

*nl*and magnetic quantum number

*m*. Eigenvalues are evaluated for $ N nl I + 1 2 $ electrons to obtain the electron affinity and $ N nl I \u2212 1 2 $ to obtain the local ionization potential.

_{l}^{47,58}It naturally follows that

*U*may be expressed as the first derivative of the representative eigenvalues,

*I*and subshell

*nl*, it is expected that a

*U*term in a model Hamiltonian or in DFT+U should strongly modulate the overall band gap or HOMO-LUMO splitting of the system. Total energy-based or eigenvalue-based estimates of

*U*may thus be obtained either through electron addition and removal or through constrained DFT

^{46}calculations.

*U*may be equivalently cast in terms of the chemical potential,

*μ*,

*δ*→ 0, where the first term is the ionization potential of the

*nl*subshell and the second term is the negative of the electron affinity of the

*nl*subshell.

^{59}This expression for

*U*clarifies that we are adjusting the local chemical potential or the preference for electrons to be occupied on a certain

*nl*subshell with respect to the rest of the system. Here, we are saying that the chemical potential of a certain set of localized manifold electrons is too high or too low with respect to the surrounding system. Electrons will not preferentially occupy this site if they can hybridize with and delocalize to other orbitals (bands). Physicists will recognize this problem as the case where the d-band width is too large and placement is too high with respect to the Fermi energy. In chemistry, this problem arises with incorrect spin state ordering and incorrect spatial symmetries in open shell 3

*d*systems. We first identified and introduced the value of a +

*U*to correct level placement and hybridization in molecular, biological, and heterogeneous catalytic materials

^{38,54,60–64}(see also schematic in Fig. 3). The chemical potential-based tuning picture of the

*U*has been also validated in the successful systematic tuning of computed battery material voltages

^{65}and in oxygen vacancy energies in ceria.

^{66}

When calculating the *U* within any approach, the upward shift of the potential before incorporating screening is typically referred^{53} to as the *bare* response and may in practice be obtained from the first iteration of a self-consistent field calculation. Many of the early developments calculating model parameters took place during the era of linear muffin-tin orbital DFT calculations,^{56} which were a predecessor to the now-popular projector-augmented wave method.^{67} Early constrained DFT calculations were also extended to consider inter-site Coulomb interactions, i.e., the hybridization between 2*p* electrons and 3*d* electrons in a transition metal oxide,^{48} but difficulties were encountered due to strong dependence on the size of the super-cell employed, an observation later echoed in “simplified” DFT+U.^{53} Such “inter-site” interactions were considered from the early days of both Anderson and extended Hubbard models.

Careful work^{47} on parameterizing Anderson impurity model interaction terms also compared the relative utility of using total energy differences versus eigenvalue definitions for the calculation of intra-site and inter-site interaction parameters as well as the use of the constraint to vary occupations (fixed valence occupation versus free) while computing the eigenvalues. Such work is likely worth revisiting when weighing the options for calculating *U*, as the typically employed approach necessitates an inversion of a response matrix^{53} (that we noted, Ref. 38 introduces numerical noise) while energy differences would not suffer from this problem.

Within all this careful work, one might begin to ask oneself whether LDA should be sufficient to describe correctly parameters for systems that have localized electrons. In fact, it was observed that inter-site parameters (i.e., the extent of hybridization between localized and delocalized states) were overestimates due to a poor ground state induced by severe self-interaction error.^{68} Closely following on this work, others observed that poor LDA ground state properties prevented calculation of *U* and that self-interaction corrections were needed to improve the LDA ground state description of the material.^{69} In my view, Ref. 69, which has only been cited 17 times (according to Google Scholar), is the true catalyst for what shortly thereafter was introduced as DFT+U. The recognition that parameters for model Hamiltonians from electronic structure would depend strongly on the density of the system was later echoed in our own self-consistent *U* calculation scheme.^{54} It is self-evident that the Coulomb energy should be calculated on the ground state density of the relevant low-lying DFT+U state. While in our original work, we considered this self-consistent calculation (see Ref. 54), my experience is that a self-consistent *U* is only necessary when strong qualitative changes are observed in the density with applied *U*, i.e., a GGA metal becomes a GGA+U insulator (see also Sec. IV).

## IV. THE CODIFICATION AND IMPLEMENTATION OF PLANE-WAVE LDA+U

Careful readers may be surprised to see how much of the groundwork of modern DFT+U^{70–73} in the 1990s was already established by work in prior decades. I know I certainly was surprised when I started to pull this perspective together. The recent review by Himmetoglu *et al.*^{74} does an excellent job of covering many of the finer details of the LDA+U implementation common in a number of plane wave codes. It was these efforts in implementation during the mid-2000s that broadened the utility of DFT+U from its origins in the LMTO formalism.

*E*

_{DFT}), (2) the Hubbard term (

*E*

_{Hub}) that explicitly models the Coulomb energy associated with the localized atomic orbitals of interest (e.g., 3

*d*electrons), and (3) removal of a double counting term (

*E*

_{DC}) that corresponds to a mean field description of these interactions obtained in the homogeneous electron gas limit with a standard DFT functional,

*nl*on site

*I*with associated magnetic quantum numbers,

*m*, upon which a Hubbard term is applied and a double-counting term that depends on total occupation of that subshell,

_{l}*N*

^{Iσ}. Most commonly, the double counting term is obtained as

^{53,73}a simplifying assumption, as I have already pointed out, is that

*U*

_{eff}=

*U*−

*J*yielding an overall expression for the DFT+U energy as

*nl*subshell on atom

*I*, σ is a spin index, and $ U nl I $ is the effective electron-electron repulsion interaction parameter that may be calculated

^{38,51,53,54,64}or, more commonly, tuned

^{75–77}and is specific to each atom and subshell. A variety of definitions are available for the occupation matrices, $ n nl I , \sigma $, that enter into the +U energy functional. In the solid-state physics community, where DFT+U is most commonly employed, occupation matrices are obtained by projecting plane-wave-based molecular orbitals or bands onto a localized basis set. When starting my studies of DFT+U for molecules, I continued to use plane-wave codes because those were the tools that had been already implemented. Using those tools, we extensively considered the role that +U term might play on molecules in order to observe trends in chemical bonding that a +U term has had. A thorough discussion of possible basis set definitions for plane-wave-based DFT+U is presented in Ref. 74.

*U*for Hubbard model parameters of the system of interest (see Sec. III). While a number of expressions were introduced in Sec. III, a variant of constrained DFT is the most commonly preferred approach

^{53,55}to compute

*U*in DFT codes. Here, a weakly perturbing potential, $ \alpha nl I $, is applied to the

*nl*subshell of atom

*I*via a term,

*E*

^{lin},

*nl*subshell of atom

*I*. We may then compute

*U*from the second derivative of the total energy minus the

*E*

^{lin}term with respect to total occupations on site

*I*. Note that this perturbation induces a natural shift in occupations often referred to as the

*bare*,

*non-interacting*response that must then be subtracted to obtain the unphysical curvature associated with the interacting electron density. The

*U*obtained in this way is then expressed as

*nl*′ subshell of atom

*J*with respect to perturbations of

*nl*levels on atom

*I*, i.e., $ d ( N n l \u2032 J ) d \alpha nl I $, is computed and inverted to obtain terms for calculating

*U*.

^{54}to compute

*U*directly on the DFT+U density, by plotting the computed

*U*

_{out}against the applied

*U*

_{in}permits extrapolation back to a so-called self-consistent U (

*U*

_{scf}) at

*U*= 0,

*U*is most useful for cases where the DFT+U electronic state is inaccessible with the starting underlying exchange-correlation functional. In my own experience, while such DFT+U electronic states may be higher in energy, it is still possible to metastably converge these states, in which case the self-consistent calculation of

*U*is not necessary. For instance, across 28 transition metal complexes involving carbide, nitride, oxide, and fluoride bonds,

^{38}self-consistent

*U*values deviated on average by only 0.23 eV from the standard linear response value, corresponding to an average of a 9% difference between the two numbers. In the absence of qualitative changes in the density, the linear response

*U*and self-consistent

*U should*be in strong agreement, otherwise we have a substantial failure of our theory. Often in my own work, I have focused instead on metastably converging multiple low-lying states through biasing initial guesses with differing magnetizations or using the initial guess from a self-consistent result obtained at a different

*U*value. Much of these efforts are what made possible our careful convergence of many low-lying iron dimer states

^{54}or of the convergence of many differing intermediates in the conversion of methane to methanol.

^{60}Some of the methods I developed to converge all of these states have been made more systematic and coded up as an approach, some may know as “U

_{ramping}.”

^{78}I am not aware, however, of anyone having established a systematic scan for low-lying electronic states (states with differing magnetization), which is a necessary tool I hope to make available soon.

It is also worthwhile to note that while *U* may be straightforwardly computed with little more than the cost of a handful of energy calculations, some observe benefit in obtaining improved accuracy using DFT+U with computed values,^{79} and many others prefer to use *U* as a parameter to tune agreement with experiment on a key observable. In my own work, I have found the calculated *U* values themselves to be illuminating descriptors for the underlying electronic structure, e.g., in comparing covalent and ionic bonding in isoelectronic systems,^{38} noting changes in local environment in a reaction coordinate,^{60} or identifying when hybridization of 3*d* and 4*s* states was critical.^{38} Nevertheless, I have also encountered challenges with computing *U* from linear response, including sensitivity to numerical noise^{38} and values that ran counter to simple hybridization arguments, motivating the call for renormalization of computed *U* values.^{80}

## V. BAND GAPS AND DFT+U

*U*offers and has broad applicability to all systems of interest. Nevertheless, let us look again at the band gap problem within DFT. The single particle Kohn-Sham orbital energies are well-known to not necessarily correspond to the fundamental gap (FG),

^{81,82}which is instead defined as

_{KS}) arises from the energy gap between single-particle Kohn-Sham orbitals and the second term (Δ

_{XC}) comes from the discontinuity of the exchange-correlation potential.

It is more reasonable to say that in the limit of DFT+U only being dominated by the *U* effect, then the band gap would be exactly equal to the value of *U* for the underlying atomic system. This value would be a huge overestimate of the true band gap of the solid-state system. However, when mixed with LDA’s underestimation of the band gap, we have a cancellation of errors that results in band gaps generally being pushed in the correct direction. Nevertheless, I would like to borrow the acronym Axel Becke proposed in his perspective “OOO” for occupied orbitals only and extend it to suggest that DFT+U is best applied only to partially occupied subshells and for ground state properties (i.e., occupied orbitals only). While +U can and will continue to be used to open gaps, one should not assume it would ever be possible to obtain quantitative gaps unless one tunes the *U* to obtain that gap. A generalization of the Δ-SCF method for solids, known as Δ-sol, is likely a more fruitful approach.^{83} Regardless, to exclusively tune gaps is missing the point—level placement and hybridization, i.e., which states mix with which—is what practical DFT is getting most wrong. If we can improve the unoccupied single-particle orbital energies, that is fortuitous but not necessarily on firm foundation. This lack of firm foundation has perhaps emboldened detractors of DFT+U. The chemical potential or local ionization potential/electron affinity picture is more attractive for understanding what adjusting *U* will do to a system.

## VI. DFT+U FOR THEORETICAL CHEMISTRY

^{60}For localized basis sets, it is most convenient to utilize elements of the Mulliken population matrix (

**q**), which are defined as

**P**) and overlap (

**S**) matrices. The occupation matrices and applied

*U*values used in +U corrections are specific to a particular subshell that is defined by

*n*and

*l*quantum numbers and localized on an atom,

*I*. If the system is closed-shell, population matrix values are reduced by a factor of two (i.e., a fully occupied orbital always corresponds to a matrix element of 1). This method may be straightforwardly

^{84}combined with either hybrids or even Hartree-Fock theory, so here I will use the notation HF/DFT+U. I emphasize the possibility of using the +U method within Hartree-Fock, even though it is formally self-interaction free and tends to suffer from electron overlocalization, as was discussed earlier in the context of metals. The most transparent utility is then to apply negative values of

*U*, which in turn discourage integer occupations and encourage delocalization. The secondary benefit stems from the fact that the +U correction is ultimately a penalty on the extent of fractionality in the occupation matrix and therefore may be used to tune relative energies, even in Hartree-Fock, as long as the two systems being compared have different relative values for the penalty factor Tr[

**n**(

**1**−

**n**)] (or replace

**n**with

**q**if using Mulliken notation). Such an approach works best for “rigid” molecular orbitals where the primary effect would be in tuning relative energetics of various electronic states rather than dramatically shifting the hybridization of the underlying molecular orbitals. In my group, we have found some use for this occupation-matrix-oriented way of thinking for treating both open shell radicals and closed shell systems on more even footing in both DFT and HF as well as for correcting for basis set incompleteness or other energetic imbalances in electronic structure methods.

*I*, and associated principle quantum numbers,

*n*and

*l*, of the atomic basis functions, and σ corresponds to the spin index. The total energy correction is additive to a result from DFT, and the +U correction is incorporated into the self-consistent calculation through direct modification of the potential,

*v*

_{μν}) is explicitly

*U*term systematically tunes the single-particle orbital energies (see example in Fig. 4).

Density functional theorists may be hesitant to interpret orbital energies, but they are both a necessity to understand and interpret hybridization and carry significant qualitative and quantitative information (see detailed discussion in Ref. 85). The projected density of states is often used in materials science and catalysis to identify the relative energy windows in which contributions coming from differing orbitals and permits interpretation of properties. When we instead deny the validity and importance of orbital energies or shapes because they are not the exact “interacting” system, we deny ourselves the opportunity to glean more information from our simulations. Regardless of personal preference, what the *U* primarily does is to shift the orbital energies or eigenvalues of the dominantly atomic-like molecular orbitals. The derivative of the potential here is in the atomic basis but would necessitate transformation back to a molecular basis, and the less-atomic-like molecular orbitals will have little shift to their potential.

## VII. OUTLOOK ON THE BROADER ROLE OF MODEL HAMILTONIANS WITHIN DENSITY FUNCTIONAL THEORY

I have been a developer and user of DFT+U theory for around 10 yr now and would like to provide newcomers to this area a brief summary of some of my experience. What I learned from working with and playing with these tools was that while a *U* was a good first order approximation, adding on a Hamiltonian that preferentially localizes electrons can always be used to over-localize these electrons and eliminate crucial, chemical bonding. Typically, calculated values of *U* do not exclude bonding, but sufficiently large values of *U* will elongate bonds and alter angular potentials leading to qualitative changes in geometry.^{80} Instead, alternative methods include constraining bonds,^{80} varying *U* along the bonding coordinate^{64} or considering intersite terms^{80} (i.e., *V*) is necessary.

In my own work, I am most interested in catalysis at single-site transition metal complexes because they provide the optimum combination of simplicity in electronic structure and catalytic control and selectivity. The challenges associated with simply predicting the ground state of unusual systems brought me to introduce a few methods that may have broader applicability to solid-state materials as well as the original organometallics for which I developed these approaches. DFT + *U*(**R**)^{64} permits the direct comparison of DFT+U energies obtained at differing values of U. While in earlier work, I had advocated for computing and applying average values of *U*^{38,54,60–63} and had even estimated the uncertainty in this approach as about 1 kcal/mol introduced for 1 eV off “optimal U” used; cases with large variation in computed values of *U* were still problematic. The original DFT + *U*(**R**)^{64} paper contains several key pieces: (1) an approach to predict a gradient of *U* using information obtained only from a geometric gradient and the typical linear response calculation, (2) identification that relative energies, once referenced, or, alternatively, forces or stresses can be compared at differing values of *U*, (3) suggestions for what that reference could be in a united atom or dissociated atom limit, and (4) interpolation schemes of both one dimensional reaction coordinates and two-dimensional potential energy surfaces. Such approaches may be straightforwardly applied to solid-state materials.

There have also been many successful applications of the standard DFT+U theory to problems in catalysis^{86} and materials science,^{87–89} but many more applications of DFT+U are possible with small extension to the method. The +U energy correction also produces some useful physical behaviors that are transferable to other problems in electronic structure theory unrelated to self-interaction error. Namely, the energy correction depends exclusively on an occupation matrix that is very sensitive to the electron configuration for a particular atom and subshell, which is in turn influenced by the bonding environment. If an electronic structure method consistently overstabilizes certain geometric or electronic configurations, the *U* parameters can be used to correct this imbalance. The +U potential is linear in the occupations of the orbitals to which it is applied. With positive *U* parameters, this correction corresponds to lowering the energies of molecular orbitals comprised of over 50% of a localized atomic orbital and destabilizing molecular orbitals with less than 50% of the localized atomic orbital. When the subshell and atom correspond to frontier molecular orbitals, as is often the case, the molecule’s ionization potential obtained from Koopmans’^{90} or Janak’s^{57} theorem is also adjusted. The +U correction may therefore be thought of as tuning the effective electronegativity of both substituent atoms, as defined by relative energies of substituent molecular orbitals and the overall molecule, making it possible to envision applying this to correct basis set incompleteness or relative imbalance between radical and non-radical species in DFT.

## VIII. CONCLUDING REMARKS

For newcomers to this area, it may be a surprise to find that DFT+U calculations are hard to converge or result in multiple possible solutions—please understand that is because DFT+U is only needed in cases where the electronic structure is complex and the system has multiple possible ground or low-lying solutions. These results are not “exotica” of DFT+U but they are in fact signs that you have chosen an interesting problem to study. In the beginning, I experienced such problems when I wrongly presumed DFT+U or two atoms of iron were both simple challenges that would be a short detour on my way to studying larger problems in materials chemistry. I enjoyed having the chance to have a Ph.D. worth of time and beyond to really start to learn why, and I have since enjoyed learning new things I still “do not know” about DFT+U in that time since.

I hope this perspective has compelled you to take a second look at a well-known acronym and the scientists who made it what it is today. I know that is what I got most for myself from writing it for you. Please let me know if anything I have said is unclear. I am always eager to receive feedback on how I can communicate science better.^{91} I cannot tell you yet where I think DFT+U is going, but I think there is still much to be learned from the not often cited papers that led to the birth of what is now a very widely used tool.

## Acknowledgments

I am grateful for my students: Mr. Tim Ioannidis, Ms. Helena Qi, Ms. Lisi Xie, and Ms. Qing Zhao. Many thanks to Richard Braatz, Catherine Drennan, Craig Fennie, Francois Gygi, Giulia Galli, Todd Martinez, and Adam Steeves for helpful discussions during this time. I am also grateful for the people who inspired the Ph.D. thesis topic that led me to bring DFT+U to quantum chemistry and transition-metal-containing molecules from the solid-state physics community (Matteo Cococcioni, Filippo DeAngelis, Stefano DeGironcoli, and Nicola Marzari, among others). I am grateful to them for giving me such a fun toy to play with, but it was a challenging puzzle; I thought it was my job to crack as soon as I possibly could. I learned DFT+U and the iron dimer before I learned larger things, such as the thousands of atoms in the proteins I studied in my postdoctorate with Todd Martinez at Stanford. I am indebted to my Ph.D. advisor at MIT, Nicola Marzari (now at EPFL) for allowing me the chance to think like a “quantum chemist” in a computational materials science group. Last but not least, thanks to our cluster for bearing with us as we stress test him (Mr. Gil Braltar) with some difficult calculations. H.J.K. holds a Career Award at the Scientific Interface from the Burroughs Wellcome Fund. Support was also provided by the Research Support Corporation (MIT), NSF ECCS-1449291, the MIT Energy Initiative, as well as startup funds provided by the department of chemical engineering at MIT.

## REFERENCES

*The Self-Consistent Field for Molecules and Solids*

*Electronic Structure: Basic Theory and Practical Methods*

**9**, 523 (2013).

*Fuel Cell Science: Theory, Fundamentals, and Bio-Catalysis*

Some of you may be most familiar with me through my electronic structure tutorials that first started in an attempt to better explain DFT+U and help people avoid some of my early challenges in applying the method broadly (http://hjklol.mit.edu/Tutorials). It is a delight for me to hear occasionally that these resources are still helpful to some of you, but unfortunately my current job keeps me busier than I expected and away from adding new tutorials or responding to queries. I hope that my students will someday pick up the torch and begin our tutorial series anew as they and I learn new codes and skills together.