DFT+\(U\)(+\(J\))
 Author:
David D. O’Regan, Trinity College Dublin
 Author:
Davide Sarpa, University of Southampton
 Date:
July 2015 (Revised June/July 2024)
DFT+\(U\) is fully and selfconsistently implemented in ONETEP, together with a number of advanced ancillary functionalities. The method is linearscaling with respect to system size, exhibiting no systematic tendency to slow convergence to the groundstate. DFT+\(U\) in its conventional fixedprojector form introduces only a small increase in computational prefactor with respect to the underlying exchangecorrelation functional [ORegan2012].
A very short introduction to DFT+\(U\)
DFT+\(U\) [Anisimov1991], [Anisimov1997], [Dudarev1998], also known as LDA+\(U\) or LSDA+\(U\), is a method used to improve the description of socalled strongly correlated materials offered by DFT within conventional approximations for exchangecorrelation (XC) such as the LSDA and \(\sigma\)GGA, quantitatively and even qualitatively. These functionals, based on the locallyevaluated density and its gradients, can sometimes fail to reproduce the physics associated with localised orbitals of \(3d\) and \(4f\) character characteristic of conventionallyclassed strongly correlated materials, a category consisting of not only firstrow transition metals and their oxides, but also lanthanoid oxide materials, and other materials such as certain magnetic semiconductors and organometallic molecules.
Typically, the LDA and its extensions have the tendency to favour highspin groundstates and to underestimate local magnetic moments and the band gap in cases where it is related to electron localisation, e.g., in chargetransfer or Mott insulators. These failures stem from the inability of the approximate xc functionals to satisfy the flatplane condition, which combines two constraints that the exact functional must meet: piecewise linearity with respect to electron number and the spin constancy constraint [Aron_J.Cohen2012]
Put simply, if the system under study comprises open \(3d\) or \(4f\) subshells, then there is a good chance that the LDA will find a minimum energy configuration by partly occupying them, leading to an underestimation of the insulating gap and any associated magnetic order.
In DFT+U, the idea is to describe the strongly correlated electronic states using the Hubbard model by projecting the KohnSham orbitals onto a localised basis, while the rest of the valence electrons are treated at the standard DFT level.
The general form for the total energy in DFT+\(U\) is given by:
In the above \(E_{\text{DFT}}[\rho(r)]\) is the energy at the standard DFT level, \(E_{\text{Hub}}\left[\{n_i\}\right]\) is the Hubbard correction and the last term \(E_{dc}\left[\{n_i\}\right]\) is a doublecounting correction for that contribution already included in the LDA term.
The DFT+\(U\) correction term is usually thought of as an explicit meanfield treatment of the exchangecorrelation energy contributed by the correlated sites (subspaces projected out with functions of \(3d\) and or \(4f\) character). There are many variants of the DFT+\(U\) method, the flavour implemented in ONETEP is the basisset independent, rotationally invariant quadratic penalty functional of Ref [Cococcioni2005], defined by the additive energy correction.
Here, \(U\) is an estimate of the scalar screened densitydensity Coulomb repulsion between localised orbitals. The occupancy matrix of the correlated atomic site \(I\), for spin channel \(\sigma\), is defined, in the case of orthonormal projector functions \(\lbrace \lvert \varphi^{(I)}_m \rangle \rbrace\), and densitymatrix \(\hat{\rho}^{(\sigma)}\), by
The DFT+\(U\) method forces the onsite occupancy matrix to be idempotent. This means that the DFT+\(U\) method penalises noninteger occupancy of these orbitals, tending to fill states with occupancy greater than \(0.5\) or occupancy less than \(0.5\), as can be seen from the expression for the DFT+\(U\) potential.
The DFT+\(U\) term may be considered as a correction which cancels the contribution to the energy arising due to the spurious selfinteraction of a partially occupied orbital [Cococcioni2005]. In this case, the \(U\) parameter is the curvature of the total energy with respect to the occupancy of the correlated manifold  which should be a piecewise linear curve were Janak’s theorem satisfied [Janak1978].
The \(U\) parameter can be computed empirically or from linearresponse theory (among other methods such as constrained DFT) according to the prescription given in Refs. [Cococcioni2005], [Kulik2006].
DFT+\(U\)+\(J\)
A corrective term to include \(J\) for antiparallel spin is described in Ref. [Himmetoglu2011] and is given by
Inclusion of this term leads to the so called DFT+\(U\)+\(J\) method, which is fully implemented in ONETEP and can be activated by assigning a value to the \(J\) term.
Computing \(U\) and \(J\) from linearresponse
The basic idea is to compare the response of the system to a perturbation in the DFT and in the DFT+\(U\) frameworks. We start by defining the response function \(\chi\), which describes how the occupation of localised orbitals changes with respect to a shift in the potential acting on these orbitals: The linear response method determines the Hubbard \(U\) parameter by comparing the response of the system to a perturbation in standard DFT and DFT+\(U\) frameworks.
We define the response function \(\chi\) as:
where \(n\) is the occupation matrix of the localised orbitals and \(\alpha\) is a potential shift applied to these orbitals.
We compute two response functions:
\(\chi_0\): the bare KohnSham (KS) response (without \(U\))
\(\chi\): the interacting response (with \(U\))
These are related by:
which allow us to compute \(U\).
In practice, we compute \(\chi_0\) and \(\chi\) by applying a small perturbation \(\alpha\) to the system:
The perturbation is applied by shifting the potential of the localised orbitals:
We then iterate until selfconsistency is achieved.
This is done in a supercell as the perturbation should not interact with its periodic images.
This is the conventional linear response but it poses practical problems: The response \(\chi_0\) is usually computed via the first iteration of the KohnSham equations during a selfconsistent field (SCF) calculation; that is, the response is to be measured following the initial charge redistribution introduced by the perturbation but before the KohnSham potential is updated. This approach is impractical to implement in codes that use a directminimization procedure of the total energy with respect to the density, KohnSham orbitals, or density matrix.
Another approach to compute \(U\) and \(J\) is known as minimum tracking method [Linscott2018].
Minimum Tracking Method
The minimum tracking method is based on a reformulation of the response matrices based on the ground state density of the perturbed system. We can identify the interacting and noninteracting response matrices as:
This allows us to work around the practical issues from the conventional linear response. This approach can also be extended to include the \(J\) exchange term (The response matrices now become rankfour tensors [Linscott2018]). In practice this is done by modifying the perturbation by including an additional term (spinsplitting):
Using NGWFs and projector selfconsistency
Any reasonable set of localised atomiclike functions may, in principle, be used for the projectors defining the correlated subspaces in DFT+\(U\); the choice is somewhat arbitrary and the description “atomic orbitals” does not uniquely define them. One possible approach is to use Wannier functions for the KohnSham orbitals, so that the correlated subspaces are proper subspaces of the KohnSham Hilbert space. Indeed, there is numerical evidence to suggest that Maximally Localised Wannier Functions (MLWFs) [Marzari1997], [Souza2001], in particular, provide a basis that maximises a particular measure of the onsite Coulomb repulsion [Miyake2008], and MLWFs are in common use as a minimal basis with which to construct tightbinding models from firstprinciples.
In ONETEP, a set of variationallyoptimised nonorthogonal generalised Wannier functions (NGWFs) are generated as a byproduct of totalenergy minimisation. NGWFs exhibit some similar properties to MLWFs and other flavours of localised Wannier functions, and, for example, can be used to calculate finitedifference response properties in a similar way [ORegan20122]. As they are conveniently available in ONETEP, we have made it possible to reuse the NGWFs from the end of a groundstate calculation as a set of Hubbard projectors with which to define the DFT+\(U\) correction. For this, it was necessary to develop a tensoriallyconsistent formulation of DFT+\(U\) in order to accommodate nonorthogonal projector functions [ORegan2011]; projector nonorthogonality for a given subspace is automatically compensated for.
In order to ensure that NGWFs with appropriate symmetry are chosen as
Hubbard projectors for a given atom, those \(n\) NGWFs
\(\lvert \phi_\alpha \rangle\) that maximise \(\sum^n_{m,\alpha }\langle \varphi_m \rvert \phi^\alpha \rangle \langle \phi_\alpha \rvert \varphi_m \rangle\), for a given set of
\(n\) hydrogenic orbitals \(\lvert \varphi_m \rangle\), defined
in the hubbard
block, are selected for the task.
The keyword hubbard_max_iter
, (defaulting to \(0\)), sets the task to
HUBBARDSCF
, which performs a selfconsistency cycle over the Hubbard
projectors, demonstrated in
Refs. [ORegan2010], [ORegan2011].
The density from one minimisation is reused at the beginning of the next, and setting
hubbard_max_iter
to \(2\) one can carry out a DFT+\(U\)
calculation using the LDA NGWFs as projectors.
The keywords hubbard_energy_tol
, hubbard_conv_win
, and
hubbard_proj_mixing
are used to manage the Hubbard projector
selfconsistency cycle. For convergence, the ground state energy must
deviate less than hubbard_energy_tol
(defaulting to
\(10^{8}\) Ha) from one HUBBARDSCF
iteration to the next, over
hubbard_conv_win
(defaulting to \(2\)) iterations. A fraction
hubbard_proj_mixing
(defaulting to \(0.0\)) of the previous
Hubbard projectors may be mixed with the new ones in order to accelerate
the procedure, although this has never been found to be necessary.
Setting hubbard_proj_mixing
to a negative value causes the
projectors to be read in from a .tightbox_hub_projs
file, for
restarting a HUBBARDSCF
calculation or for a variety of
postprocessing tasks.
How to activate DFT+\(U\) in ONETEP
In order to activate the DFT+\(U\) functionality, the hubbard
block is added to the input file. For example, in the case of a system
containing iron and cerium atoms incorrectly described by the
exchangecorrelation functional, which we suspect could benefit from the
DFT+\(U\) correction to improve the description of localisation,
we might use the hubbard
block:
%block hubbard
Fe1 2 4.0 0.0 10.0 0.00 1.0
Fe2 2 4.0 0.0 10.0 0.00 1.0
Ce1 3 6.0 0.0 10.0 0.50 0.0
%endblock hubbard
The columns of the hubbard
block are described as follows:
Species Label
The species to apply the DFT+\(U\) correction to. In this example Fe1, Fe2 and Ce1.
Angular Momentum: \(l\)
The angular momentum of the projectors which the Hubbard correction is applied to. In this example \(l=2\) for Fe1 and Fe2 and \(l=3\) for Ce1. Conventionally, the radial quantum number \(r=l+1\) is used to generate atomcentred atomic projectors, so that \(l=2\) gives \(3d\) orbitals, \(l=3\) gives \(4f\) orbitals etc.
(please get in contact if you need to use a \(r \ne l+1\) combination, or multiple subshells per atom).
Hubbard \(U\) value
The value of the Hubbard \(U\) for this subshell, in electronvolts.
Hund’s exchange \(J\) value
The value of the Hund’s exchange \(J\) for this subshell, in electronvolts. The rotationally invariant exchange corrective term described in detail in Ref. [Himmetoglu2011] (The so called DFT+\(U\)+\(J\)) is fully implemented in ONETEP (including forces etc), and activated for any \(J \ne 0\).
Effective Charge \(\mathbf{Z}\) and Projectors type
Case 1: \(\mathbf{ Z < 0}\) (Default)
A subset of the orbitals generated by solving the atomic problem subject to the pseudopotential for the species in question is chosen (in which case the projectors form a subset of the initial guesses for the ONETEP NGWFs); here the magnitude of the negative Z makes no difference.
Case 2: \(\mathbf{ Z > 0}\)
The projectors are generated in the form of solutions to the hydrogenic Schrödinger equation. In this case \(\mathbf{Z}\) is the effective charge divided by the ratio of effective masses used to generate projectors. A good guess for this number might be the ClementiRaimondi effective charge, tabulated in Refs. [Clementi1963], [Clementi1967], and the choice of radial profile does matter [ORegan2010].
In both cases, the projectors are effectively renormalised within an atomcentred sphere with the same radius as the NGWFs on that atom.
The \(\alpha\) prefactor
An additional potential acting on the subspace in question, the prefactor \(\alpha\) is here entered in electronvolts. This is needed, for example, in order to locally vary the potential in order to determine the value of \(U\) which is consistent with the screened response in the system with linearresponse theory [Cococcioni2005], [Kulik2006], or to break a spatial symmetry, such as in a mixedvalence system. In the example given, we are additionally penalising the occupancy on cerium \(4f\) atomic orbitals.
The spinsplitting factor
The spinsplitting factor, in electronvolts, which is deducted from the \(\alpha\) factor for the spinup channel and added to \(\alpha\) for the spindown channel. In the example shown here we’re promoting spinup magnetisation for iron atoms Fe1, and spindown for Fe2. This can be very useful for appropriately breaking magnetic symmetries in antiferromagnetic solids or openshell singlet molecules, or for estimating the magnetic susceptibility or exchange coupling.
N.B. Users may find the DFT+\(U\) functionality useful in cases of systems even when the DFT+\(U\) correction is not needed (setting the all \(U\) parameters to zero does not disable the functionality). The implementation offers a very inexpensive method for carrying out carefullydefined atomcentred atomic population analysis, or breaking symmetries in spin or charge ordered systems.
DFT+\(U\) keywords
Keyword 
Type 
Default 
Description 


Task 
— 
Activate a projectorselfconsistent
DFT+\(U\) calculation.


Logical 
False 
Activate a nonvariational onthefly form of
projector selfconsistency in DFT+\(U\)
or cDFT, in which the projectors are updated
whenever the NGWFs are.
task:
HUBBARDSCF is then not needed. 

Integer 

The minimum number of Hubbard projector update
steps satisfying the incremental energy
tolerance
HUBBARD_ENERGY_TOL requiredfor convergence in task :
HUBBARDSCF . 

Physical 

The maximum incremental energy change between
Hubbard projector update steps allowed for
convergence in task :
HUBBARDSCF . 

Integer 

The form of DFT+\(U\) energy term used.
Contact developers if you need to try something
beyond the default.


Integer 

The maximum allowed number of Hubbard projector
update steps taken in a projector
selfconsistent DFT+\(U\) or cDFT
calculation in task:
HUBBARDSCF 

Physical 

The incremental change in energy, in
totalenergy minimisation, at which any
spinsplitting (Zeeman) type term in DFT+U is
switched off, and the minimisation history
reset. Useful for breaking openshell,
antiferromagnetic, or chargedensity
wave symmetries.


Real 

The fraction of previous Hubbard projector to
mix with new for projector selfconsistent
DFT+\(U\) or cDFT in task :
HUBBARDSCF .Not found to be necessary.


Logical 

Read Hubbard projectors from tightbox_hub_projs
file in restart calculations involving
DFT+\(U\).


Integer 

The form of correction used to correct for any
nonorthogonality between Hubbard projectors.
Contact developers if you need to try something
other than the default “tensorial” correction.

Tutorials
 Example on the use of DFT+\(U\) for Hematite, a strongly correlated system.
 Example on how to compute \(U\) and \(J\) from linear response.
To be added
Compatibility
The DFT+\(U\) functionality is fully compatible with almost all other parts of the ONETEP code, such as listed below, since it simply involves an additional term in the Hamiltonian and ionic forces. Please get in touch first if you would like to use a more exotic combination of these functionalities:
Totalenergy minimisation and ionic forces.
Geometry optimisation, molecular dynamics and phonon calculations.
All other functionals including hybrids and Van der Waals functionals.
Implicit solvation.
The PAW formalism and ultrasoft pseudopotentials.
Constrained DFT.
Local density of states (including a correlated subspace decomposition).
Natural bond orbital calculations.
Conductionband optimisation and Fermi’s Golden Rule spectra.
Calculations of changes in electric polarisation.
Timedependent DFT.
Electronic transmission calculations.
The extension of the DFT+\(U\) implementation to cluster Dynamical meanfield theory has also been implemented in ONETEP; for an example of its capabilities see Ref. [Weber2012].
D. D. O’Regan, N. D. M. Hine, M. C. Payne and A. A. Mostofi, Phys. Rev. B 85, 085107 (2012). https://doi.org/10.1103/PhysRevB.85.085107
J. Z. V. I. Anisimov and O. K. Andersen, Phys. Rev. B 44, 943 (1991). https://doi.org/10.1103/PhysRevB.44.943
V. I. Anisimov, F. Aryasetiawan, and A. I. Liechtenstein, J. Phys.: Condens. Matter 9, 767 (1997). https://iopscience.iop.org/article/10.1088/09538984/9/4/002
S. L. Dudarev, Phys. Rev. B 57, 3 (1998). https://doi.org/10.1103/PhysRevB.57.1505
A. J. Cohen, P. MorSanchez and W. Yang, Chem. Rev. 2012, 112, 289–320. https://doi.org/10.1021/cr200107z
M. Cococcioni and S. de Gironcoli, Phys. Rev. B 71, 035105 (2005). https://doi.org/10.1103/PhysRevB.71.035105
J. F. Janak, Phys. Rev. B 18, 12 (1978). https://doi.org/10.1103/PhysRevB.18.7165
H. J. Kulik, M. Cococcioni, D. A. Scherlis and N. Marzari, Phys. Rev. Lett. 97, 103001 (2006). https://doi.org/10.1103/PhysRevLett.97.103001
B. Himmetoglu, R. M. Wentzcovitch, and M. Cococcioni, Phys. Rev. B,84, 115108 (2011). https://doi.org/10.1103/PhysRevB.84.115108
E. Clementi and D.L. Raimondi, J. Chem. Phys. 38, 2686 (1963). https://doi.org/10.1063/1.1733573
E. Clementi, D.L. Raimondi, and W.P. Reinhardt, J. Chem. Phys. 47, 1300 (1967). https://doi.org/10.1063/1.1712084
D. D. O’Regan, N. D. M. Hine, M. C. Payne and A. A. Mostofi, Phys. Rev. B 82, 081102 (2010). https://doi.org/10.1103/PhysRevB.82.081102
C. Weber, D. D. O’Regan, N. D. M. Hine, M. C. Payne, G. Kotliar and P. B. Littlewood, Phys. Rev. Lett. 108, 256402 (2012). https://doi.org/10.1103/PhysRevLett.108.256402
N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 (1997). https://doi.org/10.1103/PhysRevB.56.12847
I. Souza, N. Marzari and D. Vanderbilt, Phys. Rev. B 65, 035109 (2001). https://doi.org/10.1103/PhysRevB.65.035109
T. Miyake and F. Aryasetiawan, Phys. Rev. B 77, 085122 (2008). https://doi.org/10.1103/PhysRevB.77.085122
D. D. O’Regan, M. C. Payne, and A. A. Mostofi, Phys. Rev. B 85, 193101 (2012). https://doi.org/10.1103/PhysRevB.85.193101
D. D. O’Regan, M. C. Payne and A. A. Mostofi, Phys. Rev. B 83, 245124 (2011). https://doi.org/10.1103/PhysRevB.83.245124
E.B. Linscott, D. J. Cole, M. C. Payne, D. D. O’Regan, Phys. Rev. B 98, 235157 (2018). https://doi.org/10.1103/PhysRevB.98.235157