What is the appropriate way of determining a value for the Hubbard-like potential U for LDA+U / GGA+U calculations?

Question

It is known that neither LDA or GGA correctly account for the strong on-site Coulomb interaction of localized electrons[1]. A cost-effective way of correcting this is by using a Hubbard-like term $U$. In literature, I have often found that $U$ is introduced as an empirical parameter, taken from previous computational studies or experiments. Alternatively, I have seen studies where various values for $U$ are used (typically between 2 and 5 eV). Neither of these methods seem too convinving. Is there a standard and reliable manner of determining the (best) value for $U$?

Answer 1

This is a bit late, but I would say the short answer to your direct question is: technically no, there is no "standard" reliable way, since there are several approaches to determining U self-consistently from first principles. By self-consistently I mean a first-principles U is calculated, which when applied changes the electronic structure, so a new U is calculated and applied, and so on, until the value of U has converged. However, these methods are each fairly reliable and often give similar results. This is an area I've worked in a bit so I'll summarize a few methods here, some ongoing work in this area, some resources for further reading and practical application, and some recommended practices. This is in the context of my experience in solid-state periodic DFT but in principle it should apply to any DFT calculation where these methods are possible and appropriate.

1. Hubbard U and Determination of U from First-Principles

a. Basic Concepts

The basic idea of DFT+U is that we are replacing the treatment of some of the charge density with a Hubbard correction term, where otherwise it would be treated by the chosen DFT exchange-correlation (XC) functional. This correction, in practice, is usually a screened Hartree-Fock-like term as opposed to being more reminiscent of the actual Hubbard model. A double-counting term is subtracted to remove the estimated contribution from the original XC functional. However, this term is not uniquely defined, so the choice of double-counting correction will affect the results. Most calculations use a fully-localized limit, most appropriate for strongly localized states.

DFT+U has the effect of introducing discontinuity in the total energy as a function of electrons added to the system--the rationale behind why it can improve the prediction of band gaps (if this is unclear, read on the band gap problem in DFT). In fact, the value of U can be defined in terms of this unphysical curvature present in LDA and GGA. A more intuitive interpretation of DFT+U can be illustrated by the form of the rotationally-invariant "simplified" correction (reducing explicit U and J terms to an effective $U_{\text{eff}}=U-J$) developed by Dudarev et al. and Cococcioni et al., where the Hubbard energy is minimized when orbitals in the Hubbard manifold are either completely full or empty, i.e. no partial occupations due to hybridization. This means there is a tendency to decrease delocalization in the calculation (compensating for how GGA, for example, tends to over-delocalize).

A more recent development is DFT+U+V, which includes inter-site interaction V and can improve the description of covalently-bonded materials.

I recommend reading this excellent review article by Himmetoglu et al. for an overview of DFT+U.

b. Constrained Random Phase Approximation

The constrained random phase approximation (cRPA), developed by Aryasetiawan et al., is used to find a frequency-dependent U that can be used in DFT+DMFT, but in the static limit it can be used in DFT+U. Essentially, once you have chosen how to define the localized and delocalized states for your system, the Coulomb interaction between the localized states is calculated while including screening effects from the delocalized states. The "constrained" comes from the fact that only the Hartree term is used to calculate the dielectric function used for screening, for simplicity (as opposed to both Hartree and exchange-correlation). I personally have not used cRPA to calculate values of U, so check the reference paper and its citations for further reading on the topic.

c. Linear Response

The linear response method by Cococcioni et al. defines U such that when applied the unphysical curvature in the total energy vs. the number of electrons present in the system is eliminated. In exact DFT this is a piecewise-continuous linear function, while in approximate LDA and GGA DFT it is a smooth curved function. Constrained DFT was an approach to correct this by varying the Hubbard orbital occupations, and determining the corresponding change in energy. Linear response approaches this in a more convenient way for most DFT codes, by applying a varying perturbative potential and then measuring the resulting change in occupation. Once this is done for several perturbations, U can be calculated (see resources section for how to do this). This approach typically requires U to be calculated in a supercell, to prevent the Hubbard states from being affected by periodic images of the perturbation.

d. Density Functional Perturbation Theory

A very nice recent development in the calculation of U is from Timrov, Mazari and Cococcioni. They reformulate the linear response method from a single perturbation in a supercell to a sum of perturbations in the primitive cell. I'm not as familiar with the theory here, but it's implemented in the new hp.x code included with Quantum Espresso. Very interesting.

d. ACBN0

ACBN0 (named for the authors) takes inspiration from previous work by Mosey and Carter, and explicitly calculates U from a Hartree-Fock-like interaction between the Hubbard orbitals of interest. Some screening-like reduction of the interaction is introduced by renormalizing the occupations of the Kohn-Sham orbitals according to their projectability on the Hubbard basis--so less-localized states should have a drastically reduced magnitude of U. This method also allows calculation of many site-dependent U values from a single scf calculation. In theory this can be incorporated into the self-consistency loop of the DFT calculation, but current implementations are a post-processing step. It can be used in the PAOFLOW and AFLOW$\pi$ codes. It's also been recently demonstrated with DFT+U+V.

2. Resources

a. References for Implementations of Hubbard U Corrections

• Rotationally-invariant implementation with explicit U and J
• Simplified rotationally-invariant scheme ($U_{\text{eff}}$) (also this)
• Simplified scheme re-introducing explicit J term

b. Future Possibilities

A limitation of the simplified implementations of DFT+U is that the Hubbard orbitals are treated in a way that assumes some spherical symmetry (i.e. very close to true atomic orbitals). In compounds that have significant crystal field splitting this is not technically appropriate. Some early work used different U values for $t_{2g}$ and $e_{g}$ electrons in perovskite oxides, for example. The exchange term J is also treated in an average way, which worsens treatment of materials where the localization depends on Hund's rule magnetism. It could be very interesting to determine U, J and V for specific subsets of the Hubbard orbitals. This may provide better treatment of orbitally-ordered materials, or materials where some of the d electrons form an itinerant band (and should have a much smaller U value than the localized states, and may also participate in screening).

c. Practical Information/Tutorials

Prof. Heather Kulik at MIT has some nice tutorials and slides on DFT+U:

• Intro to DFT+U
• Calculating Hubbard U, U for multiple sites, U for periodic systems, calculating U and J
• Troubleshooting common problems, the right U for any element
• Including Hubbard U variations with DFT+U(R)

3. Recommended Practices

1. Values of U are in general non-transferable. It is poor practice to take a literature value of U when the work does not use the same functional, DFT+U implementation, material (especially drastically different chemical environments), Hubbard basis, or even pseudopotential that you are using. This is especially true for values of U chosen empirically for a different property than the one you're studying.
2. Calculate U self-consistently whenever possible. It's becoming easier and easier to do this.
3. Think carefully about which atoms and states you want to apply U to. It's sometimes useful to apply U to oxygen p states in addition to metal d states in oxides, for example.
4. DFT+U often introduces local minima in the energy, related to the occupations of your Hubbard atom states. This seems especially pronounced in spin-polarized calculations. Manually setting the starting occupation matrix to what you are trying to study can help converge to the appropriate energy minimum. This is useful when trying to study different possible spin states on a metal cation, for example. It can also help speed up convergence, in my experience.
5. If you are choosing U empirically, think about whether you are getting the "right" answer for the wrong reason. Band gaps can be correct at the expense of a realistic picture of electronic structure, not to mention other quantities like lattice parameters, formation enthalpy, etc. Compare your calculations with higher levels of theory in the literature when possible.
6. Current DFT+U implementations are not, to my knowledge, variational with respect to the total energy. So the optimal value of U will not typically result in the lowest calculated total energy. To make comparisons between different calculations, you need to be using the same values of U in both calculations--unless you specifically account for the change in the energy surface with respect to U (as in DFT+U(R) or DFT+U(V)).

Answer 2

You can calculate U in an ab initio way via linear response theory, for instance. See an example here. However, there is no guarantee that a linear response U value is empirically ideal. That is part of the reason why, in practice, many people invoke some empirical U term and hope that it holds reasonably well for their system (especially compared to U = 0 eV). In the case of using a range of U values, the goal is usually to put some sort of lower and upper bound on the computed property of interest, since the ideal U value is not known.

Answer 3

The answer is based on Chapter 1 - The DFT+U: Approaches, Accuracy, and Applications from the book Density Functional Calculations: Recent Progresses of Theory and Application

For practical implementation of DFT+U, the strength of the on-site interactions is described by a couple of parameters: the on-site Coulomb term U and the site exchange term J. These parameters “U and J” can be extracted from ab initio calculations, but usually are obtained semi-empirically. The implementation of the DFT+U requires a clear understanding of the approximations it is based on and a precise evaluation of the conditions under which it can be expected to provide accurate quantitative predictions

Different procedures have been addressed to determine the Hubbard U from first principles. In these procedures, the U parameter can generally be calculated using a self-consistent and basis set in an independent way. For each atom, the U value is found to be dependent on the material specific parameters, including its position in the lattice and the structural and magnetic properties of the crystal, and also dependent on the localized basis set employed to describe the on-site occupation in the Hubbard functional. Most programs these days use the method presented by Cococcioni et al., where the values of U can be determined through a linear response method, in which the response of the occupation of localized states to a small perturbation of the local potential is calculated. The U is self-consistently determined, which is fully consistent with the definition of the DFT+U Hamiltonian, making this approach for the potential calculations fully ab initio.

The U value is not known and practically is often tuned semi-empirically to make a good agreement with experimental or higher level computational results. However, the semi-empirical way of evaluating the U parameter fails to capture its dependence on the volume, structure, or the magnetic phase of the crystal, and also does not permit the capturing of changes in the on-site electronic interaction under changing physical conditions, such as chemical reactions. Despite the limitations of choosing the U value semiempirically for systems, where variations of on-site electronic interactions are present, it is found to be the most common practice used in literature, where the value of U is usually compared to the experimental bandgap.The semi-empirical tuning is found to be the most common practice employed by researchers due to the significant computational cost of ab initio calculations that U can have, and also, the computed U is not necessarily being better than the empirical ones.

A review on different values of U (ranging from -2 eV to 10 eV ) used for different materials are also available in the chapter