The effect of double counting, spin density, and Hund interaction in the different DFT+U functionals.

A systematic comparative study has been performed to better understand DFT+U (density functional theory + U) method. We examine the effect of choosing different double counting and exchange-correlation functionals. The calculated energy distribution and the Hund-J dependence of potential profile for representative configurations clearly show the different behaviors of each DFT+U formalism. In particular, adopting spin-dependent exchange-correlation functionals likely leads to undesirable magnetic solution. Our analyses are further highlighted by real material examples ranging from insulating oxides (MnO and NiO) to metallic magnetic systems (SrRuO3 and BaFe2As2). The current work sheds new light on understanding DFT+U and provides a guideline to use the related methods.

Introduction

Ab initio description of strongly correlated materials has been a challenge in condensed matter physics and materials science. A promising and widely-used scheme is to combine local density approximation (LDA) with Hubbard-type model Hamiltonian approach within density functional theory (DFT) framework[1]. One of the earliest attempts of this kind is DFT+U[2-4] which is now established as a standard approach. However, the calculation results of this type of methods strongly depends on the choice of double-counting energy functionals (which remove the conceptually equivalent contribution already present in LDA or GGA (generalized gradient approximation)) as well as interaction parameters (such as on-site Coulomb repulsion U and Hund interaction J). This feature severely limits the predictive power of DFT+U and its cousins such as DFT + DMFT (dynamical mean-field theory).

There have been many attempts to establish a proper double-counting scheme[5-14]. The difficulty lies in the nonlinear dependence of exchange-correlation (XC) functionals on the charge and/or spin density. It is therefore non-trivial to extract the precise portion of LDA/GGA XC energy for the correlated subspace. Ever since its first invention of DFT+U method, several phenomenological recipes have been suggested among which most widely used are so-called FLL (fully localized limit)[3,5-7] and AMF (around mean-field)[2,6]. Even though these double counting implementations have been extensively exploited, a comprehensive understanding of their working principles has not been reached. It is still unclear how and how much all these different formalisms give different results and predictions. In spite of previous analyses including some recent case studies of transition-metal systems within FLL[6-8,15-18], many functionals seem to be used often at random choice and without a proper guiding principle. As a result, it remains difficult to compare the results or predictions obtained by different DFT+U formalisms.

In this paper, we perform a comparative study of representative DFT+U functionals including FLL and AMF double countings. The effect of XC functional choice is also examined. To understand the detailed working principles of each DFT+U formalism, we first examine the simplified model systems in terms of their energetics and potentials. Special attention has been paid to the J dependence which has rarely been addressed before. Our analysis clearly shows the different behaviors of DFT+U functionals and their origins. In particular, when spin-polarized version of LDA or GGA is adopted, it can likely produce the undesirable effects. The characteristic features are further highlighted with real material examples covering strongly correlated insulating oxides (MnO and NiO) and metallic magnetic systems (SrRuO3 and BaFe2As2). Our work sheds new light on understanding DFT+U formalism and related methodology, thereby providing an useful guideline for its applications.

Formalism

In this section for the completeness and clarity of our presentation and notation, we briefly summarize DFT+U formalisms within non-collinear density functional scheme. Simplification to collinear case is straightforward. We use ‘CDFT’ and ‘SDFT’ in referring to the density functional scheme within charge (spin-unpolarized) density and spin density LDA/GGA, respectively. Hence ‘CDFT+U’ refers to LDA+U or GGA + U, and ‘SDFT+U’ to LSDA + U (local spin density approximation + U) or SGGA + U (spin-polarized GGA + U). Also, we use terms “cFLL”/“cAMF” to denote CDFT+U with FLL/AMF double counting and “sFLL”/“sAMF” to their SDFT+U versions.

DFT+U energy functionals

DFT+U total energy correction to CDFT or SDFT can be written as[4]:where and refers to the interaction energy within d- or f-shells and the double counting term, respectively for a particular atom s. From now on, we omit atom index s for simplicity. In the present study, E refers to either (FLL) or (AMF) depending on the choice of double counting term.

The FLL form of Eint reads[3,19]:where are the elements of on-site density matrix (DM) n for orbitals {m} and spins σ, σ′ ( or )[20,21]. The matrix elements of on-site Coulomb interaction can be expressed by[3,22]:where α and F refers to Racah-Wigner numbers and Slater integrals, respectively[3,22], and S is a transformation matrix from spherical harmonics to the predefined local basis sets. We follow the conventional expression of U = F0, J = (F2 + F4)/14, and F4/F2 = 0.625 for d-orbitals. The effect of using different ratio between F4 and F2 is found to be negligible (see Supplementary Information).

Expressing Edc has long been an important issue and still remains as an open problem[11,12]. Note that Edc itself should depend on the given XC energy functional. The FLL double counting based on CDFT+U (or cFLL) can be written as[5,7]:where N = Tr[n] within the correlated subspace. For SDFT+U (or sFLL), effect of spin-polarized XC energy should also be taken into account[3,6,15]:where the magentization and is Pauli matrices[15]. Note that the difference is the third term of Eq. (5). This fomulation of Eq. (5) has been widely used.

In AMF formalism[2,6,15], the energy correction is given by the fluctuation with respect to the average occupation of the correlated orbitals[2]:where are the elements of the redefined DM . In CDFT+U (or cAMF)[2],where l denotes the angular momentum quantum number for the correlated subspace (e.g., l = 2 for d-shells) and I is the identity matrix. In SDFT+U (or sAMF)[6,15],

DFT+U potentials

The matrix elements of orbital dependent potentials are given by . For FLL, the interaction potential for spin diagonal and off-diagonal part is given respectively by[3,19],and

Here, denotes the opposite spin to σ. Within CDFT, the double counting potential is[5,7]:and

Note that the off-diagonal components vanish and thus are spin independent. It is in a sharp contrast to the case of SDFT. In SDFT+U[3,6,15]:where (taking trace over orbitals m).

In AMF, the potential is given by taking derivative of Eq. (6) with respect to density fluctuation [2,6,15]:where refers to Eqs (7) and (8) for CDFT+U (or cAMF) and SDFT+U (or sAMF), respectively.

Analysis of Model Systems

To get a systematic understanding of how each DFT+U functional works, we analyze model systems in this section. We investigate the behaviors of energy functionals and potentials as a function of key parameters, which provides useful insight into their differences.

Energetics

In general, DFT+U DM is not necessarily diagonal[3]. As it can always be diagonalized, however, we below assume the diagonalized DM without loss of generality.

Total energy corrections by DFT+U in the case of collinear spins are now reduced to (The double counting form of AMF formalism is simplified by using the sum rule of Coulomb interaction matrix; and )[6,23]):wherefrom in Eq. (2). The fourth terms in Eqs (18 and 20) are responsible for the effective exchange interaction of SDFT (i.e., LSDA/SGGA). To represent the precise amount of this energy is a non-trivial task. Here we follow the conventional way of using Stoner parameter I with which SDFT contribution to the energy gain via spin polarization is represented by ΔESDFT = −IM2/4[2,24,25]. Note that in sFLL, this contribution is cancelled out when J = I (see Eq. (18)).

Now let us see how these functionals work in different conditions. Before taking real material examples in the next section, we consider some idealized model systems. With a fixed value of U = 5 eV, the energy distributions of d-shell electronic configurations are presented in Fig. 1 (see also Fig. 3 of ref.[23]). We use both J and I as control parameters. Here, all possible configurations of integer occupancy for a given electron number N are considered (e.g., 10C4 = 210 configurations for N = 4). We present the energy from DFT+U and XC functional contributions, which is defined as for sFLL (sAMF) and for cFLL (cAMF).

Figure 1

The energy distribution calculated by four different functionals; (a–d) sFLL, (e) cFLL, (f–i) sAMF, and (j) cAMF. E is defined as for sFLL (sAMF) and for cFLL (cAMF). All possible configurations with integer occupancy for given N have been considered. The value of U is fixed to 5 eV.

Figure 1(a) shows the result of J = 0 which can represent so-called ‘simplified rotationally invariant’ formalism by Dudarev et al.[26]. Note that the configurations with the same N are degenerate within EsFLL and this degeneracy is lifted by SDFT energy of −IM2/4. Therefore, the largest possible M configuration is always favored energetically. By comparing Fig. 1(a) with (d), one can clearly notice the role of J; lifting degeneracy within the same N-M configurations[23].

If the energy contribution from SDFT is negligible (i.e., I = 0 in ΔESDFT; Fig. 1(b)), the smaller M configurations are favored. Only when it becomes significant (Fig. 1(d)), the larger M states are stabilized and the Hund’s first rule is satisfied. While sFLL has been considered to be appropriate for high spin systems[23], this behavior is mainly attributed to SDFT exchange rather than to DFT+U correction, , as clearly seen by comparing Fig. 1(b,d). In sFLL, the low spin or nonmagnetic solution is favored as far as J is significantly larger than I; see Fig. 1(c).

In cFLL, the spin state is controlled solely by the term . Note that Fig. 1(e) is quite similar with Fig. 1(d). If I = J in sFLL, the third term in Eq. (18) cancels ΔESDFT contribution and sFLL becomes equivalent to cFLL. If the exchange contribution implicit in SDFT is larger than J (i.e., I > J), sFLL favors the larger M state more than cFLL (compare Fig. 1(d,e)).

The estimation of the intrinsic exchange in SDFT is not trivial and in general material dependent. Recent works reported that it is about ~1.0–1.5 eV for 3d transition metal systems such as nickelates, SrMnO3, SrVO3, and bcc Fe, which can be regarded as large[17,18]. As shown in Fig. 1(b–d), the exchange contribution from SDFT plays a major role in determining the moment formation, and therefore sFLL can prefer the unphysically large moment solutions. Further, SGGA has in general the stronger tendency toward the magnetic solution than LSDA[27], which is another source of ambiguity. It is certainly a drawback of SDFT+U especially for predicting material property.

In the case of AMF, the difference between CDFT and SDFT is more dramatic; see Fig. 1(f–j). As studied by Ylvisaker et al.[23], sAMF favors the low spin state and requires quite large value of I to recover Hund’s first rule. As shown in Fig. 1(i), sAMF still favors the lowest moment solution even for I = 1 eV, which is in a sharp contrast to cAMF favoring the moment formation as in cFLL (Fig. 1(e,j)). It is attributed to the fourth term of Eq. (20) which penalizes the larger moment formation. For example, with U = 5 eV, . Thus I should be greater than 1 + 4J/5 for exchange energy gain by SDFT. This feature can cause some practical problems in using AMF functionals.

J-dependence of potentials

To understand the effect of J on the moment formation and spectral property, here we further analyze DFT+U potentials. The J-only contribution to DFT+U potentials (separated from U contributions) for an orbital m and spin σ can be expressed as (assuming the diagonalized DM):where is obtained from Eq. (9) by taking non-monopole terms in Coulomb interaction matrix elements,and 〈m1, m2|V|m1, m2〉 is defined asIn Eqs (22–25), the second terms are double counting contributions.

One can clearly notice that sFLL and sAMF potentials have the spin-dependent double counting which causes the additional up/down spin potential difference. Namely, the spin-splitting is affected by double counting terms. For and , on the other hand, the spin-splitting is only controlled by interaction potential, .

In Fig. 2, the calculated J-induced spin-splittings for the model systems are presented (see also Table 1 for the list of configurations). The potential difference, for a given orbital α, can be estimated in the unit of J through Eqs (22–26). Noticeable is the same behavior of cFLL and cAMF, in which is quite substantial and always positive, favoring the moment (M) formation. This feature is attributed to the spin potential in Eqs (22) and (24) where the spin-splitting is only controlled by due to the exact cancellation of up- and down-spin double counting potentials. Thus, it is not specific to a particular form of double counting scheme. Note that the effect of J in CDFT+U (cFLL and cAMF) is consistent with what is expected from Hartree-Fock approximation.

Figure 2

The calculated J-induced spin splitting for the configurations and a given orbital α defined in Table 1. The results are presented in the unit of J.

Table 1

The electronic configurations considered in Fig. 2.

Configuration	M	N	Occupation	α
	1	1	|00100; 00000〉	d xy
1′	1	1	|00131313;00000〉	dxy, dzx, dyz
2	1	5	|00111; 00011〉	d xy
2′	1	5	|00111;00232323〉	dxy, dzx, dyz
3	1	7	|10111; 00111〉	d_z^2
3′	1	7	|1212111;00111〉	d_z^2,d_x^2−y^2
4	1	9	|11111; 01111〉	d_z^2
4′	1	9	|11111;1212111〉	d_z^2,d_x^2−y^2
5	2	2	|00110; 00000〉	dxy, dzx
5′	2	2	|00232323;00000〉	dxy, dzx, dyz
6	2	4	|00111; 00100〉	dzx, dyz
6′	2	4	|00111;00131313〉	dxy, dzx, dyz
7	2	8	|11111; 00111〉	d_z^2,d_x^2−y^2
8	3	3	|00111; 00000〉	dxy, dzx, dyz
9	3	7	|11111; 00011〉	d xy
9′	3	7	|11111;00232323〉	dxy, dzx, dyz
10	4	4	|01111; 00000〉	d xy
10′	4	4	|1212111;00000〉	dxy, dzx, dyz
11	4	6	|11111; 00100〉	dzx, dyz
11′	4	6	|11111;00131313〉	dxy, dzx, dyz
12	5	5	|11111; 00000〉	dxy, dzx, dyz

In the fourth column, d-shell occupations are presented in a form of where denotes the number of electrons occupied in the m orbital with spin σ. The primed configurations refer to the fractional occupations. The magnetic moment M (in the unit of μ) and the number of electrons N are given in the second and third column, respectively. The α are chosen to represent the lowest unoccupied or partially occupied down-spin orbitals assuming octahedral environment.

Very different features are found in sFLL where the sign of depends on the configuration. In particular, for configurations of M ≥ 3 (i.e., configuration 8–12), sFLL suppresses the spin-splittings, which is the case of SrMnO3 reported by Chen et al.[18] (see configuration 8). The trend of suppressing spin-splitting is most pronounced at half-filling (configuration 12), e.g., MnO. Further, it is important to note that the negative spin-splitting is not a general feature of sFLL double counting contrary to what is speculated by ref.[18]. See the positive configurations in Fig. 2. Our result clearly shows that both sFLL and sAMF can produce the positive spin-splitting potential.

We note that SDFT+U (sFLL and sAMF) behaves in a counter-intuitive way from the point of view of Hartree-Fock picture. It is because the spin-dependent double countings do not in general cancel out the exchange interaction from SDFT. To recover the Hartree-Fock behavior, it is desirable to use CDFT+U.

Application to Real Materials

Calculation detail

All calculations were performed using our new implementation of DFT+U into OpenMX software package (www.openmx-square.org), which is based on the nonorthogonal LCPAO (linear combination of localized pseudoatomic orbitals) formalism[28-30]. We adopted Troullier-Martins type norm-conserving pseudopotentials[31] with partial core correction. We used 9 × 9 × 9, 12 × 12 × 12 (8 × 8 × 6), and 14 × 14 × 7 k-points for rocksalt MnO and NiO, cubic (orthorhombic Pbnm) SrRuO3, and BaFe2As2 in the first Brillouin zone, respectively, and the energy cutoff of 500 Ry for numerical integrations in real space grid. The localized orbitals were generated with radial cutoff of 6.0 (Mn, Ni, and Fe) and 7.0 (Ru) a.u.[28,29]. Experimental lattice parameters were used for all materials. For the XC functional, L(S)DA[32] parameterized by Perdew and Zunger[33] was used. Unless otherwise specified, we adopted ‘dual’ projector[34] for on-site DM. For more discussion on local projectors in LCPAO scheme, see ref.[34].

MnO and NiO

Now we consider real materials. The first examples are MnO and NiO, corresponding to the configuration 12 and 7 in Fig. 2, respectively (see also Table 1). Although these two prototype correlated insulators have been extensively studied by using DFT+U, the systematic J-dependence of the electronic and magnetic property has rarely been addressed.

In Fig. 3, the calculated spin-splittings and magnetic moments by four different DFT+U formalisms (namely, cFLL, sFLL, cAMF, and sAMF) are compared as a function of J. First of all, we note that the calculated is consistent with our analyses presented in Fig. 2. In MnO, the splitting is rapidly increased in cFLL and cAMF as J increases, which is consistent with the positive value of in Fig. 2. On the other hand, it is gradually reduced in sFLL as a function of J, being consistent with the small and negative in Fig. 2. The results of NiO are also very well compared with the configuration 7 in Fig. 2.

Figure 3

The J dependence of spin-splitting () and M in the ground states of (a,b) MnO, (c,d) NiO, and (e,f) cubic SrRuO3. U = 3, 5, and 2 eV for MnO, NiO and cubic SrRuO3, respectively. Left and Right panels present the spin-splitting and M, respectively. In MnO and NiO, the energy level corresponding to orbital α and spin σ is quantified by its center of mass position of DOS; , where is DOS for given α and σ at the energy E. For SrRuO3, due to the strong Ru d - O p hybridization, the spin-splitting was estimated by the up- and down-spin difference of DOS peak position.

It is noted that sAMF predicts the entirely wrong magnetic ground state,  μ/Mn (see green lines in Fig. 3(a,b)). This low spin configuration is no longer represented by configuration 12 in Fig. 2. This is an outstanding example to show that sAMF can unphysically favor the low spin state due to the overestimated I. In this kind of case, the use of sAMF is highly undesirable.

The high spin ground state of MnO is well reproduced by sFLL, cFLL, and cAMF in a reasonable range of J (Fig. 3(a,b)). In sFLL, this ground state configuration is obtained even at J = 0 eV due to the intrinsic exchange within SDFT (U = 0) large enough to stabilize the high spin. The calculated density of states (DOS) in Fig. 4(a,b) clearly shows the different J dependence of cFLL and sFLL functionals. While the up/down spin state split is mainly controlled by J in cFLL, it is quite significant already at small J in the case of sFLL.

Figure 4

(a–c) The calculated Mn d DOS by (a) cFLL and (b) sFLL within high-spin configurations. (c) The total energy difference ΔE for MnO as a function of J. (d–f) The calculated Ni d DOS by (d) cFLL and (e) sFLL. (f) ΔE for NiO. The upper and lower panels in the DOS plots represent up and down spin parts, respectively.

To further elucidate the difference between CDFT+U and SDFT+U, Fig. 4(c) shows the total energy difference between antiferro- and ferro-magnetic phases (ΔE = EAF − EFM) calculated by cFLL and sFLL. The J dependence of ΔE exhibits the opposite trends; as J increases, cFLL tends to less favor the AF order while sFLL more favors it. From the superexchange magnetic coupling of Jex ~ −t2/(U + 4J) (t: Mn-site effective hopping integral), the behavior predicted by cFLL is more reasonable than sFLL.

In NiO (Fig. 3(c,d)), the M is insensitive to J,  μ/Ni while the slight increase is observed in cFLL and cAMF following the trend of the spin-splitting (see also Fig. 4(d,e)). Here we note that in this d8 case the low and high spin configuration is irrelevant to get the ground state property. The calculated ΔE change is also quite small in sFLL (Fig. 4(f)). In cFLL, ΔE = −0.320 and −0.224 eV/f.u. at J = 0 and 1 eV, respectively, being consistent with superexchange estimation.

SrRuO3

SrRuO3 is a ferromagnetic metal with a transition temperature of T ~ 160 K[35]. DFT+U has often been used to study SrRuO3[36-39] in spite of its metallic nature[40]. Therefore it will be informative to investigate the DFT+U functional dependence in this material. The configuration 6′ in Fig. 2 and Table 1 corresponds to this case. Figure 3(e,f) shows the calculated spin-splitting and magnetic moment, respectively. They are consistent with the results of Fig. 2; namely, the slight decreasing (increasing) trend of splitting and moment in sFLL (sAMF) and the large increase in cFLL and cAMF as a function of J.

It is noted that sFLL gives the fully polarized spin moment of  μ/f.u. for both cubic and distorted orthorhombic (not shown) structures. This half-metallic phase has been reported before by using sFLL version of SDFT+U[36-38], however, it is not well supported by experiments. The result of sAMF shows the smaller spin splitting and moment than those of sFLL as also reported in ref.[38]. This behavior of sAMF and sFLL are consistent with what is observed in MnO and NiO discussed above. Namely, it is attributed to the spin-dependent double counting which depends on U as well as J in sAMF (Eq. (20)). Due to its metallic nature, the magnetism of SrRuO3 can be more sensitive to the choice of double counting.

CDFT+U (i.e., cFLL and cAMF) shows notably different behaviors. The calculated magnetic moment and splitting are gradually increased as a function of J (Fig. 3(e,f)) and the half-metallic phase is observed only for large J ( eV for cubic and 0.8 eV for orthorhombic structure). In a reasonable range of J ≃ 0.4–0.6 eV (We could not find the direct estimation of J for the d-shell. Considering the exchange coupling within t2 orbitals of Slater-Kanamori Hamiltonian, J is estimated to be ~0.4–0.6 eV[41,42].), the calculated moment is and 1.6 μ/f.u. for cubic and orthorhombic structure, respectively, in good agreement with experiments[35].

As mentioned in the previous section, the exchange contribution by SGGA is expected to be greater than the LSDA[27]. This tendency is clearly shown in Fig. 5(b). In SGGA + U, the moment size is further enhanced (M = 1.96 μ/f.u.) than LSDA + U (M = 1.67 μ/f.u.). On the other hand, in the case of CDFT+U (Fig. 5(a)), GGA + U gives basically the same result with LDA+U (M = 1.41 μ/f.u.).

Figure 5

The calculated DOS of cubic SrRuO3 by (a) cFLL and (b) sFLL. Two XC functionals for CDFT and SDFT are adopted; L(S)DA (green solid lines) and (S)GGA (magenta dotted lines). PBE (Perdew-Burke-Ernzerhof) parameterization for (S)GGA[51] is used. U = 2 and J = 0.5 eV for cFLL (a), and U = 1 and J = 0 eV for sFLL (b).

BaFe2As2

The superconducting Fe pnictides have been a subject of intensive research activities. From the viewpoint of first-principles calculations, the unusually large magnetic moment by SDFT compared to experiments is a long standing issue[43-46]. Interestingly, to reproduce experimental moments, negative U values within SDFT+U[47,48] have been adopted. As pointed out in ref.[48], however, it is hard to be justified in the physics sense. Here we note that the intrinsic exchange contribution of ~IM in SDFT can be too large as discussed in the above, and SDFT may not be the right starting point to take the correlation effects into account.

We found that CDFT+U can provide much more sensible picture for magnetism in this material. Table 2 shows the calculated magnetic moment for BaFe2As2 with cRPA (constrained random phase approximation) value of U = 2.3 eV[49]. The result of McFLL is in a fairly good agreement with experiment ( μ/Fe[50]) for J = 0.3–0.5 eV whereas MsFLL always overestimate the moments. Note that the reasonable size of M is reproduced with realistic value of U and J only within CDFT+U. As shown in Table 2, the moment is also sensitive to the way of defining local DM projector since the ‘full’ projector tends to take the smaller on-site electron occupation compared to the ‘dual’[34]. The best comparison with experiment is achieved with J = 0.3 eV for ‘dual’ and  eV for ‘full’ projector.

Table 2

Calculated magnetic moment of BaFe2As2 by cFLL (McFLL) and sFLL (MsFLL).

U [eV]	J [eV]	DM	McFLL [μB/Fe]	MsFLL [μB/Fe]
2.3	0.3	dual	0.94	2.82
		full	0.29	2.63
	0.5	dual	1.78	2.77
		full	0.75	2.59
	0.7	dual	2.34	2.73
		full	1.33	2.56

Experimental crystal sructure in spin-density-wave phase is used[50]. Two different definitions of DM are used for the comparison (namely, ‘dual’ and ‘full’).

Also noticeable is the different J dependence of moment by two functionals; McFLL (MsFLL) increases (decreases) as J increases. This feature is again consistent with the behavior discussed in the previous section. The consistent result of cFLL with experiment is impressive even though the dynamic correlation beyond DFT+U certainly plays the role in this system[40].

Summary and Conclusion

We performed a comparative analysis on DFT+U functionals employing two widely-used double counting forms and their relation to standard XC functionals. The detailed investigations on each formulation as well as the real material examples provided a clear understanding of different behaviors of DFT+U functionals. The calculated energetics and spin potentials for representative model systems clearly show the role of double counting and XC functional in determining the ground state magnetic property. Competition between the effect of J and the spin density XC energy is the key to understand the SDFT + U result. Application to real materials including MnO, NiO, SrRuO3, and BaFe2As2 further clarify the different tendency between the formalisms, supporting the analyses with model systems. As a rule of thumb, CDFT+U is suggested as the desirable choice for most purposes.