Density functional theory in transition-metal chemistry: a self-consistent Hubbard U approach.
ABSTRACT Transition-metal centers are the active sites for a broad variety of biological and inorganic chemical reactions. Notwithstanding this central importance, density-functional theory calculations based on generalized-gradient approximations often fail to describe energetics, multiplet structures, reaction barriers, and geometries around the active sites. We suggest here an alternative approach, derived from the Hubbard U correction to solid-state problems, that provides an excellent agreement with correlated-electron quantum chemistry calculations in test cases that range from the ground state of Fe2 and Fe2- to the addition elimination of molecular hydrogen on FeO+. The Hubbard U is determined with a novel self-consistent procedure based on a linear-response approach.
- Physical Review B 09/2014; 90(12). · 3.66 Impact Factor
- [Show abstract] [Hide abstract]
ABSTRACT: The halogenases are a group of enzymes that have only come to the fore over the last 10 years thanks to the discovery and characterization of several novel representatives. They have revealed the fascinating variety of distinct chemical mechanisms that nature utilizes to activate halogens and introduce them into organic substrates. Computational studies using a range of approaches have already elucidated many details of the mechanisms of these enzymes, often in synergistic combination with experiment. This Review summarizes the main insights gained from these studies. It also seeks to identify open questions that are amenable to computational investigations. The studies discussed herein serve to illustrate some of the limitations of the current computational approaches and the challenges encountered in computational mechanistic enzymology.Frontiers in Chemistry 11/2014; 2:98.
- Philosophical Magazine 08/2014; 94(26):3052-3071. · 1.43 Impact Factor
arXiv:cond-mat/0608285v1 [cond-mat.soft] 12 Aug 2006
Density functional theory in transition-metal chemistry: a self-consistent Hubbard U
Heather J. Kulik, Matteo Cococcioni, Damian A. Scherlis,∗and Nicola Marzari
Department of Materials Science and Engineering,
Massachusetts Institute of Technology, Cambridge, MA 02139, USA
(Dated: February 5, 2008)
Transition-metal centers are the active sites for many biological and inorganic chemical reac-
tions. Notwithstanding this central importance, density-functional theory calculations based on
generalized-gradient approximations often fail to describe energetics, multiplet structures, reaction
barriers, and geometries around the active sites. We suggest here an alternative approach, derived
from the Hubbard U correction to solid-state problems, that provides an excellent agreement with
correlated-electron quantum chemistry calculations in test cases that range from the ground state
of Fe2 and Fe−
determined with a novel self-consistent procedure based on a linear-response approach.
2to the addition-elimination of molecular hydrogen on FeO+. The Hubbard U is
Transition metals are central to our understanding of
many fundamental reactions, as active sites in naturally-
existing or synthetic molecules that range from metallo-
porphyrins and oxidoreductases  to alkene metathesis
catalysts  to light-harvesting photosynthetic complexes
.Despite this relevance, most electronic-structure
approaches fail to describe consistently or accurately
transition-metal centers. Examples include neutral and
charged iron dimers , FeO+, Mn(salen) epoxidation
catalysts , or hemeproteins .
In this Letter, we argue that generalized gradient
approximations (GGA) augmented by a Hubbard U
term, already very successful in the solid state [10,
11], also greatly improve single-site or few-site energies,
thanks to a more accurate description of self- and intra-
atomic interactions. Nevertheless, U is not a fitting pa-
rameter, but an intrinsic response property: as shown
by Cococcioni and de Gironcoli , U measures the
spurious curvature of the GGA energy functional as a
function of occupations, and GGA+U largely recovers
the piecewise-linear behavior of the exact ground-state
energy. U is determined by the difference between the
screened and bare second derivative of the energy with
respect to on-site occupations λI
orbital, and I the atomic site). While in the origi-
nal derivation U was calculated from the GGA ground
state, we argue here that U should be consistently ob-
tained from the GGA+U ground state itself. This be-
comes especially relevant when GGA and GGA+U differ
qualitatively (metal vs. insulator in the solid state, dif-
ferent symmetry in a molecule). To clarify our approach,
we first identify in the GGA+U functional the electronic
terms that have quadratic dependence on the occupa-
i(i is the spin-
i(1 − λI
The first term represents the contribution already con-
tained in the standard GGA functional, modeled here
as a double-counting term, while the second term is the
customary “+U” correction. Therefore, Uscf represents
the effective on-site electron-electron interaction already
present in the GGA energy functional for the GGA+U
ground state when U is chosen to be Uin. Consistency
is enforced by choosing Uinto be equal to Uscf. The U
obtained from linear-response (labeled here Uout) is
also obtained by differentiating Eq. 1 with respect to λI
where m = 1/?
degeneracy of the orbitals whose population is changing
during the perturbation (to linear order, δλI
principle Uscf depends on Uin, we find it to be constant
over a broad interval, as apparent from Fig. 1: Uout
is linear in Uin for the relevant range of Uin ∼ Uscf.
Thus, from few linear-response calculations for different
Uin ground states we are able to extract the Uscf that
should be used.
We employ this formulation in the study of the Fe−
and Fe2 dimers and the addition-elimination reaction
of molecular hydrogen on FeO+:
matic cases of the challenges for first principles meth-
ods to accurately reproduce the many low-lying mul-
tiplet potential energy surfaces associated with transi-
tion metals. It has been argued that spin density func-
tional theory can describe the lowest lying state of a
given spatial and spin symmetry [13, 14], but difficul-
ties remain in obtaining accurate multiplet splittings.
Our GGA or GGA+U calculations have been per-
formed with Quantum-ESPRESSO ; coupled clus-
ter (CCSD(T)) and B3LYP calculations have been per-
formed with Gaussian03 .
The iron dimer has been investigated both theoret-
ically [4, 19, 20, 21] and experimentally [22, 23, 24].
i)2can be interpreted as an effective
j). Even if in
i= 1 and
these are paradig-
FIG. 1: Linear-response Uout calculated from the GGA+Uin
ground state of7∆uFe2, together with the extrapolated Uscf.
U0 is Uout calculated for Uin=0.
The experimental photoelectron spectrum of Fe−
low 2 eV is remarkably simple - there are only two
prominent peaks, one at 1.0 eV and a second peak
0.53 eV above it, corresponding to two allowed transi-
tions to different neutral Fe2states. A recent multi-
reference configuration-interaction (MRCI) study has
assigned the three experimental electronic states involved
CCSD(T) has been shown to be in overall agreement.
Importantly, these electronic states are consistent with
the experimental measurements for the anion (funda-
mental frequency ω0=250±20 cm−1and bond length
Re=2.10±0.04 ˚ A ), and the two neutral Fe2 states,
which display similar properties (ω0=300±15 cm−1and
Re=2.02±0.02˚ A ).
We first apply our approach to Fe2and Fe−
a U0of 2 eV (i.e. when calculated from the GGA ground
state) and a Uscf of 3 eV (since energies at different U
are not directly comparable, we average U0and Uscfover
all states considered). GGA+Uscf shows a striking and
consistent agreement with MRCI and our CCSD(T)
results, correctly identifying both the lowest anion state
cited state,8∆g, 0.38 eV above. The lowest, singly ion-
ized neutral states, which differ from Fe−
loss of the spin down or spin up σ∗
gGGA+Uscf splitting of 0.6
eV compares very well with theoretical results (MRCI
and CCSD(T)) and the experimental splitting (0.53 eV)
in Table I. The structure of these two states (see Table
II) is also consistent with experimentally observed close
similarity of Reand ω0for the two neutral states and the
modest decrease in Re (0.08˚ A ) and increase in ω0 (∼
50 cm−1) with respect to Fe−
GGA+Uscf, GGA favors the
tive to other methods. Neutral states arising from single
ionization of the8∆g state are7∆u (3d144s2) and9∆g
gfor Fe2; more recently,
2. We obtain
u) and the first ex-
2only by the
u(4s) orbital, are9Σ−
u) by as much as 0.9 eV rela-
u, 4s3: σ2
State B3LYP GGA +U0 +Uscf CCSD(T) MRCIa
0.14 -0.52 0.04
TABLE I: Multiplet splittings (in eV) for Fe−
several levels of theory. (a) Ref. .
2 and Fe2 at
GGA+Uscf CCSD(T) MRCIa
2.20, 3012.24, 276 2.23, 272 2.1, 250
2.08, 355 2.12, 321 2.4, –
8∆g 2.07, 360
2.17, 296 2.18,299
2.16, 304 2.17,310
2.00, 404 2.25,195
2.28, 220 2.35, –
7∆u 1.99, 413
9∆g 2.26, 285
TABLE II: Bond lengths (˚ A ) and harmonic frequencies,
ωe, (cm−1) for Fe−
fundamental frequencies, ω0). (a) Ref. . (b) Ref. .
2and Fe2, compared to experiment (here,
(3d134s3) which result from the loss of σ∗
electrons, respectively. In addition, these two states have
differing bond lengths (Reof 1.99 and 2.26˚ A ) and fre-
quencies (ωeof 413 cm−1and 285 cm−1), and thus are
not compatible with experiment[4, 23].
Our second test case explores the potential energy
surfaces of the highly exothermic (∆H < −1.6 eV )
addition-elimination reaction of molecular hydrogen on
bare FeO+. This spin-allowed reaction occurs with ex-
ceedingly low efficiency (1 in every 100-1000 gas-phase
collisions results in products), yet when it does proceed
it is observed to be barrierless[25, 26, 27].
parent contradiction has been explained by a two-state-
reactivity model [5, 28, 29], wherein the steep reaction
barriers along the spin surface of the reactants and prod-
ucts (sextets in both cases) preclude an efficient, exother-
mic reaction. Instead, the reaction must occur along
a shallow but excited spin surface (here, the quartet),
and the reaction bottleneck is the coupling of the two
surfaces which permits the necessary spin-inversion at
the entrance and exit channels. For several exchange-
correlation functionals, (including B3LYP)[5, 29], the re-
action coordinates have failed to agree qualitatively with
experiments[25, 26, 27], higher level correlated-electron
calculations[28, 30], or with the established paradigm of
a two-state model.
For the bare FeO+reactant, GGA predicts a
ground state and two nearly degenerate low-lying quartet
states,4∆ and4Φ, 0.84 eV above. GGA+Uscf (5.5 eV)
preferentially stabilizes4Φ FeO+and yields a6Σ+→4Φ
u(4s) and σg(3d)
GGA+U 1.66 749 432
CCSD(T) 1.66 724 434
1.62 901 328
1.56 1038 332
TABLE III: Equilibrium bond lengths, Re (˚ A), harmonic fre-
quencies, ωe (cm−1), and anharmonicities, ωexe (cm−1) for
the6Σ+and4Φ states of FeO+.
FeO++ H2 reaction using GGA (blue) as compared against
a CCSD(T) reference (black).
dashed indicates quartet. (Color online.)
Potential energy surface and geometries for the
Solid indicates sextet while
splitting of 0.54 eV in quantitative agreement with the
symmetry and splitting (0.57 eV) predicted by CCSD(T).
The U correction also reduces the 3d character of minor-
ity spin π molecular orbitals which dramatically improves
bond lengths, harmonic frequencies, and anharmonici-
ties, as shown in Table III.
We thus proceed to study the full sextet and quar-
tet potential energy surfaces (PES) for this reaction.
We stress that, as is commonly found for open-shell
transition-metal molecules, several low-lying PES exist
for each multiplicity and we present results for the lowest-
lying symmetry of each multiplicity. The Uscfapplied in
this global PES is 5 eV, very close to the average of the
Uscf (4.93 eV) calculated for the quartet (5.02 eV) and
sextet (4.84 eV) at each stationary point; the values of
U0 are similar (quartet = 4.71 eV; sextet = 4.76 eV).
Although most states possess a Uscf close to the global
average, the few deviations will be highlighted later.
Our GGA results for the intermediates (Int) and tran-
sition states (TS) along the reaction coordinate confirm
the previously noted failures. Aside from the overesti-
mate of FeO+splittings, the most notable deviations are
unusually steep barriers (0.54 eV) along the quartet sur-
face, lack of spin-crossing near the products, and a a dra-
matic underestimate in the exothermicity, as depicted in
FIG. 3: Potential energy surface and geometries for the FeO+
+ H2 reaction using GGA+U (5 eV) (blue) as compared
against a CCSD(T) reference (black). Solid indicates sextet
while dashed indicates quartet, as in Fig. 2. (Color online.)
∆E6→4 GGA GGA+U CCSD(T)
TABLE IV: Multiplet splittings (in eV) using GGA, GGA+U
(U = 5 eV except in parentheses, UInt−3,av= 3.5 eV) and
With GGA+U (5 eV), we obtain consistency with
CCSD(T), as shown in Fig. 3. The reactant FeO+split-
ting is reduced, the splitting at Int-1 increases, corre-
sponding to a shallow quartet reaction coordinate, and
the exothermicity and spin crossover near the products
are in good agreement with experiment and theoretical
paradigm. The quantitative accuracy of GGA+U be-
comes fully evident in the intermediate splittings (Table
IV), forward and back reaction barriers (Table V), and
overall mean absolute errors (MAE) in multiplet split-
tings that are reduced (with respect to CCSD(T) ref-
erence) from 0.20 eV for GGA to 0.04 eV for GGA+U.
Geometries are also improved: the MAE for bond lengths
Forward Reaction Back Reaction
∆Ea GGA GGA+U CCSD(T) GGA GGA+U CCSD(T)
TS-261.22 0.82(1.16) 1.11
TABLE V: Comparison of GGA, GGA+U (U = 5 eV except
in parentheses, U4s = U3d = 4 eV) and CCSD(T) forward
and back reaction barriers (in eV).
are reduced from 4.3 pm (GGA) to 2.2 pm (GGA+U).
. The GGA+U and CCSD(T) states also possess con-
sistent orbital occupations and symmetry.
The few examples of Uscfdeviating from 5 eV are pri-
marily at the exit channel, where large changes in hy-
bridization occur. The quartet Int-3 is the only case for
which we obtain a low Uscf(2 eV) which originates from
the reduced hybridization of Fe 3d states. We chose to
recalculate the splitting with a Uscf,av that was a lo-
cal average on the Int-3 states. With this U of 3.5 eV,
we obtain a splitting of 0.12 eV, in even closer agree-
ment with CCSD(T). While this reduced hybridization
of the 3d states is unusual, we stress that it is consis-
tently predicted in our linear-response approach. Along
the sextet surface, the iron valence occupations corre-
spond to 3d64s1, and we find that the the interplay of 3d
and 4s states to be critical for describing the second bar-
rier along the sextet reaction surface. A matrix extension
of our formalism considers also the response of the 4s
orbitals, and we obtain U4s,scf= 4.0 eV and U3d,scf=4.0
eV around the barrier (U4sis instead found to be nearly
zero elsewhere). Inclusion of the 4s response for both sex-
tets Int-2 and Int-3 increases the forward reaction barrier
to 1.16 eV while the backward barrier remains unchanged
- in accordance with CCSD(T).
In conclusion, we have shown how a self-consistent
GGA+U approach can provide a dramatic improvement
to the description of multiplet potential energy surfaces
for transition-metal complexes that are otherwise poorly
described by common exchange-correlation functionals,
while preserving the very favorable computational costs
and scaling of local density-based functionals.
improvements include spin energetics, state symmetries,
and quantitative description of complex reaction coordi-
nates. U has been treated as an intrinsic, non-empirical
property of the system considered, and never as a fit-
ting parameter, and it has been obtained through a self-
consistent extension to the linear-response formulation
of Cococcioni and de Gironcoli. Such development
will allow large-scale and accurate calculations on
transition-metal complexes, with applications in the field
of catalysis, biochemistry, and environmental science.
We thank F. de Angelis for pointing out the H2 on
FeO+reaction and S. de Gironcoli for helpful discussions
on Uscf. This work was supported by an NSF graduate
fellowship and ARO-MURI DAAD-19-03-1-0169. Com-
putational facilities were provided through NSF grant
DMR-0414849 and PNNL grant EMSL-UP-9597.
∗Current Address: Departamento de Qu´imica Inorg´ anica,
Anal´itica y Qu´imica F´isica, Universidad de Buenos Aires
 A. C. Rosenzweig et al., Nature 366, 537 (1993).
 R. R. Schrock and J. A. Osborn, J. Am. Chem. Soc. 98,
 P. Jordan et al., Nature 411, 909 (2001).
 D. G. Leopold et al., J. Chem. Phys. 88, 3780 (1988).
 S. Shaik and M. Filatov, J. Phys. Chem. A 102, 3835
 L. Cavallo and H. Jacobsen, J. Phys. Chem. A 107, 5466
 C. Rovira et al., Biophys. J. 81, 435 (2001).
 J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.
Lett. 77, 3865 (1996).
 V. I. Anisimov, J. Zaanen, and O. K. Andersen, Phys.
Rev. B 44 943(1991).
 F. Zhou et al., Phys. Rev. B 70, 235121 (2004).
 A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen Phys.
Rev. B 52, R5467 (1995).
 M. Cococcioni and S. de Gironcoli, Phys. Rev. B 71,
 T. Ziegler, A. Rauk and E. J. Baerends, Theor. Chim.
Acta. 43, 261 (1977).
 O. Gunnarsson and B. I. Lundqvist, Phys. Rev. B 13,
4274(1976); U. von Barth, Phys. Rev. A 20, 1693(1979).
 Multiplets are defined by their spin component along
the z-axis and are thus not eigenstates of the square of
the spin operator. Consequences and alternatives are dis-
cussed more fully in M. Filatov and S. Shaik, J. Chem.
Phys. 110, 116(1999) and references therein.
 S. Baroni et al., http://www.quantum-espresso.org. Cal-
culations are completed in the Perdew Burke Ernzer-
hof (GGA) approximation using ultrasoft pseudopo-
tentials with a plane wave cutoff of 40 Ry and density
cutoff of 480 Ry; transition states are obtained with the
nudged elastic band method.
 G. Henkelman, B. P. Uberuaga, and H. Jonsson, J. Chem.
Phys. 113, 9901 (2000).
 M. J. Frisch et al., Gaussian, Inc. (2004). Hybrid func-
tional B3LYP (Becke’s 3-parameter exchange and Lee,
Yang, and Parr’s correlation) and CCSD(T) calculations
use the 6-311++G(3df,3pd) basis set.
 S. Dhar and N. R. Kestner, Phys. Rev. A. 38, 1111
 O. Hubner and J. Sauer, Chem. Phys. Lett. 358, 442
(2002) and references therein.
 A. Irigoras et al., Chem. Phys. Lett. 376, 310 (2003).
 T. McNab, H. Micklitz, and P. Barrett, Phys. Rev. B 4,
 D. Leopold and W. Lineberger, J. Chem. Phys 85, 51
 E. A. Rohlfing et al., J. Chem. Phys. 81, 3846 (1984).
 D. E. Clemmer et al., J. Phys. Chem. 98, 6522 (1994).
 D. Schroeder et al., Int. J. Mass. Spec. 161, 175 (1997).
 J. M. Mercero et al., Int. J. Mass. Spec. 240, 37 (2005)
and references therein.
 A. Fiedler et al., J. Am. Chem. Soc. 116, 10734 (1994).
 D. Danovich and S. Shaik, J. Am. Chem. Soc. 119, 1773
 A. Irigoras, J. Fowler, and J. Ugalde, J. Am. Chem. Soc.
121, 8549 (1999).
 P. Giannozzi, F. De Angelis, and R. Car, J. Chem. Phys.
120, 5903 (2004).
 The ωeare obtained with a power law fit following Hollas,
High Resolution Spectroscopy, Wiley, 1998. Anharmonic-
ity prevents direct comparison to the experimental ω0.
 The PES of Figs. 2 and 3 have been aligned at6Int-1.
 CCSD(T) geometries are from an 0.01˚ A interpolation in
each degree of freedom (up to five) in total as much as
800 for a structure. The states considered are reactants
and intermediates for each quartet and sextet PES.