3.19 Hubbard corrected DFT (DFT+U)

So far this was just a brief sketch of the DFT+U approach in general. In the following we present the precise definition of DFT+U how it is implemented in FHI-aims. Without loss of generality we only show the equations with FLL as double-counting correction.

	$E_{\mathrm{U}}^{\mathrm{FLL}}\left[\left\{n_{Im{m}^{\prime }}^{\sigma }\right\}\right]$	$=E_{\mathrm{U}}^0\left[\left\{n_{Im{m}^{\prime }}^{\sigma }\right\}\right]-E_{\mathrm{dc}}\left[\left\{n_{Im{m}^{\prime }}^{\sigma }\right\}\right]$
		$={\displaystyle \frac{1}{2}}{\displaystyle \sum _{\sigma ,I}}{U}_{\mathrm{eff}}^{I}Tr\left[{𝐧}_I^{\sigma }\left(1-{𝐧}_I^{\sigma }\right)\right]$
		$={\displaystyle \frac{1}{2}}{\displaystyle \sum _{\sigma ,I}}{U}_{\mathrm{eff}}^I\left[Tr\left({𝐧}_I^{\sigma }\right)-Tr\left({𝐧}_I^{\sigma }{𝐧}_I^{\sigma }\right)\right].$		(3.61)

These functional is known as the spherically averaged form of DFT+U. It was first proposed by Dudarev et al.[79] and it is also rotational invariant. In this formulation, the effective on-site interactions enter via their spherical atomic averages. This is justified by the fact, that localized states still have atomic character and hence, spherical symmetry. In fact, for most materials this definition gives good results.

It should be pointed out, that $U_{\mathrm{eff}}$ can be seen as an effective value of the coulomb interaction that also includes exchange corrections. This parameter has to be specified by hand, so far, no possibility is implemented to calculate this parameter self-consistently. Common to all approaches is that all the calculated results sensitively depend on the applied $U_{\mathrm{eff}}$ value. This value not only depends on the atom for which DFT+U is applied. It also depends on the surroundings of the atom, the lattice parameters and physical conditions. Furthermore, it also depends on the localized basis set of the underlying quantum DFT code. This limits the comparability of different values in a strong way. In general, for each DFT+U implementation and system, one should recalculate $U_{\mathrm{eff}}$.

The most important quantity in equation 3.19.1 is the so called DFT+U occupation matrix $𝐧$. This matrix simply tells us how many electrons are in a certain shell on a certain atom. The problem here is the inability to break down the total charge density into atom specific contributions. Or in other words, there is no proper operator for counting the number of electrons on an atom. Hence, the choice of the occupation matrix will affect the outcome of a calculation. Within FHI-aims we offer two specific choices: the on-site representation of the occupation matrix

$n_{Im{m}^{\prime }}^{\sigma }(\mathrm{on}-\mathrm{site})=D_{Im,I{m}^{\prime }}^{\sigma }$

(3.62)

and

$n_{Im{m}^{\prime }}^{\sigma }(\mathrm{dual})={\displaystyle \frac{1}{2}}{\displaystyle \sum _{Jk}}\left[D_{Im,Jk}^{\sigma }{S}_{Jk,I{m}^{\prime }}+S_{Im,Jk}{D}_{Jk,I{m}^{\prime }}^{\sigma }\right].$

(3.63)

The latter is known as the dual representation [126]. Within the dual representation the occupation numbers are calculated in a similar way as in the Mulliken analysis. The main difference between both is that the on-site version only accounts for overlaps within a specific sub shell on an certain atom. The dual representation also accounts for the overlap with the surrounding atoms. It is emphasized that all general aspects of DFT+U are met by all matrix representations. Furthermore, more detailed studies regarding the performance of the occupation matrix for various transition metal oxides revealed that in principle there is no definition which is clearly the best [290]. Unfortunately, we only offer forces for the on-site representation. The on-site version is also the default occupation matrix in FHI-aims and we strongly recommend to use it.

By now, one might have noticed that DFT+U is by far not a black box method and it gets even worse if one considers in detail how the occupation matrix is constructed. In general, each occupation matrix can be expressed in terms of a local projector operator, ${\widehat{P}}_{Im{m}^{\prime }}^{\sigma }$. The ($m$,$m^{\prime }$)-th element of a occupation matrix at site $I$ is then given by

$n_{Im{m}^{\prime }}^{\sigma }={\displaystyle \sum _{\gamma }}{f}_{\gamma }⟨{\mathrm{\Psi }}_{\gamma }^{\sigma }|{\widehat{P}}_{Im{m}^{\prime }}^{\sigma }|{\mathrm{\Psi }}_{\gamma }^{\sigma }⟩.$

(3.64)

For example for the on-site projection operator this would lead to

${\widehat{P}}_{Im{m}^{\prime }}^{\sigma }(on\text{-}site)=|{\tilde{\varphi }}_{I{m}^{\prime }}^{\sigma }⟩⟨{\tilde{\varphi }}_{Im}^{\sigma }|.$

(3.65)

Here, ${\tilde{\varphi }}_{Im}$ denotes the dual basis functions which are defined in terms of the inverse overlap matrix ${𝐒}^{-1}$,

$|{\tilde{\varphi }}_m^{\sigma }⟩={\displaystyle \sum _{I^{\prime }{m}^{\prime }}}{S}_{Im,I^{\prime }{m}^{\prime }}^{-1}|{\tilde{\varphi }}_{m^{\prime }}^{\sigma }⟩.$

(3.66)

The question now is which basis functions should be used in the projection? As default we are using the atomic type basis functions of the minimal basis set in FHI-aims. Here we automatically assume, as they are atomic like basis functions, that they will contribute most to the localized states. However, in general it is not known if other basis functions should also be included in the DFT+U projection e.g. tier1 $3d$ if one deals with first row transition metals. Usually one can notice that by a strange behavior of the occupation matrix during the scf-cycle (occupation numbers drop to zero as the electrons occupy other basis functions). For that purpose we offer to include also other basis functions in the description of DFT+U (see plus_u_use_hydors). We also like to highlight the corresponding paper related to our implementation where we address fundamental issues of DFT+U in a LCAO electronic structure code. However, do not panic, for most of the systems the default settings should be sufficient enough.

So far, we presented DFT+U in quite some detail. However, we just wanted to highlight that DFT+U is far from being a black box method. However, the handling of a actual DFT+U calculation in FHI-aims is quite easy. One just have to specify the double-counting correction first via the plus_u_petukhov_mixing. Afterwards one can specify the U value and the angular momentum shell to which DFT+U should be applied for each species. Of course one can specify different U values for different species in a simulation. Only for hard cases where convergence can not be reached easily, it is quite useful to checkout the other keywords. Some of them can be quite useful such as the plus_u_matrix_control. Here, one first converges the density with help of a fixed occupation matrix which can be edited by hand. Afterwards one can use the restart information to calculate everything self-consistently. This can be quite useful as it turns out that DFT+U is quite sensitive to the initial guess of a calculation. Furthermore, it is quite useful also to start from a LDA or GGA ground state density.

There are two common schemes for dealing with the double counting problem in DFT+U: The AMF method assumes that the effect of the DFT+U term on the actual occupations remains small, so that the occupations can be assumed to be equal within each shell for the purpose of the double counting correction. The FLL method, on the other hand, assumes a maximal effect of the DFT+U term on the occupation numbers, handling double counting correctly in the case that all orbitals with in the shell are either fully occupied or empty. The self consistent mixing of both limits improves the handling of the intermediate range (see Ref. [248]).

This is extremely useful because one can simply edit the file and manipulate the matrix according to some specific spin configuration. Consider to use it with restart options.