DFT+U: a Hubbard Hamiltonian-based correction for modeling localized electronic states

DFT+U: uma correção baseada no Hamiltoniano de Hubbard para modelar estados eletrônicos localizados

Abstract

Density Functional Theory (DFT) is a widely used method to calculate the electronic structure of materials. One of the pivotal ingredients in the traditional DFT modeling strategy is the Exchange and Correlation Functional ($E_{XC}$). An inappropriate choice of the $E_{XC}$ can lead to substantially erroneous predictions. For example, the Generalized Gradient Approximation (GGA) dramatically fails to predict the band gap of strongly correlated systems. DFT$+U$, a corrective methodology based on the Hubbard model, offers a cost-effective approach to address this problem. This work provides an introductory overview of the DFT$+U$ method. First, we discuss the importance of the $E_{XC}$ functional in DFT, and then we introduce the DFT+U approach. Subsequently, we illustrate the method’s applicability by calculating the electronic, magnetic, and structural properties of FeO using PBE and PBE+U functionals. We show that GGA fails to describe the electronic character and magnetic ordering of FeO, while DFT$+U$ provides an accurate description of the FeO properties.

Keywords:
Hubbard Model; Strongly Correlated Systems; Density Functional Theory

Resumo

A Teoria do Funcional da Densidade (DFT) é um método amplamente utilizado para calcular a estrutura eletrônica de materiais. Um dos ingredientes centrais da estratégia tradicional de modelagem na DFT é o funcional de troca e correlação ($E_{XC}$). Uma escolha inadequada do $E_{XC}$ pode levar a previsões substancialmente incorretas. Por exemplo, a aproximação do Gradiente Generalizado (GGA) falha drasticamente em prever o band gap em sistemas fortemente correlacionados. O método DFT$+U$, uma metodologia corretiva baseada no modelo de Hubbard, oferece uma abordagem computacionalmente eficiente para lidar com esse problema. Este trabalho apresenta uma visão introdutória do método DFT$+U$. Primeiro, discutimos a importância do funcional $E_{XC}$ na DFT e, em seguida, introduzimos a abordagem DFT$+U$. Posteriormente, ilustramos a aplicabilidade do método por meio do cálculo das propriedades eletrônicas, magnéticas e estruturais do FeO usando os funcionais PBE e PBE$+U$. Mostramos que a GGA falha em descrever corretamente o caráter eletrônico e a ordem magnética do FeO, enquanto a DFT$+U$ fornece uma descrição precisa das propriedades do FeO. Palavras-chave: Modelo de Hubbard, Sistemas Fortemente Correlacionados, Teoria do Funcional da Densidade.

Palavras-chave:
Modelo de Hubbard; Sistemas Fortemente Correlacionados; Teoria do Funcional da Densidade

1. Introduction

Interactions between charged particles dictate a material’s properties. Therefore, electronic structure calculations play a crucial role in materials science and technology; they provide the possibility to describe, understand, and even design materials for device applications, for example [¹, ², ³]. In general, computational electronic structure methods attempt to tackle the quantum many-body problem, which can rarely be solved exactly; thus, approximations are commonly used.

Strategies based on Density Functional Theory (DFT) balance accuracy and computational effort in solving the many-body problem. The formalism underlying DFT provides an exact framework to obtain the ground state (GS) of a many-body system [⁴]; nevertheless, in practice, approximations are necessary, as discussed later in this paper [⁵]. One central aspect of DFT is using the electronic density ($n$) as the key quantity, contrasting with the wavefunction methods, which directly attempt to obtain the system’s wavefunction.

In many-body wavefunction methods, we describe a system of $N$ interacting electrons by constructing a wavefunction that depends on the coordinates of all electrons simultaneously. This wavefunction lives in a high-dimensional space: since each electron is characterized by three spatial coordinates, the full wavefunction is defined in a $3N$-dimensional space. We typically use combinations of single-particle wavefunctions to build this many-electron wavefunction.

For example, linear combinations of Slater determinants can ensure the overall wavefunction satisfies the antisymmetry required by the Pauli exclusion principle. Once this object is constructed, all physical properties of interest–such as the total energy, density, or response functions–can, in principle, be obtained by computing appropriate expectation values. However, due to the exponential scaling of the complexity with the number of electrons, this approach quickly becomes computationally demanding, limiting its practical use to relatively small systems [⁵, ⁶].

DFT takes a fundamentally different approach: although we still solve an effective set of one-body equations, the single-particle wavefunctions are used not to reconstruct a $3N$-dimensional many-body wavefunction, but rather to build the electronic density $n\left(\mathit{r}\right)$, a scalar quantity defined by three spatial variables. Observables are then evaluated as functionals of this much simpler object. Replacing the many-body wavefunction with the electronic density captures essential quantum mechanical effects and makes DFT significantly more practical and scalable for realistic materials calculations. More details can be seen in the appendix A.

From the practical standpoint, DFT requires approximations for the Exchange and Correlation (XC) energy functional [⁵, ⁷] ($E_{XC}\left[n\left(\mathit{r}\right)\right]$). Defining a universal $E_{XC}\left[n\left(\mathit{r}\right)\right]$ is an open problem in DFT; thus, we need to add the right ingredients to the $E_{XC}\left[n\left(\mathit{r}\right)\right]$ to describe a particular class of materials [⁷]. The most common methods to express $E_{XC}$ are the Local Density Approximation (LDA) and the Generalized Gradient Approximation (GGA). Both approaches are successful in describing certain sets of materials. However, LDA and GGA provide inaccurate descriptions of strongly correlated materials [⁸], systems in which the potential energy from electron-electron interactions is comparable to, or can even surpass, the kinetic energy of the electronic particles [⁹]. These materials typically present open-shell $d$ or $f$ localized states, where the electrons interact strongly, and correlation effects become crucial [²⁴, ⁹].

The failure to account for strong correlations leads to clear signatures, such as predicting Transition Metal Oxides (TMOs) to be metallic when, in fact, they are insulators due to the localized nature of the partially filled $d$ states [⁸]. TMOs are prototypical examples of Mott insulators, where the insulating behavior arises from electron correlation [⁸]. The lack of success is related to both approximations relying on mean-field strategies, which leads to an excessive delocalization of valence electrons [¹⁰]. Therefore, modeling systems with prominent electron localization requires corrective methods. In this work, we discuss DFT+$U$: an approach to addressing strongly localized valence electrons in DFT [¹¹, ¹²]. The parameter $U$, inspired by the Hubbard model, introduces an effective on-site Coulomb repulsion to adjust the GGA or LDA inherent over delocalization. We apply the DFT+$U$ to the Mott insulator FeO and show how the $U$ correction overcomes the band-gap problem.

2. The XC Functional Role in DFT

2.1 The energy density functional

We can model atoms, molecules, and solids as a group of ions and electrons interacting via Coulomb forces. Suppose we want to determine a material’s total energy ($E$). In that case, we need to compute the electrons’ total kinetic energy ($T_e$), the electron-electron ($V_{ee}$), and the ion-electron ($V_{ext}$) interaction energies. From this point on, we adopt the Born-Oppenheimer (BO) approximation, which assumes that the electronic cloud adjusts instantaneously to the much slower motion of the nuclei, such as vibrations and rotations. This separation of timescales, justified by the significant mass difference between electrons and nuclei, allows electronic and nuclear motions to be treated separately. Within this framework, the nuclei are handled semi-classically using the Hellmann–Feynman theorem, where the forces acting on them are derived from the electronic energy [⁵].

The BO approach enables obtaining equilibrium structures and structural relaxations while retaining a quantum mechanical description of the electrons. The electronic state depends parametrically on the positions of the nuclei. From a practical standpoint, the ion-ion interaction energy, being independent of the electronic wavefunction, varies parametrically with the ionic coordinates. This is essential in the Hellmann–Feynman scheme for determining the equilibrium positions of the atoms.

Within the DFT framework, $E$ is expressed as a functional of the electronic density $n\left(\mathit{r}\right)$, and the electronic density $n\left(\mathit{r}\right)$ goes from an observable to a key object or variable [⁵, ¹³, ⁷]. All properties are fully described through an energy functional relation:

(1) $$E\left[n\right]=T\left[n\right]+V_{ee}\left[n\right]+V_{ext}\left[n\right]=F\left[n\right]+V_{ext}\left[n\right],$$

the term $F\left[n\right]=T\left[n\right]+V_{ee}\left[n\right]$ is the same for all electronic systems, and it is referred to as the universal functional. The $V_{ext}$ varies for each system, and unlike $F\left[n\right]$, it does not depend exclusively on the electrons; thus, it is commonly named as the external potential.

The Kohn-Sham (KS) method obtains the universal functional by creating a non-interacting auxiliary system (KS ansatz) with the same density as the original one [⁵]. From the practical point of view, we define a KS universal functional as:

(2) $$F\left[n\right]=T_s\left[n\right]+E_H\left[n\right]+E_{XC}\left[n\right],$$

the functional $T_s\left[n\right]$ is the total electron kinetic energy of the auxiliary system. $E_H$ is the Hartree interaction energy, which comes from the solution to the classical Poisson equation [¹⁴], thus including only classical effects on the interaction between electron densities. The exchange and correlation functional $E_{XC}\left[n\right]$ contains all interactions omitted from the sum $T_s\left[n\right]+E_H\left[n\right]$, so that equation 2 provides an exact solution for the universal functional. Despite $E_{XC}\left[n\right]$ accounting for a small portion of the total energy, it is crucial to model chemical bonds as commonly referred to as nature’s glue [¹⁵]. By applying the rules of DFT, as expressed in appendix A, one arrives at the Kohn-Sham (KS) equations:

(3) $$\left[-\frac{{\hslash }^2}{2m}{\nabla }^2+v_{eff}\left(\mathit{r}\right)\right]{\varphi }_i\left(\mathit{r}\right)={ϵ}_i{\varphi }_i\left(\mathit{r}\right),$$

with

$$v_{eff}\left(\mathit{r}\right)=v_{ext}\left(\mathit{r}\right)+v_H\left(\mathit{r}\right)+v_{XC}\left(\mathit{r}\right),$$

${\widehat{h}}_{KS}=-\left({\hslash }^2/2m\right){\nabla }^2+v_{eff}$ is the KS hamiltonian, ${\varphi }_i\left(\mathit{r}\right)$ are the KS orbitals and ${ϵ}_i$ are the KS eigenvalues (see appendix A).

2.2 Exchange and correlation functional

The energy functional in equation 1 provides a path to calculate materials’ properties; this approach’s effectiveness depends on how accurately we express the terms in equation 2. All many-body quantum-mechanical effects must be included in $E_{XC}$, but the exact mathematical expression of $E_{XC}\left[n\left(\mathit{r}\right)\right]$ is unknown [¹³]; therefore, it is natural that KS-DFT misses correlation effects in practice [⁵].

There are essentially two main approaches to construct a density functional [¹⁶]. The first relies on mathematical modeling of known physical effects, aiming to capture fundamental properties of the electronic system through approximations and norms grounded in theory [¹⁷]. The second approach is data-driven. For example, machine-learning-based density functionals are trained on large datasets of reference calculations, and the underlying patterns of electronic behavior are extracted directly from data rather than from explicit physical assumptions [¹⁷, ¹⁸]. Regarding the functionals based on physical effects, there are levels of approximations for the $E_{XC}$ term, usually disposed of in the so-called Jacob’s ladder of functionals, as seen in Figure 1, inspired by reference [¹⁹]. The higher the rung, the more accurate the functional, which typically leads to a higher computational cost.

Figure 1
Jacob’s ladder illustrating the various types of approximations for the term $E_{XC}$. The lowest part of the ladder is Hartree theory, which does not account for exchange and correlation effects. The highest part represents an exact $E_{XC}$.

The simplest approximation is based on the uniform electron gas and is the $1^{st}$ Jacob rung [²⁰, ⁵]. It is called the Local Density Approximation (LDA). Its mathematical form depends only on the density as a function of position $n^{\sigma }\left(\mathit{r}\right)$, with $\sigma$ representing the spin degree of freedom for spin-polarized calculations. Due to its simplicity, it is computationally cheap but only describes a limited group of systems appropriately. Since it is formulated based on the uniform electron gas, it is expected to give reasonable results with solids well described by a homogeneous electron gas, such as the alkali metals [⁵, ¹].

The GGA approach consists of adding the dependency of $\left|\mathbf{\nabla }{n}^{\sigma }\left(\mathit{r}\right)\right|$ to the LDA functional. It represents the $2^{nd}$ Jacob rung and an increase in accuracy since the contribution to the energy is not only the density itself anymore, but also its spatial variation. Even though it is a slightly computationally more expensive approach than the LDA, its cost-benefit is better [²¹].

Meta-GGAs are in the $3^{rd}$ rung and rely on the dependency on the kinetic energy density, $\tau$, or on the laplacian of the electron density ${\nabla }^2{n}^{\sigma }\left(\mathit{r}\right)$ [¹⁷]. They are more accurate and computationally more expensive than GGA functionals. If one considers $E_{XC}$ as a Taylor expansion, it adds more terms to the series to decrease inaccuracies. In this way, the more terms one adds, the more accurately $E_{XC}$ describes many-body quantum effects.

Hybrid functionals carry this terminology because they use a portion of the exact Hartree-Fock exchange energy. Therefore, they are non-local and generally more accurate functionals for predicting e.g. band gaps and formation energies when compared to GGA or LDA [⁵, ²²]. Nonetheless, they are much more expensive computationally [⁵, ²³].

Since all $E_{XC}$ used are in practice approximations, KS-DFT can miss part of the correlation effects in some systems. This can have high implications for strongly correlated systems, such as TMOs. DFT+$U$ stands as a cost-effective approach to modeling strongly correlated systems. It introduces a corrective parameter $U$ to consider the electron localization. This discussion raises a natural question: when does one use a hybrid functional or DFT+$U$? A natural answer is that it depends on what one wants to calculate. The hybrid functional approach is computationally more expensive than the $+U$ correction [²³, ⁸, ²⁵]. On the other hand, one needs to find an effective route to tune the $U$ parameter [²⁶].

2.3 The Hubbard XC

DFT+$U$ is based on the Hubbard model, which is essentially an extension of the tight-binding approximation with an on-site Coulomb repulsion interaction [²⁷]. The model is for non-relativistic fermions with spin $1/2$, electrons in this case, and it is summarized by the following Hamiltonian [²⁸]:

(4) $$\widehat{H}=-t\sum _{⟨i,j⟩\sigma }\left({\widehat{c}}_{i\sigma }^{†}{\widehat{c}}_{j\sigma }+{\widehat{c}}_{j\sigma }^{†}{\widehat{c}}_{i\sigma }\right)+U\sum _i{\widehat{n}}_{i\uparrow }{\widehat{n}}_{i\downarrow }.$$

${\widehat{c}}_{i\sigma }^{†}$ represents the creation operator, which creates a fermion with spin $\sigma$ on the site $i$, ${\widehat{c}}_{j\sigma }$ annihilates a fermion in site $j$ with spin $\sigma$, ${\widehat{n}}_{i\sigma }={\widehat{c}}_{i\sigma }^{†}{\widehat{c}}_{i\sigma }$ is the number operator, and it counts the number of electrons on-site $i$ with spin $\sigma$. The Hamiltonian is divided into the kinetic energy term, the first term on the Right-Hand Side (RHS), and a potential term proportional to the Hubbard $U$, which is the second term on the RHS. The Hubbard $U$ is the key parameter in the DFT+$U$ description.

In equation 4, $t$ represents the nearest-neighbor hopping amplitude, while $U$ quantifies the on-site Coulomb repulsion between electrons. These two terms compete with each other, as illustrated in Figure 2: for an electron to hop to a neighboring site, it must overcome the energy cost associated with $U$. Of particular interest is the regime where $U$ is comparable, or bigger, to the bandwidth $W$, which is determined by the hopping parameter $t$. This intermediate regime is interesting for localized states–such as those involving $d$ and $f$ electrons–and is characteristic of strongly correlated systems. Because both terms are comparable in energy, perturbation theory is not a good method within this limit, and more complex approaches are needed.

Figure 2
Electrons are allowed to hop between neighboring sites with probability amplitude $t$, but double occupation of a site is energetically costly, incurring a penalty of $U$.

From the Hubbard Hamiltonian in the $U/t>>1$ limit, one can obtain a Heisenberg spin-spin interaction for an antiferromagnetic (AFM) system [²⁹]. This result is applicable to several transition-metal monoxides, which are strongly correlated systems and typically exhibit an AFM ground-state magnetic phase [⁸]. Therefore, the idea behind DFT+$U$ is not to replace the underlying exchange-correlation functional–whether it be LDA, GGA, or a more elaborate approximation–but to correct it by adding an explicit Hubbard-like term. This correction is applied selectively to the chosen manifold of localized states, typically the d and/or f electrons, while the delocalized s and p electrons remain described by the standard Kohn–Sham Hamiltonian. This approach has proven especially useful for improving the description of magnetic properties in strongly correlated materials, as discussed above.

This approach is guaranteed by the potential below [³⁰], which is added explicitly on DFT equations (more precisely on $V_{eff}$ of equation 3, shown in the appendix):

(5) $${\widehat{v}}_{XC}^{U\sigma }=U\sum _{im}\left(\frac{1}{2}-n_m^{i\sigma }\right)|m⟩⟨m|.$$

$U$ is the Hubbard parameter. The $n_m^{i\sigma }$ term is the occupation matrix eigenvalue on site $i$, localized state $|m⟩$ ($d$ and/or $f$ states) and spin state $\sigma$. This potential acts only on localized states $|m⟩$ because of the projection operator written as a function of localized states ${\widehat{P}}_m=|m⟩⟨m|$. A diagram of its action in the DFT+$U$ context is shown in Figure 3.

Figure 3
Diagram of the DFT+$U$ action. $\mathscr{H}$ is the Hilbert space of the system (spanned by the eigenstates of the KS Hamiltonian), ${\mathscr{H}}_L$ is the localized Hilbert space, where the $d$ and $f$ states live, $|m⟩$ and $|m^{\prime }⟩$ are localized states and $U$ is the Hubbard parameter. The rectangles represent the states, and the blue filling indicates the electron occupation. This diagram allows one to explicitly see the transition from the less-than-half-occupied orbital to the more-than-half-occupied orbital by the electrons.

Since DFT essentially deals with energy minimization problems, minimizing the energy eigenvalues of ${\widehat{v}}_{XC}^U$ is obtained by moving electrons from less-than-half occupied states, $|m^{\prime }⟩$ in the figure, to more-than-half occupied states, $|m⟩$, always keeping the number of particles constant, concluded by analyzing the functional form of the ${\widehat{v}}_{XC}$ potential. This tends to split the levels, opening a gap mediated by the strength of the interaction $U$, completing the DFT+$U$ description of materials.

The position of DFT+$U$ on Jacob’s ladder is dynamic–it depends on the base functional and on how well the $U$ parameter is chosen. When $U$ is close to its physical value, the correction improves the description of correlated systems and effectively shifts the method upward, like an elevator toward heaven, not necessarily reaching it. However, poorly chosen values of $U$ can degrade the results and pull the functional downward. A similar perspective is discussed in the work of Kulik et al.[³¹]. The DFT+U computational cost is slightly higher compared with standard DFT functionals such as PBE. This modest increase is negligible when compared to the substantial computational demand of hybrid functionals. Therefore, if one aims to describe strongly correlated systems accurately, DFT+$U$ offers a great cost–benefit ratio relative to standard DFT: the computational cost remains nearly the same, while the physical description is significantly improved [³²].

3. Computational Details

All calculations were performed using the projector-augmented wave (PAW) method as implemented in Quantum ESPRESSO [³³, ³⁴, ³⁵]. The calculations were performed using the spin-polarized formulation of DFT. A plane-wave cutoff energy of 65 Ry was used, and a Monkhorst-Pack grid with a $5\times 5\times 5$k-mesh was employed to sample the Brillouin zone. Ionic relaxation was considered converged when the forces between atoms were smaller than ${10}^{-7}\text{Ry/}{a}_0$, and the electronic structure was considered converged once the energy difference between electronic steps were smaller than ${10}^{-8}\text{Ry/}{a}_0$. We used the PBE functional, which is a particular implementation of GGA [³⁶], with a Hubbard $U$ correction equal to $5.4\text{eV}$ applied to the Fe $d$ electrons.

The $U$ value was obtained using the linear response approach. The method proposed by Cococcioni et al. [²⁶] is based on the observation that, in an exact description, the total energy as a function of the occupation number should be piecewise linear between integer occupations. Deviations from this linear behavior arise from self-interaction errors (SIE) inherent to approximate DFT functionals. The key idea is to determine the value of $U$ that minimizes the SIE, i.e., the value that restores piecewise linearity in the energy versus occupation curve. In this framework, the unphysical curvature caused by the SIE is effectively canceled when the numerical value of $U$ is set equal to the curvature of the LDA or GGA energy profile, which is precisely the correction provided by the method [²⁶].

4. Using DFT+U - FeO

In this session, we apply the DFT+$U$ (PBE+$U$ from now on) to the iron (II) oxide, FeO, which presents a cubic Rock Salt crystal structure with an antiferromagnetic ordering; the Fe planes alternate the magnetization along the [111] direction as illustrated in Figure 4. The strong coupling between O $2p$ and Fe $3d$ orbitals leads to strongly localized and correlated electrons, which calls for approaches beyond LDA and GGA. Therefore, FeO offers a playground to apply PBE+$U$ to model strongly correlated systems. For the Hubbard correction parameter, we obtained $5.4\text{eV}$ for the Fe atoms, which agrees with previously calculated values in the literature [⁸, ³⁷, ³⁸, ³⁹]. Since the ${\widehat{v}}_{XC}^U$ only acts on localized states, its value is applied only to Fe $d$ electrons. This step was crucial since the Hubbard $U$ is an external system-dependent parameter.

Figure 4
FeO crystal structure. On the left, the Rock Salt FCC structure, and on the right, the ABC stacked structure in the $\left[111\right]$ direction.

Furthermore, the orbital occupation symmetry is broken [²⁶]. By Crystal Field Theory (CFT), the octahedral FeO geometry splits the five-fold degenerate d states into two sets with different energies [⁴⁰]. One of these sets is triply degenerate ($t_{2g}$), and the other is doubly degenerate ($e_g$). To describe the FeO system as a Mott insulator, the $a_{1g}$ level needs to be occupied, or the only level which has a different occupation eigenvalue in the set of $t_{2g}$ states [²⁶]. This nudge on occupation numbers is sufficient to take the system off a metallic local minimum on the energy landscape. Other forms of symmetry breaking may exist, allowing the system to transition to states of even lower total energy.

Compared to the PBE description, the PBE+$U$ model successfully reproduces the FeO GS magnetic properties. We compare three different magnetic configurations for both functionals, ferromagnetic (FM), AFM, and non-magnetic (NM), as exhibited in Figure 5. The energy difference between FM and AFM configurations, predicted by PBE is $\Delta E=0.019\text{eV/atom}$. The AFM phase is separated from the NM phase by $\Delta E=0.57\text{eV/atom}$.

Figure 5
Comparison between different magnetic phases calculated with (a) PBE and (b) PBE+$U$. Each point on each curve represents a self-consistent electronic calculation. The inset shows the quasi-degenerate region between AFM and FM.

Clearly, the energetically most favorable configuration in the PBE+$U$ formalism is the AFM, which indicates that our model indeed gives the proper experimental GS magnetic configuration [⁴¹]. Furthermore, the energy difference from the AFM to the FM configuration is $\Delta E=0.32\text{eV/atom}$. This difference exceeds $k_{B}T$ at room temperature ($k_{B}T=0.026eV$), indicating that the model predicts an AFM GS, and thermal fluctuations at ambient temperature should not lead to a magnetic phase transition.

The pivotal importance of including the $U$ correction to describe FeO properties can be highlighted by comparing the results obtained with and without the $U$ correction. Table 1 summarizes lattice parameter ($a$), Fe magnetic moments ($M_{Fe}$), and band-gap values in both cases. While the lattice parameter became overestimated with the addition of the Hubbard parameter $U$, which is expected for GGA calculations, Fe magnetic moments became much more accurate, illustrating the importance of the $+U$ correction. FeO is also predicted as an insulator with the correct band-gap value, in contrast with the metallic GS predicted by PBE.

Table 1
Comparison between the PBE and PBE+$U$ methods calculated in this study, and experimental data available in the literature. One can note the better description of magnetic moments and the good agreement of the experimental and calculated band gaps when using the +$U$ correction.

Method/source	$a$ (Å)	$M_{Fe}$ (${\mu }_B$)	band gap (eV)
PBE	4.27	3.27	metallic
PBE+U	4.43	3.5	2.54
Exp. [42, 43]	4.33	3.32–4.2	2.4

Figure 6 shows the projected DOS for PBE and PBE+$U$ calculations, respectively. The lattice parameters used are consistent with Table 1, i.e.$4.27$ Å for PBE and $4.43$ Å for PBE+$U$. The DOS is projected onto the $d$ spherical harmonics of Fe atoms and the $p$ states of O atoms. One can also see the direct ${\widehat{V}}_{XC}^U$ action in Figure 6 as the gap opened, if compared with the PBE DOS [⁸]. Therefore, PBE alone is not able to capture a strong correlation effect due to the $d$ states in Fe atoms [⁴⁴]. When adequately treated with an on-site interaction $U$ on the $d$ electrons, the experimental GS is captured.

Figure 6
FeO projected DOS calculated with (a) PBE and (b) PBE+U (b). One can see the gap opening once using the Hubbard $U$ correction.

The FeO ground state corresponds to the global minimum of the DFT+$U$ energy functional, which can be reached when orbital occupational symmetry breaking is allowed. In contrast, calculations performed without breaking this symmetry often converge to a higher-energy local minimum, leading to a metallic solution. However, when orbital occupational symmetry breaking is allowed–without changing the underlying Hamiltonian–the system is able to explore a broader region of the functional landscape and reach a lower-energy configuration, which better describes the true ground state of FeO. As discussed by Matteo Cococcioni and Stefano de Gironcoli [²⁶], introducing symmetry breaking based on physical insight can guide the system away from metastable solutions and toward deeper minima of the same energy functional, potentially yielding more accurate representations of correlated materials like FeO.

We can also compare these results with hybrid functional results reported by Meng et al.[⁸]. In their work, for FeO, using 0.25, 0.15, and 0.1 fractions of exact exchange, the calculated bandgaps were 2.32, 0.38, and 0.14 eV, respectively. Hybrid functionals, however, come with a much higher computational cost due to their non-local exact Fock exchange integral.

5. Conclusions

In this study, we explored the treatment of FeO, a Mott insulator, using DFT and the DFT+$U$ approaches. Our results indicate that the PBE+$U$ method significantly improves the accuracy of magnetic moments on Fe atoms and correctly predicts the insulating ground state of FeO, both qualitatively and quantitatively. While standard PBE yielded slightly better results for the lattice parameter, this does not suggest any substantial advantage of PBE for such systems. The improved description of electronic and magnetic properties by PBE+$U$ underscores its necessity for accurately modeling Mott insulators.

Nevertheless, the use of DFT+$U$ must be approached with great caution. As discussed, an inappropriate choice of the $U$ parameter can actually degrade the quality of the calculation, even compared to standard PBE results, depending on the selected value. Therefore, while DFT+$U$ is often a good choice for describing the electronic and magnetic properties of strongly correlated materials, its reliability strongly depends on the chosen $U$ and, consequently, on the methodology used to determine this value.

Through this study, we have shown that the Hubbard $U$ correction, when applied within the framework of density functional theory, is an important method, besides others like DFT+DMFT or GW, that captures the essential physics of strongly correlated materials. Both DFT+DMFT and GW rely on post-processing methods that involve constructing Green functions from DFT output data, and are considerably more expensive in terms of computational cost [⁴⁵, ⁴⁶]. This work also serves as an educational resource, highlighting the importance of DFT+$U$ in accurately simulating systems similar to FeO and emphasizing the limitations of conventional DFT approaches for strongly correlated systems.

Acknowledgments

The authors acknowledge the computational resources provided by the Multiuser Computer Center (CCM) of the Federal University of ABC (UFABC). AC thanks the financial support from the Brazilian National Council for Scientific and Technological Development-CNPq (Grants 408397/2022–5 and 302584/2023-4) and the Sao Paulo Research Foundation – FAPESP (Grant 2022/11894-1). DR thanks CNPq for the scholarship. LB thanks FAPESP for the financial support (Grant 2018/25718-5).

A. Hohenberg and Kohn theorems and KS approach

DFT provides an efficient approach to calculate the total GS energy as stated in equation 1. It is possible to calculate all GS properties of any system through DFT formalism[⁷, ⁵, ¹³].

The success of DFT is founded on the Hohenberg-Kohn (HK) theorems. The first theorem guarantees the existence of an unambiguous relationship between the external potential “felt” by electrons and the electronic density in the GS. It can be stated as:

• $\mathit{i}\mathbf{.}$ In the GS, the external potential $v_{ext}\left(\mathit{r}\right)$ is uniquely determined by the GS electronic density $n_0\left(\mathit{r}\right)$.

Thereby, in DFT, the GS electronic density $n_0\left(\mathit{r}\right)$ goes from an observable to a key object, and all properties are fully described once you have determined $n_0\left(\mathit{r}\right)$. Albeit there is a connection between $n_0\left(\mathit{r}\right)$ and $v_{ext}\left(\mathit{r}\right)$, we still need to find how this relationship can be expressed and determined quantitatively. The second HK theorem gives a hint on how we can connect $n_0\left(\mathit{r}\right)$ and $v_{ext}\left(\mathit{r}\right)$.

• $\mathit{i}\mathit{i}\mathbf{.}$ The GS energy is the global minimal of the functional $E\left[n\left(\mathit{r}\right)\right]$, and the density $n_0\left(\mathit{r}\right)$ that minimizes $E$ is the GS electronic density.

From a classical electrodynamics perspective, one can calculate the potential energy ($E_v$) of a point charge $q$ in an electrostatic potential $v\left(\mathit{r}\right)$ by [¹⁴]:

(A.1) $$E_v=q\cdot v\left(\mathit{r}\right),$$

analogously, the potential energy of a charge density $n\left(\mathit{r}\right)$ due to an external potential $v_{ext}\left(\mathit{r}\right)$ is given by:

(A.2) $$E_{v_{ext}}\left[n\left(\mathit{r}\right)\right]=\int d^{3}rn\left(\mathit{r}\right)v_{ext}\left(\mathit{r}\right),$$

which is an energy functional of the electronic density, but it exclusively computes the energy associated with the external potential. We need to construct an energy functional that determines the total energy to use theorem $\mathit{i}\mathit{i}$. Within this framework, Kohn and Sham proposed a functional (Kohn-Sham Ansatz), which consists of introducing a non-interacting auxiliary system with the same GS electron density of the original interacting system [²⁰, ⁵], but not necessarily the same wavefunction. In this way, it is possible to separate the total energy into a universal non-interacting term $F\left[n\left(\mathit{r}\right)\right]$ plus a coupling term between the auxiliary system and the specific surroundings, $v_{ext}\left(\mathit{r}\right)$. This approach is now known as KS-DFT [⁴⁷].

The total energy functional of this system is expressed as [¹³, ⁵, ²⁰]:

(A.3) $$E\left[n\left(\mathit{r}\right)\right]=F\left[n\left(\mathit{r}\right)\right]+\int d^{3}rn\left(\mathit{r}\right)v_{ext}\left(\mathit{r}\right),$$

where $F$ represents the universal auxiliary system energy functional, and the integral is the coupling term. Expanding $F\left[n\left(\mathit{r}\right)\right]$:

(A.4) $$F\left[n\left(\mathit{r}\right)\right]=T\left[n\left(\mathit{r}\right)\right]+E_H\left[n\left(\mathit{r}\right)\right]+E_{XC}\left[n\left(\mathit{r}\right)\right],$$

$T$ is the independent-particle kinetic energy, $E_H$ is the Hartree interaction energy, and $E_{XC}$ is the Exchange and Correlation energy. KS-DFT framework is an exact approach [⁵, ¹³]. As discussed in section 2, the exact mathematical form of $E_{XC}$ is unknown, and its approximations are what make KS-DFT an effective many-body mean field theory [²⁰, ³⁶, ⁴⁸].

After building a useful functional and applying the second theorem to it, we end up with a set of coupled equations known as the Kohn-Sham (KS) equations [¹³, ²⁰, ⁵]. They are explicitly written below:

(A.5) $$\left[-\frac{{\hslash }^2}{2m}{\nabla }^2+v_{eff}\left(\mathit{r}\right)\right]{\varphi }_i\left(\mathit{r}\right)={ϵ}_i{\varphi }_i\left(\mathit{r}\right),$$

with

$$v_{eff}\left(\mathit{r}\right)=v_{ext}\left(\mathit{r}\right)+v_H\left(\mathit{r}\right)+v_{XC}\left(\mathit{r}\right)$$

where $v_{ext}$ is the external potential, $v_H$ is the Hartree potential and $v_{XC}=\delta E_{XC}/\delta n\left(\mathit{r}\right)$ is the XC potential [¹³].

Clearly, the KS equations are eigenvalues/eigenvectors coupled equations and solving them gives us the KS orbitals ${\varphi }_i$ and the Lagrange multipliers (eigenvalue spectrum) ${ϵ}_i$[¹³]. Their significance is not as straightforward as the Schrödinger equation. We can think of the KS orbitals as being a first-order approximation to the electron states in a real material, and the eigenvalues as being a first-order approximation to the energy spectrum of the real interacting many-body system [¹³, ⁵].

Data Availability

Data sets generated during the current study are available from the corresponding author upon reasonable request.

References

• ^[1] C. Kittel, Introduction to Solid State Physics (Wiley, Hoboken, 2004).
• ^[2] R. Waser, Nanoelectronics and Information Technology: Advanced Electronic Materials and Novel Devices (Wiley-VCH, Weinheim, 2012).
• ^[3] E.B. Isaacs and C. Wolverton, Phys. Rev. X 9, 021042 (2019).
• ^[4] K.A. Baseden and J.W. Tye, J. Chem. Educ. 91, 2116 (2014).
• ^[5] R.M. Martin, Electronic Structure: Basic Theory and Practical Methods (Cambridge University Press, Cambridge, 2004).
• ^[6] D.S. Sholl and J.A. Steckel, Density Functional Theory: A Practical Introduction (John Wiley & Sons, 2011).
• ^[7] K. Capelle, Braz. J. Phys. 36, 1318 (2006).
• ^[8] Y. Meng, X.W. Liu, C.F. Huo, W.P. Guo, D.B. Cao, Q. Peng, A. Dearden, X. Gonze, Y. Yang, J. Wang et al., J. Chem. Theory Comput. 12, 5132 (2016).
• ^[9] R.O. Jones, in: Understanding Correlated Materials with DMFT, edited by E. Pavarini, E. Koch, A. Lichtenstein and D. Vollhardt (Forschungszentrum Jülich, Jülich, 2025), v. 15.
• ^[10] K.R. Bryenton, A.A. Adeleke, S.G. Dale and E.R. Johnson, WIREs Comput. Mol. Sci. 13, e1631 (2022).
• ^[11] A.I. Liechtenstein, V.I. Anisimov and J. Zaanen, Phys. Rev. B 52, R5467 (1995).
• ^[12] S.L. Dudarev, G.A. Botton, S.Y. Savrasov, C.J. Humphreys and A.P. Sutton, Phys. Rev. B 57, 1505 (1998).
• ^[13] F. Giustino, Materials Modelling Using Density Functional Theory: Properties and Predictions (Oxford University Press, Oxford, 2014).
• ^[14] D.J. Griffiths, Introduction to Electrodynamics (Cambridge University Press, Cambridge, 2017).
• ^[15] S. Kurth and J.P. Perdew, Int. J. Quantum Chem. 77, 814 (2000).
• ^[16] S.M. Blinder and J.E. House (ed.), Mathematical Physics in Theoretical Chemistry (Elsevier Press, Amsterdã, 2019).
• ^[17] Y. Zhang, A. Ramasamy, K. Pokharel, M. Kothakonda, B. Xiao, J.W. Furness, J. Ning, R. Zhang and J. Sun, WIREs Comput. Mol. Sci. 15, e70007 (2025).
• ^[18] R. Nagai, R. Akashi and O. Sugino, Phys. Rev. Research 4, 013106 (2022).
• ^[19] J.P. Perdew and K. Schmidt, AIP Conf. Proc. 577, 1 (2001).
• ^[20] W. Kohn and L.J. Sham, Phys. Rev. 140, A1133 (1965).
• ^[21] P. Ziesche, S. Kurth and J.P. Perdew, Comput. Mater. Sci. 11, 122 (1998).
• ^[22] M. Liu, A. Gopakumar, V.I. Hegde, J. He and C. Wolverton, Phys. Rev. Mater. 8, 043803 (2024).
• ^[23] P. Verma and D.G. Truhlar, Theor. Chem. Acc. 135, 8 (2016).
• ^[24] A. Zunger, Nat. Comput. Sci. 2, 529 (2022).
• ^[25] N.E. Kirchner-Hall, W. Zhao, Y. Xiong, I. Timrov and I. Dabo, Appl. Sci. 11, 2395 (2021).
• ^[26] M. Cococcioni and S. Gironcoli, Phys. Rev. B 71, 035105 (2005).
• ^[27] D.P. Arovas, E. Berg, S.A. Kivelson and S. Raghu, Annu. Rev. Condens. Matter Phys. 13, 239 (2022).
• ^[28] R. Shankar, Quantum Field Theory and Condensed Matter: An Introduction (Cambridge University Press, Cambridge, 2017).
• ^[29] C.L. Cleveland and R.A. Medina, Am. J. Phys. 44, 44 (1976).
• ^[30] B. Himmetoglu, A. Floris, S. de Gironcoli and M. Cococcioni, Int. J. Quantum Chem. 114, 14 (2013).
• ^[31] H.J. Kulik, J. Chem. Phys. 142, 240901 (2015).
• ^[32] M. Capdevila-Cortada, Z. Łodziana and N. López, ACS Catal. 6, 8370 (2016).
• ^[33] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G.L. Chiarotti, M. Cococcioni, I. Dabo et al., J. Phys.: Condens. Matter 21, 395502 (2009).
• ^[34] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M.B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni et al., J. Phys.: Condens. Matter 29, 465901 (2017).
• ^[35] P. Giannozzi, O. Baseggio, P. Bonfá, D. Brunato, R.Car, I. Carnimeo, C. Cavazzoni, S. Gironcoli, P. Delugas, F.F. Ruffino et al., J. Chem. Phys. 152, 105 (2020).
• ^[36] J.P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
• ^[37] W.B. Zhang, Y.H. Deng, Y.L. Hu, K.L. Han and B.Y. Tang, Solid State Commun. 142, 6 (2007).
• ^[38] J. Zaanen and G.A. Sawatzky, J. Solid State Chem. 88, 8 (1990).
• ^[39] C. Rödl, F. Fuchs, J. Furthmüller and F. Bechstedt, Phys. Rev. B 79, 235114 (2009).
• ^[40] M.S. Dresselhaus, G. Dresselhaus and A. Jorio, Group Theory: Application to the Physics of Condensed Matter (Springer, Berlin, 2007).
• ^[41] P. Zhang, R.E. Cohen and K. Haule, J. Phys.: Conf. Ser. 827, 012006 (2017).
• ^[42] W.L. Roth, Phys. Rev. 110, 1333 (1958).
• ^[43] P.D. Battle and A.K. Cheetham, J. Phys. C: Solid State Phys. 12, 337 (1979).
• ^[44] N. Yahi, Y. Azzaz, M. Ameri, M. Benouis, D. Bensaid, O. Arbouche, M. Yamani and N. Moulay, Phys. Solid State 62, 472 (2020).
• ^[45] G. Kotliar, S.Y. Savrasov, K. Haule, V.S. Oudovenko, O. Parcollet and C.A. Marianetti, Rev. Mod. Phys. 78, 865 (2006).
• ^[46] D. Golze, M. Dvorak and P. Rinke, Front. Chem. 7, 1 (2019).
• ^[47] P. Verma and D.G. Truhlar, Trends Chem. 2, 302 (2020).
• ^[48] J. Sun, A. Ruzsinszky and J.P. Perdew, Phys. Rev. Lett. 115, 036402 (2015).

Edited by

• Editor-in-Chief:
  Marcello Ferreira https://orcid.org/0000-0003-4945-3169

Publication Dates

• Publication in this collection
  06 Mar 2026
• Date of issue
  2026

History

• Received
  31 July 2025
• Reviewed
  02 Nov 2025
• Accepted
  03 Jan 2026