| Author: | David D. O'Regan, Trinity College Dublin |
|---|---|
| Author: | Davide Sarpa, University of Southampton |
| Date: | July 2015 (Revised June/July 2024) |
Table of Contents
DFT+U is fully and self-consistently implemented in ONETEP, together with a number of advanced ancillary functionalities. The method is linear-scaling with respect to system size, exhibiting no systematic tendency to slow convergence to the ground-state. DFT+U in its conventional fixed-projector form introduces only a small increase in computational pre-factor with respect to the underlying exchange-correlation functional [O-Regan2012].
DFT+U [Anisimov1991], [Anisimov1997], [Dudarev1998], also known as LDA+U or LSDA+U, is a method used to improve the description of so-called strongly correlated materials offered by DFT within conventional approximations for exchange-correlation (XC) such as the LSDA and \sigma-GGA, quantitatively and even qualitatively. These functionals, based on the locally-evaluated density and its gradients, can sometimes fail to reproduce the physics associated with localised orbitals of 3d and 4f character characteristic of conventionally-classed strongly correlated materials, a category consisting of not only first-row 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 high-spin ground-states and to underestimate local magnetic moments and the band gap in cases where it is related to electron localisation, e.g., in charge-transfer or Mott insulators. These failures stem from the inability of the approximate xc functionals to satisfy the flat-plane 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 sub-shells, 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 Kohn-Sham 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:
E_{\text{DFT+U}}[\rho(r)] = E_{\text{DFT}}[\rho(r)] + E_{\text{Hub}}\left[\{n_i\}\right] - E_{\text{dc}}\left[\{n_i\}\right]
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 double-counting correction for that contribution already included in the LDA term.
The DFT+U correction term is usually thought of as an explicit mean-field treatment of the exchange-correlation 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 basis-set independent, rotationally invariant quadratic penalty functional of Ref [Cococcioni2005], defined by the additive energy correction.
E_{\text{DFT+U}} \left[ n^{(I) (\sigma)} \right] = \sum_{I \sigma} \frac{U^{(I)}}{2} \rm{Tr}
\left[ n^{(I) (\sigma)} \left( 1 - n^{(I) (\sigma)} \right)\right].
Here, U is an estimate of the scalar screened density-density 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 density-matrix \hat{\rho}^{(\sigma)}, by
n^{(I)(\sigma)}_{m m'} = \langle \varphi_m^{(I)} \rvert \hat{\rho}^{(\sigma)}
\lvert \varphi_{m'}^{(I)} \rangle.
The DFT+U method forces the on-site occupancy matrix to be idempotent. This means that the DFT+U method penalises non-integer 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.
\hat{V}^{(\sigma)}_{\text{DFT+U}} = \sum_{I} U^{(I)}
\lvert \varphi_m^{(I)} \rangle
\left( \frac{1}{2} \delta_{m m'} - n^{(I) (\sigma)}_{m m'} \right) \langle
\varphi_{m'}^{(I)} \rvert .
The DFT+U term may be considered as a correction which cancels the contribution to the energy arising due to the spurious self-interaction 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 piece-wise linear curve were Janak’s theorem satisfied [Janak1978].
The U parameter can be computed empirically or from linear-response theory (among other methods such as constrained DFT) according to the prescription given in Refs. [Cococcioni2005], [Kulik2006].
A corrective term to include J for anti-parallel spin is described in Ref. [Himmetoglu2011] and is given by
E_J[n^\sigma] = \sum_I \frac{J^I}{2} \text{Tr}[n^{I\sigma} n^{I-\sigma}]
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.
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:
\chi = \frac{dn^{I\sigma}}{d\alpha}
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 Kohn-Sham (KS) response (without U)
- \chi: the interacting response (with U)
These are related by:
U = \chi^{-1} - \chi_0^{-1}
which allow us to compute U.
In practice, we compute \chi_0 and \chi by applying a small perturbation \alpha to the system:
\chi_0 \approx \frac{\Delta n}{\Delta \alpha}
\text{(computed without U)}
\chi \approx \frac{\Delta n}{\Delta \alpha}
\text{ (computed with a trial U)}
The perturbation is applied by shifting the potential of the localised orbitals:
V_{\text{ext}}^{p} = V_{\text{ext}} + \alpha \sum_{m,m'}\lvert\varphi_{m'}^{(I)}\rangle\langle\varphi_m^{(I)}\rvert
We then iterate until self-consistency 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 Kohn-Sham equations during a self-consistent field (SCF) calculation; that is, the response is to be measured following the initial charge redistribution introduced by the perturbation but before the Kohn-Sham potential is updated. This approach is impractical to implement in codes that use a direct-minimization procedure of the total energy with respect to the density, Kohn-Sham orbitals, or density matrix.
Another approach to compute U and J is known as minimum tracking method [Linscott2018].
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:
\chi_{IJ} = \frac{dn^I}{dv_\text{ext}^J},
(\chi_0)_{IJ} = \left[\frac{dn}{dv_\text{KS}}\left(\frac{dv_\text{KS}}{dv_\text{ext}}\right)^{-1}\right]_{IJ}
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 rank-four tensors [Linscott2018]). In practice this is done by modifying the perturbation by including an additional term (spin-splitting):
V_{\text{ext}}^{p} = V_{\text{ext}} + \beta \sum_{m,m'}\lvert\varphi_{m'}^{(I\uparrow)}\rangle\langle\varphi_m^{(I\uparrow)}\rvert-\lvert\varphi_{m'}^{(I\downarrow)}\rangle\langle\varphi_m^{(I\downarrow)}\rvert
Any reasonable set of localised atomic-like 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 Kohn-Sham orbitals, so that the correlated subspaces are proper subspaces of the Kohn-Sham 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 on-site Coulomb repulsion [Miyake2008], and MLWFs are in common use as a minimal basis with which to construct tight-binding models from first-principles.
In ONETEP, a set of variationally-optimised nonorthogonal generalised Wannier functions (NGWFs) are generated as a by-product of total-energy minimisation. NGWFs exhibit some similar properties to MLWFs and other flavours of localised Wannier functions, and, for example, can be used to calculate finite-difference response properties in a similar way [O-Regan2012-2]. As they are conveniently available in ONETEP, we have made it possible to re-use the NGWFs from the end of a ground-state calculation as a set of Hubbard projectors with which to define the DFT+U correction. For this, it was necessary to develop a tensorially-consistent formulation of DFT+U in order to accommodate nonorthogonal projector functions [O-Regan2011]; 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 self-consistency cycle over the Hubbard
projectors, demonstrated in
Refs. [O-Regan2010], [O-Regan2011].
The density from one minimisation is re-used 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
self-consistency 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
post-processing tasks.
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
exchange-correlation 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 atom-centred 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 sub-shells per atom).
Hubbard U value
The value of the Hubbard U for this sub-shell, in electron-volts.
Hund's exchange J value
The value of the Hund’s exchange J for this sub-shell, in electron-volts. 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 Clementi-Raimondi effective charge, tabulated in Refs. [Clementi1963], [Clementi1967], and the choice of radial profile does matter [O-Regan2010].
In both cases, the projectors are effectively renormalised within an atom-centred 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 electron-volts. 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 linear-response theory [Cococcioni2005], [Kulik2006], or to break a spatial symmetry, such as in a mixed-valence system. In the example given, we are additionally penalising the occupancy on cerium 4f atomic orbitals.
The spin-splitting factor
The spin-splitting factor, in electron-volts, which is deducted from the \alpha factor for the spin-up channel and added to \alpha for the spin-down channel. In the example shown here we’re promoting spin-up magnetisation for iron atoms Fe1, and spin-down for Fe2. This can be very useful for appropriately breaking magnetic symmetries in antiferromagnetic solids or open-shell 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 carefully-defined atom-centred atomic population analysis, or breaking symmetries in spin or charge ordered systems.
| Keyword | Type | Default | Description |
|---|---|---|---|
HUBBARDSCF |
Task | — | Activate a projector-self-consistent
DFT+U calculation.
|
HUBBARDSCF_ON_THE_FLY |
Logical | False | Activate a non-variational on-the-fly form of
projector self-consistency in DFT+U
or cDFT, in which the projectors are updated
whenever the NGWFs are.
task:
HUBBARDSCF is then not needed. |
HUBBARD_CONV_WIN |
Integer | 2 |
The minimum number of Hubbard projector update
steps satisfying the incremental energy
tolerance
HUBBARD_ENERGY_TOL requiredfor convergence in task :
HUBBARDSCF. |
HUBBARD_ENERGY_TOL |
Physical | 1.0E-8 Ha |
The maximum incremental energy change between
Hubbard projector update steps allowed for
convergence in task :
HUBBARDSCF. |
HUBBARD_FUNCTIONAL |
Integer | 1 |
The form of DFT+U energy term used.
Contact developers if you need to try something
beyond the default.
|
HUBBARD_MAX_ITER |
Integer | 10 |
The maximum allowed number of Hubbard projector
update steps taken in a projector
self-consistent DFT+U or cDFT
calculation in task:
HUBBARDSCF |
HUBBARD_NGWF_SPIN_THRESHOLD |
Physical | 2.0E-5 Ha |
The incremental change in energy, in
total-energy minimisation, at which any
spin-splitting (Zeeman) type term in DFT+U is
switched off, and the minimisation history
reset. Useful for breaking open-shell,
antiferromagnetic, or charge-density
wave symmetries.
|
HUBBARD_PROJ_MIXING |
Real | 0.0 |
The fraction of previous Hubbard projector to
mix with new for projector self-consistent
DFT+U or cDFT in task :
HUBBARDSCF.Not found to be necessary.
|
HUBBARD_READ_PROJECTORS |
Logical | False |
Read Hubbard projectors from tightbox_hub_projs
file in restart calculations involving
DFT+U.
|
HUBBARD_TENSOR_CORR |
Integer | 1 |
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.
|
- 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
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:
- Total-energy 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.
- Conduction-band optimisation and Fermi’s Golden Rule spectra.
- Calculations of changes in electric polarisation.
- Time-dependent DFT.
- Electronic transmission calculations.
The extension of the DFT+U implementation to cluster Dynamical mean-field theory has also been implemented in ONETEP; for an example of its capabilities see Ref. [Weber2012].
| [O-Regan2012] | 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 |
| [Anisimov1991] | J. Z. V. I. Anisimov and O. K. Andersen, Phys. Rev. B 44, 943 (1991). https://doi.org/10.1103/PhysRevB.44.943 |
| [Anisimov1997] | V. I. Anisimov, F. Aryasetiawan, and A. I. Liechtenstein, J. Phys.: Condens. Matter 9, 767 (1997). https://iopscience.iop.org/article/10.1088/0953-8984/9/4/002 |
| [Dudarev1998] | S. L. Dudarev, Phys. Rev. B 57, 3 (1998). https://doi.org/10.1103/PhysRevB.57.1505 |
| [Aron_J.Cohen2012] | A. J. Cohen, P. Mor-Sanchez and W. Yang, Chem. Rev. 2012, 112, 289–320. https://doi.org/10.1021/cr200107z |
| [Cococcioni2005] | (1, 2, 3, 4) M. Cococcioni and S. de Gironcoli, Phys. Rev. B 71, 035105 (2005). https://doi.org/10.1103/PhysRevB.71.035105 |
| [Janak1978] | J. F. Janak, Phys. Rev. B 18, 12 (1978). https://doi.org/10.1103/PhysRevB.18.7165 |
| [Kulik2006] | (1, 2) 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 |
| [Himmetoglu2011] | (1, 2) B. Himmetoglu, R. M. Wentzcovitch, and M. Cococcioni, Phys. Rev. B,84, 115108 (2011). https://doi.org/10.1103/PhysRevB.84.115108 |
| [Clementi1963] | E. Clementi and D.L. Raimondi, J. Chem. Phys. 38, 2686 (1963). https://doi.org/10.1063/1.1733573 |
| [Clementi1967] | E. Clementi, D.L. Raimondi, and W.P. Reinhardt, J. Chem. Phys. 47, 1300 (1967). https://doi.org/10.1063/1.1712084 |
| [O-Regan2010] | (1, 2) 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 |
| [Weber2012] | 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 |
| [Marzari1997] | N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 (1997). https://doi.org/10.1103/PhysRevB.56.12847 |
| [Souza2001] | I. Souza, N. Marzari and D. Vanderbilt, Phys. Rev. B 65, 035109 (2001). https://doi.org/10.1103/PhysRevB.65.035109 |
| [Miyake2008] | T. Miyake and F. Aryasetiawan, Phys. Rev. B 77, 085122 (2008). https://doi.org/10.1103/PhysRevB.77.085122 |
| [O-Regan2012-2] | 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 |
| [O-Regan2011] | (1, 2) 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 |
| [Linscott2018] | (1, 2) 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 |