The role of spin in the calculation of Hubbard $U$ and Hund's $J$ parameters from first principles

TL;DR

This paper revises the calculation of Hubbard U and Hund's J parameters in DFT+U methods by explicitly incorporating spin effects, leading to more accurate predictions for strongly-correlated materials.

Contribution

It introduces spin-aware definitions for Hubbard and Hund's parameters within the linear response framework, aligning them with the current DFT+U functional and enabling easier computation.

Findings

Accurately reproduces experimental band gaps and magnetic moments in manganese oxide.

Shows the necessity of oxygen Hubbard corrections for bond length preservation.

Demonstrates numerical stability even for near-filled electronic subspaces.

Abstract

The density functional theory (DFT)+$U$ method is a pragmatic and effective approach for calculating the ground-state properties of strongly-correlated systems, and linear response calculations are widely used to determine the requisite Hubbard parameters from first principles. We provide a detailed treatment of spin within this linear response approach, demonstrating that the conventional Hubbard $U$ formula, unlike the conventional DFT+$U$ corrective functional, incorporates interactions that are off-diagonal in the spin indices and places greater weight on one spin channel over the other. We construct alternative definitions for Hubbard and Hund's parameters that are consistent with the contemporary DFT+$U$ functional, expanding upon the minimum-tracking linear response method. This approach allows Hund's $J$ and spin-dependent $U$ parameters to be calculated with the same ease as for the standard Hubbard $U$. Our methods accurately reproduce the experimental band gap, local magnetic moments, and the valence band edge character of manganese oxide, a canonical strongly-correlated system. We also apply our approach to a complete series of transition-metal complexes [M(H$_2$O)$_6$]$^{n+}$ (for M = Ti to Zn), showing that Hubbard corrections on oxygen atoms are necessary for preserving bond lengths, and demonstrating that our methods are numerically well-behaved even for near-filled subspaces such as in zinc. However, spectroscopic properties appear beyond the reach of the standard DFT+$U$ approach. Collectively, these results shed new light on the role of spin in the calculation of the corrective parameters $U$ and $J$, and point the way towards avenues for further development of DFT+$U$-type methods.

Figures

Tables

Table 1. Table 1: Values of U 𝑈 U and J 𝐽 J (eV) for hexahydrated transition metals and a spin-up manganese atom of MnO, calculated using the various linear response schemes introduced in subsection  II.1 . The linear response calculations for the fully-filled 3 ​ d 3 𝑑 3d subspace in Fe 3+ were poorly behaved (two different pseudopotentials were tested) and have consequently been excluded.

metal	scalar		averaged $1\times 1$		$1\times 1$				simple $2\times 2$				scaled $2\times 2$
	$U$		$U$		$U^{\uparrow }$		$U^{\downarrow }$		$U$		$J$		$U$		$J$
Ti3+	$3.88$ ±	0.00	$1.66$ ±	0.00	$1.85$ ±	0.01	$1.47$ ±	0.00	$3.90$ ±	0.01	$0.34$ ±	0.00	$3.89$ ±	0.01	$0.34$ ±	0.00
V2+	$4.00$ ±	0.00	$2.78$ ±	0.00	$3.29$ ±	0.00	$2.28$ ±	0.00	$4.07$ ±	0.01	$0.34$ ±	0.00	$4.00$ ±	0.01	$0.35$ ±	0.00
Cr3+	$3.90$ ±	0.00	$1.78$ ±	0.00	$1.86$ ±	0.00	$1.70$ ±	0.00	$4.04$ ±	0.01	$0.40$ ±	0.00	$3.90$ ±	0.01	$0.42$ ±	0.00
Cr2+	$3.20$ ±	0.00	$2.39$ ±	0.00	$2.75$ ±	0.00	$2.04$ ±	0.00	$3.34$ ±	0.01	$0.33$ ±	0.00	$3.20$ ±	0.01	$0.35$ ±	0.00
Mn3+	$5.40$ ±	0.00	$2.00$ ±	0.00	$1.51$ ±	0.00	$2.50$ ±	0.00	$5.86$ ±	0.01	$0.50$ ±	0.00	$5.40$ ±	0.01	$0.53$ ±	0.00
Mn2+	$4.36$ ±	0.00	$4.05$ ±	0.08	$4.28$ ±	0.15	$3.82$ ±	0.00	$4.90$ ±	0.06	$0.37$ ±	0.06	$4.35$ ±	0.01	$0.52$ ±	0.01
Fe3+	$5.88$ ±	0.01	—		—		$5.45$ ±	0.02	—		—		$5.92$ ±	0.02	$0.81$ ±	0.02
Fe2+	$4.58$ ±	0.00	$5.07$ ±	0.09	$6.28$ ±	0.18	$3.86$ ±	0.00	$6.06$ ±	0.09	$0.43$ ±	0.06	$4.58$ ±	0.01	$0.63$ ±	0.01
Co3+	$6.25$ ±	0.00	$1.19$ ±	0.00	$1.19$ ±	0.00	$1.19$ ±	0.00	$6.25$ ±	0.00	$0.75$ ±	0.00	$6.25$ ±	0.00	$0.75$ ±	0.00
Co2+	$4.95$ ±	0.02	$6.19$ ±	0.02	$8.17$ ±	0.03	$4.22$ ±	0.02	$7.15$ ±	0.02	$0.48$ ±	0.01	$4.96$ ±	0.02	$0.65$ ±	0.01
Ni2+	$5.26$ ±	0.00	$9.84$ ±	0.02	$15.41$ ±	0.05	$4.27$ ±	0.00	$12.35$ ±	0.03	$0.75$ ±	0.02	$5.26$ ±	0.01	$0.78$ ±	0.01
Cu2+	$4.62$ ±	0.00	$-2.54$ ±	0.03	$-9.11$ ±	0.05	$4.04$ ±	0.00	$-4.99$ ±	0.02	$0.85$ ±	0.02	$4.63$ ±	0.01	$0.90$ ±	0.01

Table 2. Table 3: Values of U 𝑈 U and J 𝐽 J (eV) for the 3 ​ d 3 𝑑 3d subspace of Zn in hexahydrated zinc, calculated using the various linear response schemes and two alternative sets of Hubbard projectors (as defined by the net charge configuration of the Zn atom in a pseudoatomic solver).

PAO charge		+0		+2
scalar	$U$	$10.05$±	0.03	$34.77$±	0.01
averaged $1\times 1$	$U$	$11.60$±	0.04	$44.64$±	0.02
$1\times 1$	$U^{\uparrow }$	$11.67$±	0.06	$44.65$±	0.03
	$U^{\downarrow }$	$11.53$±	0.06	$44.63$±	0.02
simple $2\times 2$	$U$	$10.08$±	0.03	$34.79$±	0.02
	$J$	$1.75$±	0.05	$1.47$±	0.03
scaled $2\times 2$	$U$	$10.08$±	0.03	$34.79$±	0.02
	$J$	$1.75$±	0.05	$1.47$±	0.03

Table 3. Table 4: Spin flip energies (eV) for various hexahydrated transition metal systems. The quantum chemistry results are from Ref.  Yang  et al. , 2014 , and the experimental results are from Ref.  Jørgensen , 1962 (and the references therein).

metal	DFT	DFT + $U$ (+$J$)				CASSCF	CASPT2	MRCI	exp
		scalar	av $1\times 1$	$1\times 1$	scaled $2\times 2$
V2+	$1.06$	$1.10$	$1.09$	$1.28$	$1.10$	$2.01$	$1.89$	$1.98$	$1.62$
Cr3+	$1.11$	$1.04$	$1.11$	$1.33$	$1.04$	$2.41$	$2.23$	$2.35$	$2.60$
Mn2+	$2.16$	$2.41$	$2.39$	$2.41$	$2.42$	$3.42$	$2.91$	$3.25$	$2.34$
Co2+	$1.60$	$1.85$	$1.85$	$1.86$	$1.85$	$1.96$	$1.95$	$1.76$	$1.98$
Ni2+	$1.23$	$1.44$	$1.48$	$1.50$	$1.44$	$2.30$	$2.03$	$2.23$	$1.91$

Table 4. Table 5: Kohn-Sham transition energies (eV) for spin-conserving d 𝑑 d - d 𝑑 d excitations. In all cases, corrective terms were applied to both the metal 3 ​ d 3 𝑑 3d and oxygen 2 ​ p 2 𝑝 2p subspaces.

metal	final	DFT	DFT + $U$ (+ $J$)				CASSCF	CASPT2	MRCI	exp
	symmetry		scalar	av $1\times 1$	$1\times 1$	scaled $2\times 2$
Ti3+	1B2g	$0.27$	$3.16$	$1.48$	$1.63$	$2.79$	$0.00$
	1B3g	$0.28$	$3.25$	$1.51$	$1.66$	$2.86$	$0.00$
	1Ag	$1.94$	$3.98$	$2.81$	$2.92$	$3.76$	$1.69$	$1.71$	$1.76$	$2.16$
	1Ag	$2.38$	$4.60$	$3.33$	$3.44$	$4.34$	$1.70$	$1.72$	$1.77$	$2.52$
V2+	3B1g	$1.97$	$4.92$	$3.98$	$4.35$	$4.61$	$1.19$	$1.26$	$1.28$	$1.53$
	3B2g	$1.97$	$4.92$	$3.98$	$4.35$	$4.61$	$1.19$	$1.26$	$1.28$	$1.53$
	3B3g	$1.97$	$4.92$	$3.98$	$4.35$	$4.61$	$1.19$	$1.26$	$1.28$	$1.53$
Cr3+	3B1g	$2.24$	$3.98$	$3.15$	$3.17$	$3.56$	$1.69$	$1.77$	$1.79$	$2.16$
	3B2g	$2.24$	$3.98$	$3.15$	$3.17$	$3.56$	$1.69$	$1.77$	$1.79$	$2.16$
	3B3g	$2.24$	$3.98$	$3.15$	$3.17$	$3.56$	$1.69$	$1.77$	$1.79$	$2.16$
Cr2+	4Ag	$0.38$	$2.06$	$1.60$	$1.80$	$1.83$	$0.62$	$0.69$	$0.64$	$1.17$
	4B2g	$1.37$	$3.28$	$2.77$	$2.99$	$3.04$	$1.18$	$1.27$	$1.19$
	4B3g	$1.53$	$3.44$	$2.93$	$3.15$	$3.22$	$1.23$	$1.30$	$1.23$
	4B1g	$1.95$	$3.88$	$3.36$	$3.58$	$3.61$	$1.34$	$1.44$	$1.36$	$1.75$
Mn3+	4Ag	$0.21$	$1.28$	$0.68$	$0.62$	$0.97$	$0.69$	$0.77$	$0.72$	$1.11$
	4B2g	$0.97$	$5.08$	$3.63$	$1.86$	$4.87$	$1.72$	$1.96$	$1.78$	$2.53$
	4B3g	$2.64$	$5.28$	$3.99$	$3.41$	$4.95$	$1.76$	$1.99$	$1.82$	$2.53$
	4B1g	$3.00$	$5.72$	$4.38$	$3.88$	$5.42$	$1.91$	$2.21$	$2.00$	$2.53$
Fe2+	4B2g	$1.28$	$5.62$	$6.03$	$4.92$	$4.95$	$0.00$
	4B3g	$1.28$	$5.63$	$6.04$	$4.93$	$5.04$	$0.00$
	4Ag	$1.88$	$5.73$	$6.07$	$5.07$	$5.24$	$0.75$	$0.80$	$0.83$	$1.29$
	4Ag	$3.12$	$6.99$	$7.34$	$6.38$	$6.33$	$0.85$	$0.89$	$0.91$	$1.29$
Co2+	3B2g	$3.02$	$7.85$	$8.65$	$7.03$	$7.16$	$0.00$
	3B3g	$3.03$	$7.85$	$8.67$	$7.04$	$7.17$	$0.00$
	3B2g	$3.64$	$7.91$	$8.75$	$7.15$	$7.23$	$0.67$	$0.81$	$0.65$	$1.02$
	3B3g	$3.64$	$7.91$	$8.78$	$7.16$	$7.24$	$0.68$	$0.81$	$0.65$	$1.02$
	3B3g	$5.17$	$9.41$	$10.08$	$8.76$	$8.69$	$2.82$	$2.69$	$2.62$	$2.41$
	3B2g	$5.17$	$9.42$	$10.10$	$8.77$	$8.70$	$2.85$	$2.69$	$2.62$	$2.41$
Ni2+	2B3g	$4.30$	$8.74$	$13.39$	$7.72$	$7.91$	$0.75$
	2B1g	$4.30$	$8.74$	$13.39$	$7.72$	$7.91$	$0.76$	$0.89$	$0.85$	$1.05$
	2B2g	$4.30$	$8.74$	$13.39$	$7.72$	$7.91$	$0.76$
	2B1g	$4.30$	$8.74$	$13.39$	$7.72$	$7.91$	$1.31$	$1.48$	$1.45$	$1.67$
	2B2g	$4.31$	$8.74$	$13.39$	$7.72$	$7.91$	$1.31$
	2B3g	$4.31$	$8.74$	$13.39$	$7.72$	$7.91$	$1.31$
Cu2+	1Ag	$1.79$	$5.06$	$0.49$	$4.56$	$4.47$	$0.51$	$0.61$	$0.52$	$1.17$
	1B2g	$2.31$	$6.19$	$1.12$	$5.49$	$5.33$	$0.84$	$1.08$	$0.85$
	1B3g	$2.82$	$6.68$	$1.81$	$5.91$	$5.82$	$0.89$	$1.12$	$0.89$
	1B1g	$3.34$	$6.85$	$1.99$	$5.95$	$5.88$	$0.97$	$1.23$	$0.99$	$1.56$

Table 5. Table 7: Spin flip energies (eV) for various hexahydrated transition metal systems with U 𝑈 U (and J 𝐽 J for 2 × 2 2 2 2\times 2 ) updated following the flip. The quantum chemistry results are from Ref.  Yang  et al. , 2014 , and the experimental results are from Ref.  Jørgensen , 1962 (and the references therein). The surprisingly small (and even negative) spin-flip energies for 2 × 2 2 2 2\times 2 are strongly reminiscent of the findings of Millis and coworkers, who showed that the contemporary + J 𝐽 +J functional wrongly disfavours ferromagnetism. Chen and Millis ( 2016 )

metal	DFT	DFT + $U$ (+$J$)				CASSCF	CASPT2	MRCI	exp
		scalar	av $1\times 1$	$1\times 1$	scaled $2\times 2$
V2+	$1.06$	$1.11$	$1.26$	$1.28$	$0.80$	$2.01$	$1.89$	$1.98$	$1.62$
Cr3+	$1.11$	$0.94$	$1.32$	$1.33$	$-0.15$	$2.41$	$2.23$	$2.35$	$2.60$
Mn2+	$2.16$	$2.40$	$2.73$	$2.74$	$2.00$	$3.42$	$2.91$	$3.25$	$2.34$
Co2+	$1.60$	$1.83$	$1.44$	$1.62$	$0.72$	$1.96$	$1.95$	$1.76$	$1.98$
Ni2+	$1.23$	$1.81$	$1.03$	$1.50$	$0.41$	$2.30$	$2.03$	$2.23$	$1.91$

Equations

$$E_{\text{U}}[{\hat{n}}^{I\sigma }]=\sum _{I\sigma }\frac{U^I}{2}\mathrm{Tr}\left[{\hat{n}}^{I\sigma }(1-{\hat{n}}^{I\sigma })\right],$$

$${\hat{V}}_{\text{U}}=\sum _{I\sigma mn}{U}^I|{\phi }_m^I⟩\left(\frac{1}{2}-{n^{I\sigma m}}_n\right)⟨{\phi }^{In}|.$$

$$E_{\text{J}}[{\hat{n}}^{\sigma }]=\sum _{I\sigma }\frac{J^I}{2}\mathrm{Tr}\left[{\hat{n}}^{I\sigma }{\hat{n}}^{I-\sigma }\right],$$

$$\delta \hat{\rho }=\hat{\chi }\delta {\hat{v}}_{\mathrm{ext}}.$$

$$d{n}^I=\mathrm{Tr}\left[{\hat{P}}^I\delta \hat{\rho }\right]=\mathrm{Tr}\left[{\hat{P}}^I\hat{\chi }{\hat{P}}^J\right]d{v}_{\mathrm{ext}}^J$$

$${\chi }_{IJ}\equiv \frac{d{n}^I}{d{v}_{\mathrm{ext}}^J}=\mathrm{Tr}\left[{\hat{P}}^I\hat{\chi }{\hat{P}}^J\right].$$

$$U^I={\left({\chi }_0^{-1}-{\chi }^{-1}\right)}_{II}$$

$$\frac{d{v}_{\mathrm{KS}}^I}{d{v}_{\mathrm{ext}}^J}=\frac{d{v}_{\mathrm{ext}}^I}{d{v}_{\mathrm{ext}}^J}+\frac{d{v}_{\mathrm{Hxc}}^I}{d{v}_{\mathrm{ext}}^J}=\frac{d{v}_{\mathrm{ext}}^I}{d{v}_{\mathrm{ext}}^J}+\sum _K\frac{d{v}_{\mathrm{Hxc}}^I}{d{n}^K}\frac{d{n}^K}{d{v}_{\mathrm{ext}}^J},$$

$${\epsilon }^{-1}=\Omega +f\chi ⟹f=\left({\epsilon }^{-1}-\Omega \right){\chi }^{-1}.$$

$$U^I={\left[\left(\frac{d{v}_{\mathrm{KS}}}{d{v}_{\mathrm{ext}}}-1\right){\left(\frac{dn}{d{v}_{\mathrm{ext}}}\right)}^{-1}\right]}_{II}.$$

$${\chi }_{IJ}=$$

$${\left({\chi }_0\right)}_{IJ}=$$

$${\chi }_{IJ}^{\sigma {\sigma }^{\prime }}=\frac{d{n}^{I\sigma }}{d{v}^{J{\sigma }^{\prime }}},$$

$$\chi =\left(\begin{array}{cccc}{\chi }_{11}^{\uparrow \uparrow }& {\chi }_{11}^{\uparrow \downarrow }& {\chi }_{12}^{\uparrow \uparrow }& {\chi }_{12}^{\uparrow \downarrow }\\ {\chi }_{11}^{\downarrow \uparrow }& {\chi }_{11}^{\downarrow \downarrow }& {\chi }_{12}^{\downarrow \uparrow }& {\chi }_{12}^{\uparrow \downarrow }\\ {\chi }_{21}^{\uparrow \uparrow }& {\chi }_{21}^{\uparrow \downarrow }& {\chi }_{22}^{\uparrow \uparrow }& {\chi }_{22}^{\uparrow \downarrow }\\ {\chi }_{21}^{\downarrow \uparrow }& {\chi }_{21}^{\downarrow \downarrow }& {\chi }_{22}^{\downarrow \uparrow }& {\chi }_{22}^{\downarrow \downarrow }\end{array}\right)=\left(\begin{array}{cc}({\chi }^{\sigma {\sigma }^{\prime }}{)}_{11}& ({\chi }^{\sigma {\sigma }^{\prime }}{)}_{12}\\ ({\chi }^{\sigma {\sigma }^{\prime }}{)}_{21}& ({\chi }^{\sigma {\sigma }^{\prime }}{)}_{22}\end{array}\right).$$

$$f^{\sigma \sigma }=\frac{d{v}_{\mathrm{KS}}^{\sigma }}{d{n}^{\sigma }}-\frac{d{v}_{\mathrm{ext}}^{\sigma }}{d{n}^{\sigma }}=\frac{d{v}_{\mathrm{Hxc}}^{\sigma }}{d{n}^{\sigma }}.$$

$$U=\frac{1}{2}\left(f^{\uparrow \uparrow }+f^{\downarrow \downarrow }\right).$$

$$f^{\sigma {\sigma }^{\prime }}={\left[\left(\frac{d{v}_{\mathrm{KS}}}{d{v}_{\mathrm{ext}}}-1\right){\left(\frac{dn}{d{v}_{\mathrm{ext}}}\right)}^{-1}\right]}^{\sigma {\sigma }^{\prime }}$$

$$U=$$

$$\approx$$

$$U=\frac{1}{4}\left(f^{\uparrow \uparrow }+f^{\uparrow \downarrow }+f^{\downarrow \uparrow }+f^{\downarrow \downarrow }\right)$$

$$U=$$

$${\lambda }_U=$$

$$J=-\frac{1}{2}\frac{d{v}_{\mathrm{Hxc}}^{\uparrow }-d{v}_{\mathrm{Hxc}}^{\downarrow }}{d(n^{\uparrow }-n^{\downarrow })}$$

$$J=-\frac{1}{4}\left(f^{\uparrow \uparrow }-f^{\uparrow \downarrow }-f^{\downarrow \uparrow }+f^{\downarrow \downarrow }\right),$$

$$J=$$

$${\lambda }_J=$$

$$dn=d{n}^{\uparrow }+d{n}^{\downarrow }\approx$$

$$⟹\chi \approx$$

$$d{v}_{\mathrm{KS}}=\frac{1}{2}\left[d{v}_{\mathrm{KS}}^{\uparrow }+d{v}_{\mathrm{KS}}^{\downarrow }\right]\approx$$

$$⟹{\epsilon }^{-1}\approx$$

$$U=$$

$$=$$

$$\delta n(\mathbf{r})=\int d{\mathbf{r}}^{\prime }\chi (\mathbf{r},{\mathbf{r}}^{\prime })\delta v_{\text{ext}}({\mathbf{r}}^{\prime })$$

Full text

I Introduction

Over the past few decades, density functional theory (DFT) has played a key role in the simulation of many-body atomistic systems.Jones (2015); Jain et al. (2016) DFT makes such systems tractable via the Hohenberg-Kohn theoremsHohenberg and Kohn (1964) and the Kohn-Sham construction,Kohn and Sham (1965) but exchange and correlation must be approximated in the form of an exchange-correlation (xc)-functional.Perdew and Zunger (1981); Perdew et al. (1996)

One of the most prominent failures of many xc-functionals is that they do not properly correct for the self-interaction in the Hartree term. Self-interaction error (SIE) — or more generally “delocalization error”Cohen et al. (2008, 2012) — manifests itself as a spurious curvature in total energies with respect to total electron number, where instead there should be a derivative discontinuity at integer numbers of electrons and linear behavior at fractional numbers.Perdew et al. (1982) This failure is closely related to approximate DFT’s well-documented underestimation of the band gap.Cococcioni and de Gironcoli (2005); Kaduk et al. (2012); Mori-Sánchez et al. (2008)

While the origins of the SIE are well understood, it remains a challenge to avoid its introduction when constructing xc-functionals, even if exact exchange is incorporated.Cohen et al. (2008) A simple yet remarkably successful alternative is Hubbard-augmented DFT (LDA + $U$ or more generally DFT + $U$). In this scheme, Hubbard model terms are incorporated into the DFT framework. This approach was originally designed to capture Mott-Hubbard physics in transition-metal oxides,Anisimov et al. (1991, 1991); Dudarev et al. (1998) but it has subsequently gained a transparent interpretation as a corrective method for SIE due to the work of Cococcioni, Kulik, and co-workers.Cococcioni and de Gironcoli (2005); Kulik et al. (2006) They observed that Hubbard corrections can counteract the spurious SIE curvature — in other words, DFT + $U$ calculations may be constructed to cancel the SIE that is present (although this is not guaranteed).Zhao et al. (2016)

I.1 The DFT + U correction

In the DFT + $U$ scheme, one adds to the energy a corrective term (here we use the rotationally-invariant, simplified form),Dudarev et al. (1998); Anisimov and Gunnarsson (1991); Anisimov et al. (1991, 1993, 1997); Pickett et al. (1998) given by

[TABLE]

where the density operators $\hat{n}^{I\sigma}=\hat{P}^{I}\hat{\rho}^{\sigma}\hat{P}^{I}$ are projections of the (spin-dependent) Kohn-Sham density operator onto subspaces (indexed $I$) in which the SIE is to be addressed. The projectors ${\hat{P}^{I}=\sum_{m}|\varphi^{I}_{m}\rangle\langle\varphi^{Im}|}$ are typically constructed from atom-centered, fixed, spin-independent, localized, and orthonormal orbitals $\varphi^{I}_{m}$ (although they may be non-orthogonalO’Regan et al. (2011) and self-consistentO’Regan et al. (2010)). The $U^{I}$ are externally-defined parameters that determine the strength of the energy corrections. If they are well-chosen, the term that is quadratic in $\hat{n}$ can partially correct the spurious energy curvature arising from the SIE.Zhao et al. (2016) In the basis of localized orbitals $\psi_{m}^{I\sigma}$ that diagonalize the subspace occupancy matrices such that $\hat{n}^{I\sigma}\psi_{m}^{I\sigma}=\lambda_{m}^{I\sigma}\psi_{m}^{I\sigma}$, the Hubbard correction becomes $\sum_{I\sigma m}U^{I}\lambda_{m}^{I\sigma}(1-\lambda_{m}^{I\sigma})/2$, which penalizes non-integer occupancies of these orbitals $\psi_{m}^{I\sigma}$. The explicit correction to the total energy vanishes at integer occupancy matrix eigenvalues, where the xc-functional is assumed to be correct.

The corresponding correction to the Kohn-Sham potential is given by

[TABLE]

This is attractive or repulsive for occupancy matrix eigenvalues greater than or less than one-half, respectively. In the absence of any significant self-consistent response, this will penalize non-integer occupancies of the subspaces, opening an energy gap of order $U$ between any occupied and unoccupied Kohn-Sham orbitals with significant overlap with the Hubbard projectors.

In order to correct interactions between unlike spins, DFT + $U$ can be extended to become DFT + $U$ + $J$ .Solovyev et al. (1994); Mosey and Carter (2007); Himmetoglu et al. (2011, 2014) This involves a second correction to the total energy,

[TABLE]

where this correction is parameterized by the additional Hund’s coupling constants $J^{I}$. Additionally, the $U$ in Eqs. 1 and 2 becomes $U_{\mathrm{eff}}=U-J$.

A substantial advantage of DFT + $U$ (+ $J$) over other methods that address the SIE (for example, SIC-LSDA,Lüders et al. (2005); Däne et al. (2009); Hughes et al. (2007) Fermi orbital self-interaction correction,Pederson et al. (2014); Pederson (2015); Hahn et al. (2015) and Koopman’s compliant functionalsDabo et al. (2010); Borghi et al. (2014)) is its small computational cost: once any Hubbard parameters have been determined, the overhead for incorporating the additional potential and energy terms is insignificant compared to the cost of the DFT calculation itself.O’Regan et al. (2012)

I.2 Conventional linear response

In order to apply a Hubbard correction, one must select an appropriate value for the parameters $U^{I}$. This can be done pragmatically by picking values on empirical grounds — that is, chosen so that certain system characteristics are reproduced (for example, ionic geometries,Loschen et al. (2007); Castleton et al. (2007); Morgan and Watson (2007); Franchini et al. (2007) band gaps,Castleton et al. (2007); Morgan and Watson (2007); Rohrbach et al. (2004); Franchini et al. (2007) and formation enthalpiesFranchini et al. (2007); Wang et al. (2006); Ong et al. (2011)). While this approach has seen some success,Fennie and Rabe (2006); Finazzi et al. (2008) it does not guarantee that the chosen $U$ will correct the SIE energy curvature to the greatest extent achievable, or result in an improved description of other system properties, and is not even possible where there is a lack of reliable experimental or higher-level computational data.

An alternative approach is the linear response method developed by Cococcioni and de Gironcoli,Cococcioni and de Gironcoli (2005) which built upon the earlier linear response scheme of Pickett and co-workers,Pickett et al. (1998) and shares many aspects with the constrained LDA approach of Aryasetiawan and co-workers.Aryasetiawan et al. (2006) In this approach, DFT calculations are performed subject to a perturbing potential $\delta\hat{v}_{\mathrm{ext}}=dv^{J}_{\mathrm{ext}}\hat{P}^{J}$ confined to the $J$th Hubbard subspace, for a range of values of $dv^{J}_{\mathrm{ext}}$. The density operator’s response to these perturbations is given by the response operator $\hat{\chi}$:

[TABLE]

The occupancy of the $I$th Hubbard subspace will change by

[TABLE]

and thus we can define the projected response matrixO’Regan and Teobaldi (2016)

[TABLE]

A value for $U$ that corresponds to the screened response of the system is given by

[TABLE]

where $\chi_{0}$ is the response of the non-interacting system, which must be separately measured and removed from the Hubbard correction.Cococcioni (2002); Dederichs et al. (1984); Himmetoglu et al. (2014) There is also scope here for calculation of off-diagonal terms $V_{IJ}=\left(\chi^{-1}_{0}-\chi^{-1}\right)_{IJ}$, which gives rise to DFT + $U$ + $V$.Leiria Campo Jr and Cococcioni (2010); Kulik and Marzari (2011a)

Satisfyingly, the determination of $U$ via linear response removes any possible arbitrariness of the Hubbard correction: the $U$ parameter is a well-defined property of the system that can be unambiguously measured in theory.Kulik et al. (2006); Himmetoglu et al. (2014)

Recently the idea of calculating $U$ and $J$ to best emulate subspace-projected Kohn-Sham exact exchangeMosey and Carter (2007) has been further advanced.Agapito et al. (2015) But because we wish to cancel the systematic errors of approximate DFTPerdew et al. (1982); Cohen et al. (2012) to the extent possible using functionals of the DFT + $U$ form, we choose to instead develop the linear-response formalism that has been shown to successfully achieve this,Pickett et al. (1998); Cococcioni and de Gironcoli (2005); Kulik et al. (2006); Moynihan et al. (2017) and that does not incorporate any theory or model (e.g. Fock exchange) beyond what is already ordinarily present.

I.3 Problem and paper outline

There are some aspects of the linear response methodology that pose issues. Firstly, delocalization error is associated with fractional total charge, but the DFT + $U$ functional of Eq. 1 corrects fractional occupation for each spin channel separately. Conventional linear response, meanwhile, perturbs both spin channels simultaneously. These discrepancies in how we treat spin channels warrant investigation.

Secondly, measuring the non-interacting response $\chi_{0}$ is not straightforward. The common practice is to follow the example of Ref. Cococcioni and de Gironcoli, 2005, and calculate $\chi_{0}$ via the first iteration of the Kohn-Sham equations during a self-consistent field calculation — that is, the response is to be measured following the initial charge redistribution introduced by the perturbation but before the Kohn-Sham potential is updated. This approach is impractical to implement in codes that use a direct-minimization procedure of the total energy with respect to the density, Kohn-Sham orbitals, or density-matrix. This represents a substantial number of packages, including ONETEP,Skylaris et al. (2005) CONQUEST,Bowler et al. (2002); Gillan et al. (2007) SIESTA,Soler et al. (2002); Artacho et al. (2008) BigDFT,Genovese et al. (2011) OpenMX,Ozaki and Kino (2005) and CP2KWeber et al. (2008) (albeit that in some of these the self-consistent field technique is also available). In direct-minimization, updating the density and potential are not nested separately, so $\chi_{0}$ cannot be calculated in the manner prescribed above.

Finally, in the case of closed-shell systems, linear response approaches tend to dramatically overestimate $U$ and suffer from numerical instabilities.Hu et al. (2006); Kulik and Marzari (2010); Yu and Carter (2014); Huang (2017) We would like to overcome, or at least to better understand, these failures.

In Ref. Moynihan et al., 2017 one of us presented an alternative to the established SCF linear response approach for calculating the Hubbard parameters. This minimum-tracking linear response approach is suited for calculating $U$ in direct-minimization codes. In section II, we expand upon this formalism. To resolve the discrepancies between conventional linear response and the contemporary DFT + $U$ functional, we pay particular attention to spin and associated screening, proposing revised definitions for Hubbard and Hund’s parameters (section II.1). By comparing scalar linear response to our spin-specific theory, we demonstrate that the treatment of inter-spin screening in conventional linear response is somewhat inconsistent with the DFT + $U$ functional as it is most commonly employed (section II.2). While we do not claim here to arrive at an ultimate solution to this inconsistency, we do provide a simple technique by which inter-spin screening of the Hubbard $U$ may be suppressed. This results in spin-dependent $U$ parameters that are generally lower in value than the canonical $U$ for the partially-filled spin channel of a localized subspace (the spin channel that usually harbours the strong correlation effects) and that, in principle, could be applied to that spin channel alone. This hints at a possible solution to the widespread finding that first-principles $U$ parameters can be rather too large, in practice, leading to over-correction by DFT + $U$.

In the latter half of the paper (section III) we apply our theoretical developments to a complete set of hexahydrated transition metal complexes from Ti to Zn. We calculate Hubbard and Hund’s parameters using conventional and novel approaches (section III.2), and then perform DFT + $U$ calculations using these parameters to predict structural and spectroscopic properties (Sections III.4 and III.5). The numerical stability of the minimum-tracking formalism (in which Hubbard parameters are a strictly ground state property) allows us to investigate closed-shell cases with confidence. The Hubbard corrections to oxygen $2p$ subspaces are far from negligible, and help to obtain sensible structural predictions. Spectroscopic simulations of coordination complexes using DFT + $U$ see only mixed success, whereas our indirect band gap results for the long-standing challenge material MnO are very promising when compared against a wide range of more computationally demanding approximations.

II The minimum-tracking approach for calculating U via linear response

The minimum-tracking linear response approach is largely equivalent to SCF linear response, but its derivation centers on the ground-state density for each value of the perturbing potential.

As with the SCF approach, a perturbing potential $d\hat{v}_{\mathrm{ext}}=dv^{J}_{\mathrm{ext}}\hat{P}^{J}$ is applied to the $J$th Hubbard subspace. The response of the projected Kohn-Sham potential is given by the chain rule

[TABLE]

where the final step follows because while the external potential acting on site $J$ will change the density matrix everywhere, the $N$-site Hubbard model only sees the $N$ subspace density matrices. Screening due to the residual bath is incorporated within the total derivatives. The projections of one-body operators are given by $O^{I}=\mathrm{Tr}[\hat{P}^{I}\hat{O}]/\mathrm{Tr}[\hat{P}^{I}].$

Defining $f_{IJ}\equiv dv^{I}_{\mathrm{Hxc}}/dn^{J}$, $(\varepsilon^{-1})_{IJ}\equiv dv^{I}_{\mathrm{KS}}/dv^{J}_{\mathrm{ext}}$, and $\Omega_{IJ}\equiv dv^{I}_{\mathrm{ext}}/dv^{J}_{\mathrm{ext}}$, Eq. 8 becomes

[TABLE]

Finally, $U$ can be equated with the projected Hartree-plus-exchange-correlation kernel, with the residual bath screening in the background.Himmetoglu et al. (2014) This yields

[TABLE]

From hereon in, we will assume that $\Omega=\delta_{IJ}$. When Hubbard projectors from different atoms overlap this may become an approximation. We will also reserve $f$ for the matrix measured via linear response, and $U$ for the parameter to be subsequently used in a DFT + $U$ calculation. This distinction will become important.

Equation 10 is nothing more than a reformulation of Eq. 7. We can identify the interacting and non-interacting response matrices

[TABLE]

In this framework, we can see that the removal of the non-interacting response can be rigorously justified as a consequence of the Dyson equation, with $U$ being a measure of net interaction.

These definitions are nothing but a special case of standard linear response theory for DFT (Appendix A). It is crucial that the non-interacting response is calculated as the product of $\chi$ and $\varepsilon$, rather than $dn^{I}/dv^{J}_{\mathrm{KS}}$ directly. $dn^{I}/dv^{J}_{\mathrm{KS}}$ is both conceptually and numerically arbitrary with respect to the choice of external potential, and so its direct use must be circumvented.

Figure 1 demonstrates the calculation of elements of $\chi$ and $\varepsilon^{-1}$ from a typical set of linear response calculations.

Both minimum-tracking and conventional SCF linear response rely on the same external perturbation, and both make use of the Dyson equation. They only differ in their definition of the non-interacting response and the set of densities used in its calculation. In the minimum-tracking procedure, $\chi_{0}$ is constructed from ground-state densities of the perturbed system, and thus the resulting $U$ is strictly a ground-state property. This is obviously not the case for the SCF approach; there, $\chi_{0}$ is calculated in reference to an unconverged density and thus the resulting $U$ is not a local property of the ground-state density landscape (but still is a well-defined property of the ground-state Kohn-Sham eigen-system). This distinction is intriguing and worthy of further investigation, and possibly numerically inconsequential in practice.

Already, the minimum-tracking construction reveals an interesting property of the projected $\chi_{0}$ (and hence $U$): it is not necessarily symmetric. This is because $\chi_{0}$ as defined in Eq. 12 incorporates the total derivative of the potential, which is itself a partial derivative. While the bare $\chi_{0}$ is certainly symmetric, the response matrices that we deal with here are always screened by the background, and the screening depends on the subspace being perturbed. (In general, $\chi_{0}$ should not be symmetrized before inversion, even if the resulting $U$ matrix will be.) This observation will also hold for SCF linear response, since it also correctly goes beyond the symmetric result of first-order perturbation theory.

II.1 Accounting for spin

In the Hubbard energy functional (Eq. 1) spin and sites are treated on the same footing, with the corresponding indices being totally interchangeable. This raises the question: what happens to the response and interaction parameters if we further fine-grain linear response down to the level of spin?

In the minimum-tracking formulation it is straightforward to consider spin degrees of freedom. Response matrices become rank-four tensors

[TABLE]

and to measure these elements via linear response, we must perturb spin channels individually. (Practically, this is implemented as a combination of two potentials: a uniform shift applied to both spin-channels and a spin-splitting potential.)

This extension has several consequences. Spin-specific response functions can be visualized by flattening rank-four tensors down to rank-two ones: for example, a two-site system would have response matrices of the form

[TABLE]

This is not simply aesthetic: it means we are treating spin and atom indices on the same footing, like the DFT + $U$ functional does.

We can construct different models based on how we perform the inversion of this matrix (such as in Eq. 10): either (1) point-wise inversion, which decouples both sites and spin; (2) atom-wise inversion, with each $2\times 2$ block inverted individually, decoupling sites but not spins; or (3) invert the full matrix, leaving all sites and spins coupled. We will work through each of them in turn.

II.1.1 Point-wise inversion

The Hubbard parameters in this case are screened by the opposite spin on the same site (Fig. 2a). In this case, equation 10 separates into an independent equation for each atom:

[TABLE]

We have dropped the atomic indices for brevity. This simplification affords some numerical cancellation of errors, since inversion is no longer performed. The off-diagonal components of the matrix $f^{\sigma\sigma^{\prime}}$ are not meaningful in this case. The conventional DFT + $U$ functional requires a spin-independent $U$; for this we must average the spin-up and spin-down components:

[TABLE]

This will henceforth be referred to as “averaged 1$\times$1”. There is also the option to avoid this approximation and apply a different value of $U$ to each spin channel: $U^{\sigma}$ = $f^{\sigma\sigma}$ (“1$\times$1”).

It is interesting to note that Shishkin and SatoShishkin and Sato (2016) have previously advocated removing the off-diagonal components of site-indexed response matrices. This was motivated by the fact that these components were negligible so removing them did not alter the resulting Hubbard parameters. Here, however, the off-diagonal components components correspond to coupling between spin channels on the same atom. These components are sizeable and neglecting them appreciably alters Hubbard and Hund’s parameters, as we will see.

II.1.2 Atom-wise inversion

In atom-wise inversion, screening from both the bath and other sites is present in the response matrices, but the resulting $f=\chi_{0}^{-1}-\chi^{-1}$ is bare with respect to inter-spin interactions on the same atom as it is removed by the inversion of the spin-indexed response (Fig. 2b). Employing this approach amounts to assuming inter-spin interactions will be corrected separately i.e. with a + $J$ functional. (This is because in the absence of such a correction, a spin-screened $U$ would be necessary.)

Equation 10 reduces to

[TABLE]

where each term is a two-by-two matrix indexed by spin channel, and if there are $N$ atoms there are $N$ such equations. For practical use in DFT + $U$ + $J$, $f$ can be related to the scalar Hubbard parameter $U$ that, in the minimum-tracking linear-response formalism, is defined by

[TABLE]

There are two alternative approximations we can make here. The first, more naïve approach, is to further approximate this as

[TABLE]

which we will refer to as “simple 2$\times$2”. A more sophisticated approach (scaled “$2\times 2$”) is

[TABLE]

The derivations of Eqs. 19 and 20 involve varying levels of approximation, which are discussed in detail in Appendix B.

With atom-wise inversion, Hund’s parameters $J$ can be directly calculated in an analogous manner to $U$: in place of Eq. 18 we instead define, within the spin-polarized minimum-tracking linear response formalism,

[TABLE]

For $J$, simple $2\times 2$ yields

[TABLE]

while scaled $2\times 2$ gives

[TABLE]

II.1.3 Full inversion

Finally, in the case of full matrix inversion, the result is bare with respect to both inter-spin and inter-site interactions by the same logic. This implies that inter-atom interactions require, and are subject to, correction via a + $V$ term. This $V$ term would be doubly spin-dependent, and it may need to be symmetrized with respect to the site indices to retain a Hermitian Kohn-Sham Hamiltonian for each spin. We will not explore this approach further in this work.

We emphasize that including each of these successive terms ($J$ and $V$) should not be viewed as systematic improvements. In the limit that corrective parameters are introduced within and between every single subspace (such that the corresponding screened interactions are removed) the entire system becomes effectively non-interacting. Corrective terms are only appropriate where the corresponding interactions dwarf all others.

II.2 Comparisons with the conventional scalar approach

Conventional linear response calculations do not treat spin channels separately (Fig. 2c); for a single-site system $\chi$, $\varepsilon^{-1}$ and $f$ would all be scalars. It is straightforward to relate the spin-indexed response matrices of the previous section to these scalars:

[TABLE]

Likewise

[TABLE]

These two relations allow us to examine the role of spin-screening in scalar linear response. The Hubbard parameter obtained via spin-indexed, atom-wise inversion (scaled $2\times 2$; Eq. 20) can be rewritten as

[TABLE]

This is nothing less than the scalar expression $U=\chi_{0}^{-1}-\chi^{-1}$, which is used in scalar linear response. We may conclude that the conventional scalar approach and scaled $2\times 2$ are entirely equivalent.

Therefore, Hubbard parameters obtained by spin-aggregated approaches are not screened by the opposite spin channel on the same site. Since they combine both like and unlike spin interactions (c.f. Eq. 20), they do not correspond to the like-spin-only interaction $U_{\text{eff}}=U-J$ (as implied elsewhere).Kulik and Marzari (2008) We could have anticipated this result: during a scalar linear response calculation there is no shift in the external potential difference between the two spin channels, so (to first order) there is no external driver for changes in subspace spin polarization.

We noted earlier that atom-wise inversion formally necessitates a Hund’s correction, but such a correction is not usually included when the conventional linear response approach is employed. Given that these methods are equivalent, we argue that it is more consistent to include a Hund’s exchange correction term (e.g. calculated using Eq. 22) if using a Hubbard correction calculated in the conventional manner.

The precise functional form of the + $J$ correction needed is, however, the subject of ongoing research. Recently, for example, Millis and co-workers demonstrated that spin-polarized DFT already possesses some degree of intrinsic exchange splitting, and they have argued convincingly that the contemporary form of the + $J$ correction can overestimate exchange splitting.Chen and Millis (2016) This finding is corroborated by our own results discussed later in this paper (e.g. Table 7).

III Application to a complete series of hexahydrated transition metals and manganese oxide

In the second half of this work, we explore the ramifications of our theoretical developments on two test systems: hexahydrated transition metals, and manganese oxide.

In these systems, all of the metal atoms have partially filled $3d$ sub-shells. Electrons within these sub-shells are in such close proximity to one another that the interplay of their spin, charge, and orbital moment are too pronounced to be well described by local or semi-local xc-functionals.Harvey (2006); Cramer and Truhlar (2009); Yu et al. (2016); Johnson and Becke (2017) DFT + $U$ may provide a more accurate description of these systems.Kulik and Marzari (2010, 2011a); Himmetoglu et al. (2014); Zhao et al. (2016)

Manganese oxide (MnO) has a rock salt structure. At low temperatures it is antiferromagnetic,Shull et al. (1951) and has a band gap of approximately 4 eV that is substantially underestimated by semi-local functionals.Anisimov et al. (1997) Conventional linear-response calculations on MnO yield an excessively large Hubbard parameter ($U>7$ eV).Lim et al. (2016)

Meanwhile, hexahydrated transition metals comprise of a central first-row transition metal ion surrounded by six water ligands in a tetragonal arrangement (Fig. 3c). Such systems bear some resemblance to a fundamental unit of transition metal oxides, as well as organometallic systems such as the oxygen evolving complex of photosystem II.Umena et al. (2011); Suga et al. (2015)

Depending on the electronic structure of the metal, these systems may exhibit Jahn-Teller distortion, resulting in an elongated tetragonal structure with two axial waters being slightly more distant than their four equatorial counterparts (Fig. 3).

III.1 Computational details

All calculations were performed using ONETEPSkylaris et al. (2005); Hine et al. (2011a); O’Regan et al. (2012, 2010); Ruiz-Serrano et al. (2012); Fox et al. (2011); O’Regan et al. (2011) (Order-$N$ Electronic Total Energy Package, version 4.3) using the Perdew-Burke-Ernzerhof (PBE) xc-functional.Perdew et al. (1996)

For MnO, a square super-cell containing 512 atoms was simulated under periodic boundary conditions without explicit $k$-point sampling. This is a non-diagonal super-cellLloyd-Williams and Monserrat (2015) of the four-atom primitive cell, and gives an equivalent $k$-point sampling scheme that includes both $Z$ and $\Gamma$. (This is crucial because the band gap of MnO is known to be $Z$ to $\Gamma$.) The lattice parameter was set to the experimental value of 4.445 Å.Wyckoff (1963) The calculations were spin-polarized, with an energy cut-off of 1030 eV. ONETEP uses a basis of non-orthogonal generalized Wannier functions (NGWFs) that are variationally optimized in situ. Each Mn atom had ten NGWFs; O atoms, four. All NGWFs had a cutoff radius of 11.0 $a_{0}$.

For the hexahydrated metals, all calculations were spin-polarized, with an energy cut-off of 897 eV. Depending on the species, there were 9, 10, or 13 NGWFs on the transition metal atom, four on each oxygen, and one on each hydrogen. All NGWFs had 14 $a_{0}$ cutoff radii. An Elstner dispersion correctionElstner et al. (2001); Hill and Skylaris (2009) was applied, and electrostatics were treated using a padded cell and a Coulomb cut-off.Hine et al. (2011b)

For all the calculations, the Hubbard projectors were constructed from solving the neutral atomic problem subject to the pseudopotential of the species in question.Ruiz-Serrano et al. (2012) Most pseudopotentials were taken from the Rappe group pseudopotential libraryRap although those for Co and Fe were generated in-house using OPIUM.opi ; Kerker (1980); Kleinman and Bylander (1982); Hamann (1989); Rappe et al. (1990); Gonze et al. (1991); Ramer and Rappe (1999) These were scalar relativistic pseudopotentialsGrinberg et al. (2000) with non-linear core corrections.Louie et al. (1982). All DFT + $U$ + $J$ calculations used a $+J$ correction to the energy, potential, and ionic forces. We used the energetic correction shown in Eq. 3 (following the example of Ref. Himmetoglu et al., 2011 we have omitted the “$n_{\mathrm{min}}$” term that appears in that paper).

Example input and output files can be found at www.repository.cam.ac.uk/.

III.2 Calculating Hubbard parameters