The DFT+U: Approaches, Accuracy, and Applications

Abstract

This chapter introduces the Hubbard model and its applicability as a corrective tool for accurate modeling of the electronic properties of various classes of systems. The attainment of a correct description of electronic structure is critical for predicting further electronic-related properties, including intermolecular interactions and formation energies. The chapter begins with an introduction to the formulation of density functional theory (DFT) functionals, while addressing the origin of bandgap problem with correlated materials. Then, the corrective approaches proposed to solve the DFT bandgap problem are reviewed, while comparing them in terms of accuracy and computational cost. The Hubbard model will then offer a simple approach to correctly describe the behavior of highly correlated materials, known as the Mott insulators. Based on Hubbard model, DFT+U scheme is built, which is computationally convenient for accurate calculations of electronic structures. Later in this chapter, the computational and semiempirical methods of optimizing the value of the Coulomb interaction potential (U) are discussed, while evaluating the conditions under which it can be most predictive. The chapter focuses on highlighting the use of U to correct the description of the physical properties, by reviewing the results of case studies presented in literature for various classes of materials.

1. Introduction

Density functional theory (DFT) is one of the most convenient computational tools for the prediction of the properties of different classes of materials [1, 2]. Although its accuracy is acceptable as long as structural and cohesive properties are concerned, it dramatically fails in the prediction of electronic and other related properties of semiconductors up to a factor of two [3]. However, reaching a correct description of electronic structure is critical for predicting further electronic-related properties, including intermolecular interactions and formation energies. In order to solve this problem, computationally heavier jobs must be employed, using either larger basis sets or hybrid functionals, which include the solution of the exact Hartree-Fock (HF) equations, in order to reach relatively higher accuracies [4]. Nevertheless, in some cases, even solving exact HF equations can fail in correctly predicting the bandgap for a certain class of semiconductors that possess strong correlations between electrons, such as Mott insulators [5, 6]. Consistent research efforts have been employed in order to formulate more accurate functionals, by using corrective approaches or alternatives to the density functionals. The applicability of these alternatives and corrections has large dependence on the type of the system studied, its size and complexity, and the computational cost required. One of the corrective approaches employed to relieve the DFT electronic bandgap problem is the DFT+U correction method, which is the focus of this chapter. Compared to the alternative approaches, such as the hybrid functionals and the post-Hartree-Fock methods, DFT+U correction has proved to be as reliable as the other methods, but with a critical advantage of considerably lower computational cost. By successfully correcting the electronic structure of the studied system using the U correction, further accurate predictions of intermolecular interactions and formation energies can be reached [6]. In addition, the U correction can further enhance the description of physical properties, other than the electronic structure, including magnetic and structural properties of correlated systems, the electron transfer energetics, and chemical reactions. However, one of the drawbacks of the Hubbard method is that it fails in predicting the properties of systems with more delocalized electrons, such as metals. The relative success of the DFT+U method is related to its straightforward approach to account for the underestimated electronic interactions by simply adding a semiempirically tuned numerical parameter “U” [5]. This interaction parameter can be easily controlled, making the DFT+U method a tool to give a qualitative assessment of the influence of the electronic correlations on the physical properties of a system.

One of the mostly implemented methods in the DFT+U realm is the LDA+U method. It is widely used due to its simple implementation on the existing LDA codes, which makes it only slightly computationally heavier than the standard DFT computations [6]. In this chapter, we discuss the fundamental formulation of the LDA+U method and examine its applicability for practical implementations for different classes of materials, where DFT is usually found to be impractical. Popular cases of DFT shortage are discussed including materials with strong correlations, defective solid-state materials, and organometallics, while reviewing literature case studies that studied these classes of materials with DFT+U calculations. The methodology of optimizing the U correction is inspected, where it can be either formulated from first principles or achieved empirically by tuning the U value, while seeking an agreement with experimental results of the system’s physical properties. In this chapter, we also present a review of the practical implementation of U, while assessing its corrective influence on improving the description of a variety of physical properties related to certain classes of materials. In addition, the effect of the calculation parameters on the chosen U value is discussed, including the choice of the localized basis set and the type of DFT functional employed.

2. Theoretical formulation

2.1. Standard DFT problem

Using exact HF or DFT solutions, the aim is always to reach, as close as possible, the exact description of the total energy of the system. Unluckily, reaching this exact energy description is impossible and approximations have to be employed. In DFT, electronic interaction energies are simply described as the sum of classical Columbic repulsion between electronic densities in a mean field kind of way (Hartree term) and an additive term that is supposed to encompass all the correlations and spin interactions [1]. This additive term, namely the exchange and correlation (xc), is founded on approximations that have the responsibility to recover the exact energy description of the system. This approximated xc functional is a function of the electronic charge density of the system, and the accuracy of a DFT calculation is strongly dependent on the descriptive ability of this functional of the energy of the system [2]. It is generally difficult to model the dependence of the xc functional on electronic charge density, and thus, it can inadequately represent the many-body features of the N-electron ground state. For this reason, systems with physical properties that are controlled by many body electronic interactions (correlated systems) are poorly described by DFT calculations. For these systems, incorrect description of the electronic structure induces the so-called “bandgap problem,” which in turn, imposes difficulties in utilizing DFT to predict accurate intermolecular interactions, formation energies, and transition states [7].

The problem of DFT to describe correlated systems can be attributed to the tendency of xc functionals to over-delocalize valence electrons and to over-stabilize metallic ground states [5, 6]. That is why DFT fails significantly in predicting the properties of systems whose ground state is characterized by a more pronounced localization of electrons. The reason behind this delocalization is rooted to the inability of the approximated xc to completely cancel out the electronic self-interaction contained in the Hartree term; thus, a remaining “fragment” of the same electron is still there that can induce added self-interaction, consequently inducing an excessive delocalization of the wave functions [5]. For this reason, hybrid functionals were formulated to include a linear combination of a number of xc explicit density and HF exact exchange functionals, that is self-interaction free, by eliminating the extra self-interaction of electrons through the explicit introduction of a Fock exchange term. However, this method is computationally expensive and is not usually practical when larger, more complex systems are studied. Nonetheless, HF method, which describes the electronic structure with variationally optimized single determinant, cannot describe the physics of strongly correlated materials such as the Mott insulators. In order to describe the behavior of these systems, full account of the multideterminant nature of the N-electron wave function and of the many-body terms of the electronic interactions is needed [6]. Therefore, it is predicted that applying DFT calculations using approximate xc functionals, such as LDA or GGA, will poorly describe the physical properties of strongly correlated systems.

2.2. Mott insulators and the Hubbard model

According to the conventional band theories, strongly correlated materials are predicted to be conductive, while they show insulating behavior when experimentally measured. This serious flaw of the band theory was pointed out by Sir Nevil Mott, who emphasized that interelectron forces cannot be neglected, which lead to the existence of the bandgap in these falsely predicted conductors (Mott insulators) [8]. In these “metal-insulators,” the bandgap exists between bands of like character i.e., between suborbitals of the same orbitals, such as 3d character, which originates from crystal field splitting or Hund’s rule. The insulating character of the ground state stems from the strong Coulomb repulsion between electrons that forces them to localize in atomic-like orbitals (Mott localization). This Coulomb potential, responsible for localization, is described by the term “U,” and when electrons are strongly localized, they cannot move freely between atoms and rather jump from one atom to another by a “hopping” mechanism between neighbor atoms, with an amplitude t that is proportional to the dispersion (the bandwidth) of the valence electronic states. The formation of an energy gap can be settled as the competition between the Coulomb potential U between 3d electrons and the transfer integral t of the tight-binding approximation of 3d electrons between neighboring atoms. Therefore, the bandgap can be described by the U, t and an extra z term that denotes the number of nearest neighbor atoms as [6]:

$${\mathrm{E}}_{\mathrm{gap}}=\mathrm{U}-\text{2zt}$$E1

Since the problem is rooted down to the band model of the systems, alternative models have been formulated to describe the correlated systems. One of the simplest models is the “Hubbard” model [9]. The Hubbard model is able to include the so-called “on-site repulsion,” which stems from the Coulomb repulsion between electrons at the same atomic orbitals, and can therefore explain the transition between the conducting and insulating behavior of these systems. Based on this model, new Hamiltonian can be formulated with an additive Hubbard term that explicitly describes electronic interactions. The additive Hubbard Hamiltonian can be written in its simplest form as follows [6]:

$$H_{\mathit{Hub}}=t\sum _{⟨i,j⟩,\sigma }\left(c_{i,\sigma }^{†}\right.c_{j,\sigma }+\left.h.c.\right)+U\sum _i{n}_{i,\uparrow }{n}_{i,\downarrow }$$E2

As predicted, the Hubbard Hamiltonian should be dependent on the two terms t and U, with ⟨i.j⟩ denoting nearest-neighbor atomic sites and c_i^†, c_j, and n_i are electronic creation, annihilation, and number operators for electrons of spin on site i, respectively. The hopping amplitude t is proportional to the bandwidth (dispersion) of the valence electrons, while the on-site Coulomb repulsion term U is proportional to the product of the occupation numbers of atomic states on the same site [6]. The system’s insulating character develops when electrons do not have sufficient energy to overcome the repulsion potential of other electrons on neighbor sites, i.e., when t « U. The ability of the DFT scheme to predict electronic properties is fairly accurate when t » U, while for large U values, DFT significantly fails the HF method, which describes the electronic ground state with a variationally optimized single determinant, that cannot capture the physics of Mott insulators.

2.3. DFT+U

Inspired by the Hubbard model, DFT+U method is formulated to improve the description of the ground state of correlated systems. The main advantage of the DFT+U method is that it is within the realm of DFT, thus does not require significant effort to be implemented in the existing DFT codes and its computational cost is only slightly higher than that of normal DFT computations. This “U” correction can be added to the local and semilocal density functionals offering LDA+U and GGA+U computational operations. The basic role of the U correction is to treat the strong on-site Coulomb interaction of localized electrons with an additional Hubbard-like term. The Hubbard Hamiltonian describes the strongly correlated electronic states (d and f orbitals), while treating the rest of the valence electrons by the normal DFT approximations. For practical implementation of DFT+U in computational chemistry, 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 semiempirically. 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 [5, 6].

The LDA+U method is widely implemented to correct the approximate DFT xc functional. The LDA+U works in the same way as the standard LDA method to describe the valence electrons, and only for the strongly correlated electronic states (the d and f orbitals), the Hubbard model is implemented for a more accurate modeling. Therefore, the total energy of the system (E_LDA+U) is typically the summation of the standard LDA energy functional (E_Hub) for all the states and the energy of the Hubbard functional that describes the correlated states. Because of the additive Hubbard term, there will be a double counting error for the correlated states; therefore, a “double-counting” term (E_dc) must be deducted from the LDA’s total energy that describes the electronic interactions in a mean field kind of way [5].

$${\mathrm{E}}_{\mathrm{LDA}}+\mathrm{U}\left[\mathrm{\rho }\left(\mathrm{r}\right)\right]={\mathrm{E}}_{\mathrm{LDA}}\left[\mathrm{\rho }\left(\mathrm{r}\right)\right]+{\mathrm{E}}_{\mathrm{Hub}}\left[\left\{n_{\mathit{mm}}^{\mathit{I\sigma }}\right\}\right]-{\mathrm{E}}_{\mathrm{dc}}\left[{\mathrm{n}}^{\mathrm{I\sigma }}\right]$$E3

Therefore, it can be understood that the LDA+U is more like a substitution of the mean-field electronic interaction contained in the approximate xc functional. Nonetheless, the E_dc term is not uniquely defined for each system and various formulations can be applied to different systems. The most dominant of these formulations is the FLL formulation [10, 11, 12]. It is based on the implementation of fully localized limit (FLL) on systems with more localized electrons on atomic orbitals. The reason for this formulation popularity is due to its ability to expand the width of the Kohn Sham (KS) orbitals and to effectively capture Mott localization. Based on this formulation, the LDA+U can be written as:

$$E_{\mathit{LDA}+U}\left[\rho \left(r\right)\right]+\sum _l\left[\frac{U^l}{2}\sum _{m,\sigma \ne m^{\prime },{\sigma }^{\prime }}{n}_m^{\mathit{l\sigma }}{n}_{m^{\prime }}^{{\mathit{l\sigma }}^{\prime }}-\frac{U^l}{2}{n}^l\left(n^l-1\right)\right]$$E4

where $n_m^{\mathit{l\sigma }}$ m are the localized orbitals occupation numbers identified by the atomic site index I, state index m, and by the spin $\sigma$. In Eq. (4), the right-hand side second and third terms are the Hubbard and double-counting terms, specified in Eq. (3). The dependency on the occupation number is expected as the Hubbard correction is only applied to the states that are most disturbed by correlation effects. The occupation number is calculated as the projection of occupied KS orbitals on the states of a localized basis set:

$$n_{m,m^{\prime }}^{\mathit{I}\mathit{\sigma }}={\sum }_{k,v}{f}_{\mathit{kv}}^{\mathrm{ }\sigma }⟨{\mathrm{\Psi }}_{\mathit{kv}}^{\sigma }|{\varphi }_{m^{\prime }}^I⟩⟨{\varphi }_m^I|{\mathrm{\Psi }}_{\mathit{kv}}^{\sigma }⟩$$E5

where the coefficients $f_{\mathit{kv}}^{\mathrm{ }\sigma }$represent the occupations of KS states (labeled by k-point, band, and spin indices), determined by the Fermi-Dirac distribution of the corresponding single-particle energy eigen values. According to this formulation, the fractional occupations of localized orbitals is reduced, while assisting the Mott localization of electrons on particular atomic states [5].

Although the above approach described in Eq. (4) is able to capture Mott localization, it is not invariant under rotation of the atomic orbital basis set employed to define the occupation number of n in Eq. (5). This variation makes the calculations performed unfavorably dependent on the unitary transformation of the chosen localized basis set. Therefore, “rotationally invariant formulation” is introduced, which is unitary-transformation invariant of LDA+U [12]. In this formulation, the electronic interactions are fully orbital dependent, and thus considered to be the most complete formulation of the LDA+U. However, a simpler formulation that preserves rotational invariance, which is theoretically based on the full rotationally invariant formulation, had proved to work as effectively as the full formulation for most materials [11]. Based on the simplified LDA+U form, it has been customary to utilize, instead of the interaction parameter U, an effective U parameter: U_eff = U−J, where the “J” parameter is known as the exchange interaction term that accounts for Hund’s rule coupling. The U_eff is generally preferred because the J parameter is proven to be crucial to describe the electronic structure of certain classes of materials, typically those subject to strong spin-orbit coupling.

3. Practical implementations of the Hubbard correction

DFT+U is applicable for all open shell orbitals, such as d and f orbitals for transition metal elements with localized orbitals existing in extended states, as in the case of many strongly correlated materials and perovskites, where localized 3d or 4f orbitals are embedded in elongated s-p states. A complicated many-electron problem is made of electrons living in these localized orbitals, where they experience strong correlations among each other and with a subtle coupling with the extended states. Isolating a few degrees of freedom relevant to the correlation is the idea in the Hubbard model, where screened or renormalized Coulomb interaction (U) is kept among the localized orbitals’ electrons [13]. In other word, the localized orbitals in the bandgap, which are present as localized states (d- and f-states), are too close to the Fermi energy. From that aspect, the U value should be used to push theses states away from the Fermi level, such as that provided by the GGA+U theory, which adds to the Hamiltonian a term that increases the total energy preventing the unwanted delocalization of the d- or f-electrons, when two d- or f-electrons are located on the same cation [14]. It is worth mentioning that using too large values of U will over-localize the states and lead to an unphysical flattening of the appropriate bands, which unlike fitting to many other properties, will make the fit worse. Also, the increase in the U value can cause an overestimation of the lattice constants as well as a wrong estimation of the ground state energy due to the electronic interaction error. Therefore, applying Hubbard correction to solve the bandgap problem is necessary for predicting the properties of transition metal oxides. Figure 1 shows the effect of U potential on correcting the failure of DFT to predict correct bandgaps for strongly correlated materials. Note the underestimation of the bandgap in case of MnO and the incorrect prediction of the metallic behavior of FeO [15].

Figure 1.

Comparison of theoretical DOS calculated by LDA and LDA+U for (a) MnO and (b) FeO. Adapted with permission from [15].

3.1. Optimizing the U value

From the case studies and examples presented within this chapter, one can intuitively conclude that corrective functional LDA+U is particularly dependent on the numerical value of the effective potential U_eff, which is generally referred to in literature as “U” for simplicity. However, the U value is not known and practically is often tuned semiempirically to make a good agreement with experimental or higher level computational results. However, the semiempirical 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. In order to get full advantage of this method, different procedures have been addressed to determine the Hubbard U from first principles [13]. In these procedures, the U parameter can generally be calculated using a self-consistent and basis set in an independent way. These different ab initio approaches for calculating U have been applied to different material systems, where the U value is calculated for individual atoms. 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. Therefore, the value of effective interactions should be re-computed for each type of material and each type of LDA+U implementation (e.g., based on augmented plane waves, Gaussian functions, etc.). Most programs these days use the method presented by Cococcioni et al. [16], where the values of U can be determined through a linear response method [17], 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 value of U implemented by Cococcioni et al. is U_eff = U−J, where J is indirectly assumed to be zero in order to obtain a simplified expression [17]. Nonetheless, J can add some additional flexibility to the DFT+U calculations, but it may yield surprising results including reversing the trends previously obtained in the implemented DFT+U calculations [18].

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. This semiempirical trend in practical implementation of U is present because of the significant computational cost of ab initio calculation of U, and in the cases of studying static physical properties, the results of computed U are not necessarily found to be better than the empirical ones. Within this practice, however, caution should be taken while pursuing the semiempirical method [19]. If it will be possible to describe all the relevant aspects of a system, except the bandgap, with a reasonable U, one might then look into using a scissor operator or rigid shift to the bandgap [20, 21]. However, in particular cases, where calculations aim at understanding catalysis, it is natural to choose U to fit the energy of the oxidation-reduction, as catalysis is controlled by energy differences [14]. Conversely, one of the possible solutions is to venture into a negative value for the Hubbard U parameter, there is no obvious physical rationale for that yet, but the results may match with experiment for both the magnetic moment and structural properties, as illustrated later in this chapter.

To elaborate the numerical U tuning procedure, three quick examples are presented below that can show the correlation between the value of U and the predicted physical properties:

• The compilation of the correlated nature of cobalt 3d electrons in the theoretical studies of Co₃O₄ gives a good picture of the significant difference in the U value with the difference in most of the calculated properties. The variety of U values have been used ranging from 2 to 6 for the properties including bandgap [22], oxidation energy [23], and structural parameters [24], which affect the choice of the value of U for each of these properties uniquely. The calculated bandgap at the generalized gradient approximation (GGA)+U agrees well with the experimental value of 1.6 eV. On the other hand, the calculated value using the PBE0 hybrid functional (3.42 eV) is highly overestimated, due to neglecting the screening problem of the Hartree-Fock approximation [25].

• From the study made by Lu and Liu [26] on cerium compounds presented some characteristics for U values for Ce atoms in different configurations as isolated atoms and ions. They illustrated that the ion charge (Ce atoms, Ce in Ce₃H_xO₇ clusters, or CeO₂) does not significantly affect the value of U and that when ions are isolated, the values are much larger (close to 15 for Ce^2.5+ and 18 eV for Ce^3.5+).

• Within the study of BiMnO₃ that has strongly distorted perovskite structure with the GGA+U method, calculations show that distortions of the MnO₆ octahedral, which is considered the main unit in the crystal structure, are very sensitive to the value of the Coulomb repulsion U. The study showed that large U value decreases the 3d-2p hybridization, and therefore decreases the bonding effects, which in turn distortion increases the short Mn-O distances, and thus overly expands the MnO₆ octahedral [27].

3.2. Variation of U with calculation methods and parameters

The parameters assigned for DFT calculations can significantly affect the choice of the optimum U value. These parameters include pseudopotentials, basis sets, the cutoff energy, and k-point sampling. As pseudopotentials are used to reduce computational time by replacing the full electron system in the Columbic potential by a system only taking explicitly into account the “valence” electrons [28], the pseudopotential will strongly affect the U value. Thus, calculations have to be converged very well with respect to the cutoff energy and k-point sampling, while taking into account that the symmetry used in DFT+U calculations, because adding the U parameter often lowers the crystallographic symmetry, thereby the number of k-points needs to be increased.

Not only is the U value affected by the parameters applied, but it is also strongly dependent on the DFT method used. In a published review [29], a comparison of different calculated U values using different approaches was highlighted for several transition metal oxides. It was reported that with small U values, the electrons were still not localized, and that the U value depends on the used exchange-correlation functionals (LDA or GGA), the pseudopotential, the fitted experimental properties, and projection operators [30]. In the computational study of strongly correlated systems, it can be usually found in literature that researchers refer to the utilization of (DFT+U) method in their calculations, which may include generalized gradient approximation (GGA+U) [8], local density approximation (LDA+U) [31], or both [32]. To be able to choose the proper method of calculation for a studied system, one should know the limitation of each of the two methods for that specific system and to what extent is each method approved to be closer to the experimental values. Knowing the optimum U value can be reached empirically by applying different values of U for either GGA or LDA. From the following list of examples from literature, an assessment of the performance of different values of U when applied to both GGA and LDA for the same system can be realized:

• Griffin et al. [33] studied the FeAs crystal comparing the GGA+U and LDA+U levels of accuracy, using U_eff = −2 to 4 eV. The results showed that for the bond distances and angles in the crystal, the GGA+U gave results close to the experimental values when U ≤ 1 eV, whereas using LDA, the structural properties were poorly predicted. It was observed that increasing the value of U in the GGA+U increased the stabilization energy for antiferromagnetic ordering. Both GGA+U and LDA+U overestimated the value of magnetic moment. However, only the GGA+U could attain the experimental values of magnetic moment for negative U_eff [5].

• Cerium oxides (CeO₂ and Ce₂O₃) were tested by Christoph et al. [34] comparing GGA+U and LDA+U level of theory meanwhile studying the effect of the U_eff value on the calculated properties. It was found that the value of U_eff is dependent on the property under examination. The sensitivity toward U_eff values was especially high for properties of Ce₂O_3; because it has an electron in the 4f orbital, which is sensitive to the change in the effective on-site Coulomb repulsion due to the strong localization, in contrary to the CeO₂ that has an empty 4f orbital. The GGA+U showed an acceptable agreement with experiment at lower energies of U_eff than LDA did, with values of U_eff 2.5–3.5 eV for LDA+U and 1.5–2.0 eV for GGA+U, which can be due to the more accurate treatment of correlation effects within the GGA potential. On the other hand, the structural properties were better represented by the LDA+U method for CeO₂. Regarding Ce₂O₃ electronic structure, both LDA+U and GGA+U results showed a similarly good accuracy, while for the calculated reaction energies, LDA+U results showed better accuracy. [34]

• Sun et al. [35] studied PuO₂ and Pu₂O₃ oxides using both GGA+U and LDA+U methods. Although PuO₂ is known to be an insulator [36], its ground state was reported experimentally to be an antiferromagnetic phase [37]. For PuO₂, at U = 0, the ground state was a ferromagnetic metal, which is different from experimental results. Upon increasing the amplitude of U to 1.5 eV, the LDA+U and GGA+U calculations correctly predicted the antiferromagnetic insulating ground-state characteristics. For the lattice parameters, it was found that higher values of U (U = 4 eV) were needed with the LDA+U than for GGA+U. At U = 4 eV, it is expected that both LDA+U and the GGA+U would show a satisfactory prediction of the ground-state atomic structure of Pu₂O₃. However, the study showed that above the metallic-insulating transition, the reaction energy decreases with increasing U for the LDA and the GGA schemes. Therefore, for both Pu₂O₅ and PuO₂, the LDA+U and GGA+U approaches, with U as large as 6 eV, failed to describe the electronic structure correctly. When the energy gap increases, the electrons gain more localization that causes a difficulty of making any new reactions, consequently increasing the reaction energy. When U exceeds 4 eV, the conduction band electrons approximately considered to be ionized; thus, the atoms (cores or ions) have got a better chance to react with other atoms which, resulting in a reduction in the reaction energy.

As noticed in the previous studies, the U value is material dependent, besides being variable among the level of theory used. In general, the more localized the system is, the more sensitive it is to the value of U. The estimated value of U for a system of material using a specific level of calculation should not be extended to another system; rather, it should be recomputed each time for each material and even upon change of the level of calculation. Researchers will need to perform calculations using different U values within different xc functionals to get the best prediction of the calculated properties in comparison with the experimental measurements or with the other computational results as benchmark.

3.3. The effect of U on pure and defected systems

The chemical properties of transition-metal systems with localized electrons, mainly within d or f orbitals, are typically governed by the properties of the valence electrons. Experimentally, these electrons are observed to be localized in their orbitals due to strong correlations [38], whereas computationally, DFT approximated xc functionals tend to overly delocalize them while over-stabilizing metallic ground states, and thus underestimating the bandgap for semiconductors, and may reach false prediction of metallic behavior for systems like the Mott insulators. U can induce electronic localization due to the explicit account for the on-site electronic interactions. Another common problem that DFT calculations can impose is the prediction of the properties of materials with defects, as the underestimation of the bandgap by DFT can cause the conduction band (CB) or the valence band (VB) to kind of mask the true defect states. This is because defects can cause unpaired electrons and holes to form, which are overly delocalized by DFT, as it attempts to reduce the Coulomb repulsion due to self-interaction error.

Ref. [39] discusses an example of this problem studying the anatase TiO₂, where they showed that the description of the distribution of electrons in the unit cell, created from oxygen vacancies and hydrogen impurities, is wrongly predicted using GGA-PBE scheme of DFT calculations. In the case of oxygen vacancies, their calculations predicted a 2.6 eV bandgap, which is about 0.6 eV smaller than that reported experimentally. The electrons left in the system upon vacancy formation are completely delocalized over the entire cell. These electrons are incorrectly shared over all the Ti atoms of the cell, and as a result, the atomic displacements around the vacancies are predicted to be symmetric. All these findings indicate the difficulty of DFT methods to describe the properties of defects in wide bandgap metal oxides. Also, the accuracy of the description of the electronic structure of the partially reduced oxide systems was reviewed and discussed within the first principle methods [40].

The electronic structure of TiO₂–both pristine and doped–is one of the examples that is frequently studied in literature. Typically, in the anatase and rutile phases, computational studies encountered the problem of considerable underestimation of the bandgap, which presented a barrier in the prediction of further related properties. Titania is widely studied in various photoelectrochemical applications, and accurate theoretical assessment is required to be able to enhance its catalytic properties. In addition, to further improve the properties of TiO₂ as a photocatalyst, an optimization of the band structure is required, including narrowing the bandgap (Eg) to improve visible light absorption, and proper positioning of the valence band (VB) and the conduction band (CB) [41]. Efforts on narrowing the bandgap of the TiO₂ have been done through doping with metallic and nonmetallic elements that typically replace Ti or O atoms, and thus change the position of the VB and or the CB leading to a change in the bandgap [42]. In the following subsections, titania will be used as an example to assess the effect of U correction by presenting results from literature for both pristine and doped cases. We will monitor the behavior of the materials before and after U correction, while assessing the significance of the U correction for correct prediction of the material’s properties.

3.3.1. The Bandgap problem: pristine TiO₂ with U correction

Regarding the electronic structure of titania, the bandgap was underestimated by the standard DFT, while found to be overestimated when the hybrid functional Heyd-Scuseria-Ernzerhof (HSE06) was applied. However, the bandgap prediction was markedly improved by adding the Hubbard U correction. The obtained band structures using GGA-PBE showed bandgaps of 2.140 and 1.973 eV for anatase and rutile, respectively. However, upon applying the localization of the excess electronic charge using +U correction, the predicted bandgaps are accurate and in a good agreement with the experimental and the computationally expensive hybrid functional (HSE06) results [43], Figure 2. In another study, for rutile TiO2, the prediction of the experimental bandgap is achieved with a U value of 10 eV, whereas the crystal and electronic structures were better described with U < 5 eV [19].

Figure 2.

Total DOS of pure TiO₂ (anatase and rutile) [43]. Copyright (2014)—American Chemical Society.

Dompablo et al. [19] compared the effect of the U parameter value (0 < U < 10 eV) within the LDA+U and the GGA+U on the calculated properties of anatase TiO₂. Both LDA+U and GGA+U required a small value of U (3 and 6 eV, respectively) to reproduce the experimental measurements, Figure 3. However, using very large U values leads to mismatching, where the lattice parameters (a and c) and the volume of unit cell are increasing with increasing U, due to the Coulomb repulsion increase. Note that standard DFT and the hybrid functional HSE06 failed to calculate the crystal lattice.

Figure 3.

Calculated c/a ratio vs. U value. Reproduced from [19], with the permission of AIP Publishing.

On the other hand, the calculated bandgap within the GGA+U and LDA+U methods was found to be in better agreement with experiments compared to the conventional GGA or LDA, with small difference in the required U value. The bandgap was shown to increase by increasing the U value till 8.5 eV, which gave a result close to the experimental bandgap and in agreement with those obtained in previous DFT studies [44]. For values of U larger than 8.5, the bandgap was overestimated. It is worth to mention that this value (8.5 eV) is considered high when compared to other U values for other transition metal oxides [29]. In all these calculations, the Hubbard U parameter was used for the d or f orbitals of transition metals. However, when the Ti-O bonding is considered, while applying the correction only to the 3d-states, it can be estimated that this correction might have an influence over the Ti−O covalent bonding, where the Ti states are shifted and the 2p states of oxygen are not changed [45]. In this regard, several first principle calculations were derived to study the electronic, structural, and optical properties of TiO₂ polymorphs by applying the U correction for the oxygen’s 2p orbitals and titanium’s 3d orbitals [46]. In order to correct the bandgap, while avoiding the use of large U values and the bonding problem, Ataei et al. [47] reported that with values of 3.5 eV for both O 2p and Ti 3d-states, the results for the lattice constants, bandgap, and gap states are in good agreement with the experimental reports.

3.3.2. Doped-TiO₂ with U correction

In a recent study [40], a comparison was performed to elucidate the effect of different U values in representing the bandgap states produced by interstitial hydrogen atom and oxygen vacancy within the bulk Ti anatase structure. The dependence on the method used was observed, beside the value of U within GGA+U scheme, see Figure 4. When the U correction was not applied, the bandgap is underestimated, as expected, and the electrons caused by the oxygen vacancy or the hydrogen impurity are fully delocalized and have conduction band character. Upon applying the U correction, the states start to localize and are became deeply localized in the gap with increasing the value of U. In all these calculations, the Hubbard U correction was applied for the Ti 3d orbitals only; by applying the correction for the 2p oxygen orbitals with U = 3.5, the results were in agreement with previous results [47].

Figure 4.

Total and partial density of states DOS for anatase TiO₂ doped with (oxygen vacancies and interstitial H atom) obtained with GGA-PBE+U. Adapted with permission from [39]; (a) U = 0, (b) U = 2, (c) U = 3, (d) U = 4.

The intrinsic defects in TiO₂ (vacancies) have been computationally studied, providing a fast-cheap method to guide researchers in choosing the defect position in the solid crystal. The oxygen vacancy in the rutile crystal was investigated [48] using the DFT+U with U value of 4.0 eV, indicating that oxygen vacancy in the rutile crystal introduces four local states with two occupied and two unoccupied states, with no change in the bandgap (2.75 eV) [48]. The Ti vacancy effect on the bandgap (E_g) was also studied [49] using the GGA+U with U values of 7.2 eV. It was found that Ti vacancy caused ferromagnetism besides widening the valence band, and switching the TiO₂ from n-type to p-type semiconductor with higher charge mobility [49].

4. Modeling of organometallics using the (+U)

Hubbard correction is a computational tool that can be applied widely not only to crystals but also to the strongly correlated metals attached to other noncorrelated systems such as organic moieties. One of these important systems is metal organic framework (MOF).

4.1. Metal organic frameworks (MOFs)

MOFs are crystalline nanoporous materials where a centered transition metal is linked to different types of ligands, which provide a very large surface area [50] that can allow their use in supercapacitors and water splitting applications. Most of the MOFs have open metal sites, which are coordinative unsaturated metal sites with no geometric hindrance. While the whole material remains as a solid, the structure allows the complex framework to be used in gas capturing and storage, and the binding energy between the MOFs and the gas or water molecules allows the prediction of the capturing mechanism. The cage shape of the MOFs and the organic moiety allow their use in many applications such as drug delivery and fertilizers, while the magnetic behavior of MOFs allows the researchers to correctly predict how it can be used in applications. Quantum mechanics frame of work is usually used to describe the full interaction between the centered metal ion and the surrounding ligands, due to the fact that the synthesis of these materials is both time and money consuming. The complex geometry resulted from the computational calculations is important to predict the small change in electronic structure upon application of external stimuli [51].

Density functional theory (DFT) has been used to model the MOFs as it allows the “mapping” of a system of N interacting electrons onto a system of N noninteracting electrons having the same ground state charge density in an effective potential. However, DFT fails to describe electrons in open d- or f-shells [8]. The pure DFT calculations usually wrongly estimate the bandgap and the ferromagnetic (FM) or antiferromagnetic (AFM) coupling for the centered metal in the MOFs. The reason for this wrong estimation is the localized spin and itinerant spin density that are coupled via the Heisenberg exchange interaction [52, 53]. In this interaction, the ferromagnetic sign is assumed if the hybridization of the conduction electrons (dispersive LUMO band) with a doubly occupied or empty d orbital of the magnetic center is sufficiently strong. Owing to Hund’s rule, in the d shell, it is energetically favorable to induce spin polarization parallel to the d-shell spin. The itinerant spin density, however, forms at an energy penalty determined by the dispersion of the conduction band; the larger the density of states at the Fermi level, the easier is for the itinerant spin density to form. The addition of an extra interaction term that accounts for the strong on-site coulomb U correction has proved to lead to good results [54]. One more advantage of the DFT+U is that it can be used to model systems containing up to few hundred atoms [55]. The U parameter affects the predicted electronic structure and magnetic properties; in the following paragraphs, we will discuss some of the MOF applications and how to fit a proper magnitude of U in DFT+U calculations:

• The magnetic properties of the MOF of the complex dimethyl ammonium copper format (DMACuF) were predicted correctly [56] using the (GGA+U) with convenient U values (U = 4–7 eV) for Cu 3d-states to describe the effect of electron correlation associated with those states. Also, the magnetic properties of MOFs of TCNQ (7,7,8,8-tetracyanoquinodimethane) and two different (Mn and Ni) 3d transition metal atoms were predicted correctly without synthesizing. But in this case to properly describe the d electrons in Ni and Mn metal centers, spin-polarized calculations using the DFT+U with U value of (U = 4 eV) were performed [53]. It can be claimed that the varying of U in the range of 3 to 5 eV does not appreciably change the values of the Ni and Mn magnetic moments, nor the corresponding 3d level occupations, in particular, that of the Ni (3dxz) orbital that crosses the Fermi level [53].

• The binding energy of CO₂ to a Co-MOF-74 was predicted [57] using DFT+U with U values (U = 0–6 eV), and it was found that the value of U between 2 and 5 eV gives lattice parameters matching with experiment due to the fact that the Co-O bond length decreases with U, since U localizes the Co d-states, which allows the CO₂ molecule closer to the charged open metal site, increasing the electrostatic contribution to the binding energy [57].

• The Cu-BTC [58], a material consisting of copper dimers linked by 1,3,5-benzenetricarboxylate C₆O₉H₃ (BTC) units, was studied for its ability to absorb up to 3.5 H₂O per Cu as the Cu binds to the closest oxygen of the water molecule [59]. The U parameter in the meta-GGA+U calculation of the Cu-BTC was adjusted with the experimental crystallographic structure and the bandgap by minimizing the absorption at 2.3 eV. The U values gave the best results at 3.08 eV for Cu and 7.05 eV for O because those values reduced the calculated root mean square residual forces on the ions at their experimental fixed positions to its minimum value. The nonzero U of oxygen greatly reduces the residual forces, while the value of U for Cu ions controls the splitting in the Cu d levels, which have a great effect on calculated bandgap [59].

4.2. Spin-crossover (SCO)

Spin-crossover (SCO) is a unique feature in which the centered transition metal ion linked to the surrounding ligand has the ability to attain different spin states with different total spin quantum numbers (S), while keeping the same valence state [47]. This property allows MOFs and organometallics generally to reversibly switch between spin states upon application of temperature, pressure, light, or magnetic field, such as changing between low spin and high spin [60, 61]. The SCO can be predicted effectively using the U correction as well as the effect of temperature on the SCO. The use of DFT+U to model SCO was first done by Lebègue et al. [62]. SCO is generally appealing for metals that have availability to change between high spin and low spin due to the small difference between the HOMO and LUMO levels. Iron (Fe) is one of the most important examples for this property. Fe is important since it can be found in many ores and can be used in many applications such as solar cells. Besides, Fe can be called the source of life, which is the heme molecule in the blood and which is responsible for the transfer of oxygen and carbon dioxide to and from the cell, respectively, consist of Fe-porphyrin molecule. Modeling of those molecules and their reaction mechanisms provides information about drug reactions inside the blood stream. Unfortunately, the common exchange-correlation functional fails to predict the properties of the deoxygenated active site of hemoglobin and myoglobin and Fe-porphyrin molecules [63, 64]. Some of the examples of SCO in Fe complexes are listed below:

• The SCO of [TiFe(CN)₆]²⁻, [CrFe(CN)₆]²⁻, [MnFe(CN)₆]²⁻, and [CoFe(CN)₆]²⁻ frameworks has been studied [65] using the DFT+U. It was found that high U values > 8 eV should be applied to the low spin Fe site, while low U value should be applied to the high spin ion. The results showed a great agreement with other DFT calculations. The generally used DFT-GGA failed to predict the high spin of the five coordinate Fe complexes [68], but it could be obtained by the DFT+U with U~ 4 eV and J~ 1 eV. The complexes Fe(phen)₂(NCS)₂ and Fe(btr)₂(NCS)₂ were tested using a U value of 2.5 eV [62], with the energy difference between the low spin state and high spin state is in agreement with the experimental values and proved that the U coulomb term was needed. The study showed the importance of magneto elastic couplings through the correlation between the spin state and the structure [62].

• Another study on the complex [Fe(pmd)-(H₂O)M₂(CN)₄].H₂O (pmd = pyrimidine and M = Ag or Au) showed an interesting SCO behavior according to temperature [66]. This complex forms chain polymers that contain two different Fe(II) ions Fe1 and Fe2. Through hydration/dehydration, temperature changes between 130 and 230 K for the Ag-based coordination polymer changing the SCO reversibly and this change is due to the structure change caused by the water molecules in the network. For the Au-based complexes, only the SCO transition was different in the hydrated framework [66]. Such behavior could be explained using the DFT+U calculations [67]. The low spin-high spin transition was found to occur only on the sixfold nitrogen coordinate Fe₁ ion, while the Fe₂ ion coordinate with four nitrogen and two oxygen from the water molecules. For the dehydrated compounds, the effect of the Au atom caused a difference in the degree of the covalent bonding, which resulted in a distinct behavior of the Au network as compared to the Ag network. The hydrated and dehydrated Ag networks were predicted to exhibit a low spin-high spin transition, whereas the dehydrated Au network was predicted to remain in a high spin state [55].

• Fe-porphyrin molecules were found to have an intermediate spin state. The ground-state configuration was indicated to be (d_XY)²(d_π)²(d_z²)² using Mossbauer [68, 69, 70], magnetic [71] and NMR [72, 73] measurements. However, Raman spectroscopy predicted a ground state with a configuration (d_XY)²(d_π)³(d_z²)¹ [74]. Therefore, a computational calculation was important to predict the reason for those results. The DFT+U was used to predict the electronic structure and magnetic properties of Fe molecules for a range of Coulomb U parameters (U = 2–4 eV), which is reasonable for iron [10, 75], and then compared to available data in literature. It was found that GGA+U with U value of 4 eV provided an overall better comparison of the structural, electronic, magnetic properties, and energy level diagram of these systems [76, 77].

To summarize, DFT+U are good to predict the correlation in the centered metal in organometallics. The spin change between FM and AFM states or in SCO can all be well predicted by the Hubbard correction, while the pure DFT fails due to the correlation in the d or f orbitals of the centered metal.

5. Solving the CO adsorption puzzle with the U correction

Studying surface chemistry is of great significance for enhancing the overall efficiency of many electrochemical applications [78, 79, 80]. In catalysis, for example, understanding the adsorption mechanism of species on catalytic surfaces—mainly electrodes—is essential in order to formulate a design principle for the prefect catalyst that can reach the optimum efficiency for a desired electrochemical process [81, 82, 83]. Typically, the adsorption of CO on metal surface is widely acknowledged as the prototypical system for studying molecular chemisorption [84, 85, 86, 87]. Despite the extensive experimental studies, grasping the complete theoretical description of the “bonding model” has not yet been reached, due to the inability of experimental tools to fully describe the details of molecular orbital interactions and to make a profound population analysis, which is based on studying the electronic structures of the substrate and surface particles [88, 89]. To this end, DFT can be utilized to explicitly describe electronic structures of the system particles in greater details, which can help in extending the conceptual model of CO chemisorption [90, 91, 92, 93, 94]. Unfortunately, due to the inherent wrong description of the electronic structure by DFT, wrong predictions of CO preferred adsorption sites are observed that contradict experimental results, especially for the (111) surface facets of transition metals, leading to the so-called “CO adsorption Puzzle” [95, 96]. The root of this DFT problem resides on the fact that both local density and generalized gradient approximation functionals underestimate the CO bandgap, predicting wrong positions of the CO frontier orbitals, which results in an overestimated bond strength between the substrate and surface molecules [97].

One of the popular solutions that has been exploited by researchers to resolve the adsorption site prediction puzzle is the DFT+U correction [97, 98]. In this approach, the position of the 2π* orbital is shifted to higher values, by adding the on-site Coulomb interaction parameter. By doing so, the interaction of CO 2π* orbital with the metallic d-band will no longer be overestimated, bringing the appropriate estimation of the CO adsorption site. Kresse et al. [99] first implemented this method and successfully obtained a site preference in agreement with experiment, emphasizing that the use of such a simple empirical method is able to capture the essential physics of adsorption. DFT calculations utilizing GGA functionals predict adsorption on the threefold hollow site for Cu(111) and in the bridge site on Cu(001), instead of the experimental on-top site preference. Reference [98] implanted Kresse’s method to investigate the adsorption of CO on Cu(111) and (001) surfaces with 1/4 monolayer (ML) coverage on different adsorption sites. In that study, the HOMO-LUMO gap of the isolated CO molecule was demonstrated to be increased by increasing the value of U. Also, upon changing the U value, the corresponding adsorption energies of the CO over the different adsorption sites were calculated.

Reviewing the Cu (111) surface results, five different U values (0.0, 0.5, 1.0, 1.25, and 1.5 eV) were used in the calculations. It was observed that only 20 meV changes in the adsorption energy (higher coordinated hollow sites) for U = 1.25 and 0.03 eV for U = 1.5 eV. Nonetheless, the absolute value of adsorption energy decreases linearly with increasing U, where the rate of reduction is found to be larger for higher coordinated sites. It was observed that the site preference between top and bridge sites to be reversed around the U value of 0.05 eV, while between the top and hollow sites around U = 0.45 eV. Concerning the adsorbate (surface) description in the study, the calculated interlayer relaxations were the same as that calculated using the GGA (PW91) functional without the U correction. Not only does the U correction help in solving the adsorption puzzle dilemma, but it can also enhance the description of other related properties, such as the calculated work function and the vibrational spectra for the CO-metal complexes, which are also demonstrated in Ref. [98] (Figure 5).

Figure 5.

A schematic sketch of the molecular eigenstates of the CO molecule. The DFT+U technique shifts the LUMO orbitals to higher energies, but the energies of the occupied orbitals remain the same.

6. Summary and outlook

In this chapter, the corrective capability of the DFT+U is overviewed and evaluated for a number of different classes of materials. Generally, the addition of the on-site Coulomb interaction potential (U) to the standard DFT Hamiltonian proved to provide significant changes to the predicted electronic structures, which can solve the inherent DFT bandgap prediction problem. The value of U can either be theoretically calculated or semiempirically tuned to match the experimental electronic structure. For the various case studies and applications reviewed, the criticality of correcting the electronic structure predictions was manifested, as it leads to significant improvements for the prediction of further electronic-related properties. Prior to the practical assessment, the theoretical foundation of the DFT+U method is briefly discussed and is verified to be rather simple adding only marginal computational cost to the standard DFT calculations. Compared to other corrective approaches, the DFT+U formulation demonstrated to be simpler in terms of theoretical formulation and practical implementations with considerably lower computational cost, while having nearly the same predictive power; it can even capture properties of certain materials that cannot be captured by other higher level or exact calculations. One of the most popular implementations of the U correction is the description of the electronic structure for strongly correlated materials (Mott insulators). The behavior of these types of insulators cannot be captured by applying Hartree-Fock, band theory based, calculations, as the root of this problem resides on the deficiency of the band theory to capture such behavior, as it neglects the interelectron forces. One of the simple models, which explicitly accounts for the on-site repulsion between electrons at the same atomic orbitals, is the Hubbard model. Based on this model, the DFT+U method is formulated to improve the description of the ground state of correlated systems.

The theoretical and semiempirical techniques of the U optimization are discussed. The semiempirical 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. However, the semiempirical evaluation of U does not permit the capturing of changes in the on-site electronic interaction under changing physical conditions, such as chemical reactions. The practical implementations of U correction are discussed, while assessing the effect of the DFT scheme employed and the calculation parameters assigned on the numerical value of the optimum U utilized. The corrective influence of the U correction is validated by reviewing different examples and case studies in literature. Starting with the transition metal oxides, the effect of adding the U parameter to correctly describe the electronic structure of pure and defected TiO2 is reviewed, showing the different optimum values of U utilized for each level of calculation. Then, the implementation of the Hubbard correction to the systems that comprises molecules with solid-state crystals is reviewed, such as organometallics. The addition of U to the DFT calculation provides a better understanding of the behavior of the metals inside the organometallic systems. One of the most importantly studied organometallic systems is the metal organic frameworks (MOFs). Different examples in literature are reviewed, showing the effect of the U correction and how it can significantly improve the prediction of the magnetic properties of such systems. Also, one of the unique features of organometallics, which can be influenced the U correction, is the spin crossover (SCO). This property allows the MOFs and the organometallics generally to reversibly switch between spin states upon changing the external parameters. The SCO is proved to be predicted more effectively by applying the U correction, as demonstrated in the results presented in literature. Finally, the significance of the DFT+U method is manifested upon describing the adsorption mechanism of CO on transition metal systems. The influence of U correction on solving the so-called adsorption Puzzle is demonstrated, which leads to the correct prediction of CO adsorption site preference, which was an unresolved problem when DFT calculations are applied alone.

Upon reviewing the presented applications and different case studies, where the U correction significantly improved the estimated results without changing the essential physics of the systems, we can estimate the potential of the Hubbard correction to gain a greater weight in the future of computational chemistry. Despite the convenience of the semiempirical tuning of U, the capabilities of the Hubbard correction in this way cannot be fully exploited, as it cannot be used to study systems with variations of on-site electronic interactions. On the other hand, despite the availability of theoretical U calculation methods, their computational costs are considerably large, compared to the semiempirical methods. Therefore, further improvements to the ab initio calculation of U is still required, with lower computational costs, in order to conceive full potential of the U correction that is able to capture phase changes and chemical reactions for the studied physical systems.

Acknowledgments

This work was made possible by NPRP Grant no. NPRP 6-569-1-112 from the Qatar National Research Fund (a member of Qatar Foundation). The statements made herein are solely the responsibility of the authors.