Skip to content

DFT+U (Hubbard Correction)

The DFT+U method adds an on-site Hubbard correction to selected atomic subshells, penalizing fractional occupation and partially correcting the self-interaction error of (semi-)local functionals. It is most useful for systems with localized \(d\) or \(f\) electrons (transition-metal and rare-earth compounds) where standard DFT over-delocalizes the open-shell density.

Theory

Approximate (semi-)local exchange–correlation functionals suffer from self-interaction error: an electron spuriously repels itself, which biases the description toward over-delocalized, partially-occupied orbitals. For localized \(d\)/\(f\) shells this shows up as qualitatively wrong physics — metallic instead of insulating ground states, too-small magnetic moments, mis-ordered oxidation states. DFT+U adds an explicit energy penalty for fractional on-site occupation, nudging the corrected subshells back toward integer filling and restoring much of the localized character at negligible cost.

TeraChem implements the rotationally-invariant, single-parameter (Dudarev-type) form of the correction, which can be layered on top of any Hartree–Fock or exchange–correlation functional. The corrected energy is

\[ E^{\mathrm{HF/DFT}+U} = E^{\mathrm{HF/DFT}} + \frac{1}{2} \sum_{I,\sigma} \sum_{nl} U_{nl}^{I} \left[\operatorname{Tr}\big(\mathbf{q}_{nl}^{I,\sigma}\big) - \operatorname{Tr}\big(\mathbf{q}_{nl}^{I,\sigma}\,\mathbf{q}_{nl}^{I,\sigma}\big)\right], \]

where \(I\) runs over atoms, \(\sigma\) over spin, and \(nl\) over the subshell to which \(U_{nl}^{I}\) is applied. The penalty \(\operatorname{Tr}[\mathbf{q}(1-\mathbf{q})]\) vanishes when the subshell occupation eigenvalues are exactly \(0\) or \(1\) (empty or fully occupied) and is largest at half-filling, so a positive \(U\) energetically favors integer occupations.

The occupation matrix \(\mathbf{q}_{nl}^{I,\sigma}\) is spin-resolved and, by default, is built from Mulliken populations — the entrywise (Hadamard) product of the density and AO-overlap matrices, \(q_{ab} = P_{ab}\,S_{ab}\), restricted to the basis functions of the corrected subshell. The dft+u_project options instead build it from a Löwdin-orthogonalized density, \(\mathbf{q} = \mathbf{S}^{1/2}\mathbf{P}\,\mathbf{S}^{1/2}\) (optionally projected onto a separate basis). The correction enters the Fock matrix and is solved self-consistently, so it works for RHF, ROHF, and UHF references alike. \(U\) (and \(\alpha\)) are entered in eV.

A linear potential \(\alpha_{nl}^{I}\) can be applied to the same subshell, in lieu of or together with \(U\),

\[ E^{\mathrm{lin}} = \sum_{I,\sigma} \sum_{nl} \alpha_{nl}^{I}\,\operatorname{Tr}\big(\mathbf{q}_{nl}^{I,\sigma}\big). \]

Its main use is the linear-response determination of \(U\) following Cococcioni and de Gironcoli: \(U\) is recovered as the difference between the bare (first-SCF-iteration) and self-consistent responses of the subshell occupation to the perturbation \(\alpha\), \(U = d\alpha/d N_0 - d\alpha/dN\).

Analytic gradients and minimal basis sets

The dominant part of the +U nuclear gradient — the +U potential acting through the density matrix — is captured automatically by the self-consistent SCF. The only additional contribution is an overlap-derivative term (a correction to the energy-weighted density matrix) that vanishes whenever the same-subshell basis functions on an atom are orthonormal, as they are in a minimal basis set. For the minimal-basis incompleteness correction described below, analytic gradients are therefore complete and geometry optimizations and dynamics are well defined — this is how Kulik et al. optimize protein geometries with STO-3G+U. With larger, non-minimal basis sets that extra term is nonzero and is not currently added, so forces may be slightly inconsistent with the +U energy in that regime.

Alternative use: correcting minimal basis set incompleteness

The conventional motivation above is self-interaction error in localized \(d\)/\(f\) shells. The same dft+u machinery can also be repurposed for an unrelated problem — the incompleteness error of minimal basis sets — and this application was developed for and validated in TeraChem (Kulik, Seelam, Mar, and Martínez, 2016).

Minimal basis sets such as STO-3G are attractive for very large systems (proteins of \(\sim\!10^3\) atoms on GPUs), but their incompleteness produces qualitatively wrong chemistry: most notably spurious proton transfer from nitrogen to oxygen along amide backbones, and incorrect tautomer orderings, because the minimal basis imbalances the relative electron/proton affinities of nitrogen and oxygen substituents. Applying a tuned \(U\) to the 2p shells of nitrogen and oxygen alters their effective hybridization and orbital filling — i.e. it modifies matrix elements to compensate for the missing basis flexibility — and restores the balance, so that a cheap minimal-basis Hartree–Fock calculation reproduces the geometries and relative energies of a much more expensive polarized double-\(\zeta\) (e.g. 6-31G*) calculation.

Because this is an empirical correction rather than an SIE correction, the \(U\) values are tuned rather than physically derived, and — unusually — can be negative. The representative values below were fit to cytosine tautomers and shown to transfer well across formamide, guanine, thymine, and a 10-protein test set:

Atom (2p shell) \(U\) (eV)
Nitrogen \(+6\)
Oxygen \(-6\)

The positive \(U\) on nitrogen and negative \(U\) on oxygen apply opposing shifts that amplify the correction. With these values, "RHF/STO-3G+U" reduces anomalous proton-transfer events by roughly 90% (from 40 to 4 across the 10-protein set) and brings minimal-basis relative energies into qualitative — and often near-quantitative — agreement with RHF/6-31G*, at minimal-basis cost.

To set this up, apply \(U\) to the 2p shell of every nitrogen and oxygen through a Hubbard input file on top of an RHF/STO-3G calculation:

N
2P
6
O
2P
-6

coordinates    protein.xyz
basis          sto-3g
method         rhf
run            minimize
dft+u          yes
dft+u_default  zero
hubbard_input  hub.txt
end

A different question, the same keywords

This basis-set-incompleteness use of dft+u is distinct from the usual self-interaction / localized-electron application. The two share the same energy functional and keywords but answer different questions — do not conflate the published \(-6/+6\) eV values (a basis-set correction for N/O 2p shells) with \(U\) values appropriate for transition-metal \(d\) shells. If you need quantitative energetics, re-tune the parameters for your own chemistry; the published values were optimized for amide and nucleobase systems.

Turning DFT+U on

Set dft+u yes. By itself this enables the correction; the \(U\) values themselves are supplied in one of the two ways below. A linear potential \(\alpha\) can additionally (or instead) be applied — see Linear potential terms.

dft+u           yes
dft+u_default   zero        # how unspecified shells are treated
hubbard_input   hub.txt     # optional: custom per-shell U values

Specifying values of \(U\)

Values of \(U\) (in eV) can be assigned in one of two ways.

1. From the xyz file (default). Add a fifth column to the right of each atom's coordinates in the .xyz file giving that atom's \(U\) value. The value is applied to the default shell for that element. Atoms with no value listed either receive no \(U\) correction (default) or are set to a built-in default value when dft+u_default yes is used.

2. From a Hubbard input file. Set hubbard_input <file> to read custom shell assignments from a text file. Each atom type is listed sequentially, followed by the shell(s) to which a \(U\) value should be applied and the corresponding value(s) in eV. To give different values to different atoms of the same element, number the element names (e.g. C1) in both the Hubbard file and the xyz file.

The following Hubbard file and xyz file describe a water dimer in which \(U\) values are assigned to all hydrogen 1s shells and to the 2s and 2p shells of one of the two oxygen atoms:

O1
2S 2P
2 3
H
1S
1

6
  water dimer
O1   -0.702196054  -0.056060256   0.009942262
H    -1.022193224   0.846775782  -0.011488714
H     0.257521062   0.042121496   0.005218999
O     2.268880784   0.026340101   0.000508029
H     2.645502399  -0.412039965   0.766632411
H     2.641145101  -0.449872874  -0.744894473

Linear potential terms

A linear potential \(\alpha\) may be applied to a subshell — in lieu of, or in addition to, \(U\) — by setting dft+u_alpha yes. The linear potential must be given in a text file named by the alpha_input parameter, using the same layout as the Hubbard file: each atom type is listed sequentially, followed by the shell(s) and the corresponding \(\alpha\) value(s) in eV.

No defaults for \(\alpha\)

Unlike \(U\), there are no built-in defaults for the linear potential. If you enable dft+u_alpha yes you must supply alpha_input, or no linear potential will be applied.

Worked example

A 13-atom DFT+U single point (from the DftPlusU/dft_plus_U test). Here the master switch is on, unspecified shells default to zero, and custom shells are read from hub.txt:

coordinates     A.xyz
charge          0
run             energy
spinmult        1
basis           sto-3g
method          RHF
maxit           150
gpus            all
dft+u           yes
dft+u_default   zero
hubbard_input   hub.txt
end

O
2S 2P
0.0 0.00
N
2S 2P
0.0 4.00

TeraChem echoes the resolved per-atom shell assignments near the top of the output, e.g.:

DFT+U orbital information
3   N    ( 2S: 0.000000 )  ( 2P: 4.000000 )
5   O    ( 2S: 0.000000 )  ( 2P: 0.000000 )
6   N    ( 2S: 0.000000 )  ( 2P: 4.000000 )
12   N    ( 2S: 0.000000 )  ( 2P: 4.000000 )

and prints the resulting DFT+U occupation matrix for each corrected subshell.

Keyword reference

Keyword Values Default Description
dft+u yes / no no Master switch for the Hubbard correction.
dft+u_default no / zero / yes / built-in no Treatment of shells with no value specified: no/zero apply no correction; yes/built-in apply a built-in default \(U\).
hubbard_input filename — Text file with custom per-shell \(U\) values (eV). Overrides the default-shell behavior.
dft+u_alpha yes / no no Apply a linear potential \(\alpha\) to selected subshells.
alpha_input filename — Text file with the per-shell \(\alpha\) values (eV). Required when dft+u_alpha yes.
dft+u_project integer — Projection scheme used to build the occupation matrix. Values > 2 require dft+u_project_basis.
dft+u_project_basis basis name — Basis set used for the projection when dft+u_project > 2.
u_ramping yes / no / integer no Ramp \(U\) on over several SCF macro-iterations to aid convergence. yes uses 5 steps; an integer sets the number of steps.

References

  • Minimal-basis-incompleteness application and the TeraChem implementation (the alternative use of DFT+U described above — not the conventional SIE correction): H. J. Kulik, N. Seelam, B. D. Mar, and T. J. Martínez, "Adapting DFT+U for the Chemically Motivated Correction of Minimal Basis Set Incompleteness," J. Phys. Chem. A 120, 5939–5949 (2016). doi:10.1021/acs.jpca.6b04527
  • Original LDA+U method: V. I. Anisimov, J. Zaanen, and O. K. Andersen, "Band theory and Mott insulators: Hubbard \(U\) instead of Stoner \(I\)," Phys. Rev. B 44, 943 (1991). doi:10.1103/PhysRevB.44.943
  • Rotationally-invariant, single-parameter (\(U_\mathrm{eff}\)) functional implemented here: S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton, "Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+U study," Phys. Rev. B 57, 1505 (1998). doi:10.1103/PhysRevB.57.1505
  • Linear-response determination of \(U\) (the \(\alpha\) perturbation): M. Cococcioni and S. de Gironcoli, "Linear response approach to the calculation of the effective interaction parameters in the LDA+U method," Phys. Rev. B 71, 035105 (2005). doi:10.1103/PhysRevB.71.035105