# Orbital-resolved DFT+U for molecules and solids

Eric Macke  Faculty of Production Engineering, Bremen Center for Computational
Materials Science and MAPEX Center for Materials and Processes, Hybrid
Materials Interfaces Group, University of Bremen, Am Fallturm 1, 28359 Bremen,
Germany<EMAIL_ADDRESS>Iurii Timrov  Theory and Simulation of Materials
(THEOS) and National Centre for Computational Design and Discovery of Novel
Materials (MARVEL), École Polytechnique Fédérale de Lausanne, CH-1015
Lausanne, Switzerland  Nicola Marzari  Theory and Simulation of Materials
(THEOS) and National Centre for Computational Design and Discovery of Novel
Materials (MARVEL), École Polytechnique Fédérale de Lausanne, CH-1015
Lausanne, Switzerland  Lucio Colombi Ciacchi  Faculty of Production
Engineering, Bremen Center for Computational Materials Science and MAPEX
Center for Materials and Processes, Hybrid Materials Interfaces Group,
University of Bremen, Am Fallturm 1, 28359 Bremen, Germany

###### Abstract

We present an orbital-resolved extension of the Hubbard $U$ correction to
density-functional theory (DFT). Compared to the conventional shell-averaged
approach, the prediction of energetic, electronic and structural properties is
strongly improved, particularly for compounds characterized by both localized
and hybridized states in the Hubbard manifold. The numerical values of all
Hubbard parameters are readily obtained from linear-response calculations. The
relevance of this more refined approach is showcased by its application to
bulk solids pyrite (FeS2) and pyrolusite (β-MnO2), as well as to six Fe(II)
molecular complexes. Our findings indicate that a careful definition of
Hubbard manifolds is indispensable for extending the applicability of DFT+$U$
beyond its current boundaries. The present orbital-resolved scheme aims to
provide a computationally undemanding yet accurate tool for electronic
structure calculations of charge-transfer insulators, transition-metal (TM)
complexes and other compounds displaying significant orbital hybridization.

###### keywords:

American Chemical Society, Hubbard, DFT+U, DFT+U+V, LDA+U, Quantum Espresso,
linear response, constrained DFT, manifold, low-spin, orbital,

Present address: Laboratory for Materials Simulations (LMS), Paul Scherrer
Institut (PSI), CH-5232 Villigen PSI, Switzerland University of Bremen
Excellence Chair, Bremen Center for Computational Materials Science
DFT,DFT+U,DFT+U+V,DFT+U+J,cDFT,LR-
cDFT,LS,HS,FCE,SIE,PWL,KS,MO,TM,NAO,OAO,PAW,MLWF,DFPT,VBM,CBM,CASPT2/CC,AFM,KI,LDA,GGA,DMFT

##  1 Introduction

Hubbard corrections are among the most widely used improvements to approximate
Kohn-Sham (KS) density-functional theory (DFT) 1, 2. The combination of DFT
with the Hubbard model 3, referred to as DFT+$U$, was inspired by this widely
studied model of electron correlations and introduced to improve the
description provided by local or semi-local exchange-correlation (xc)
functionals, such as the local-density approximation (LDA) and the
generalized-gradient approximation (GGA), for the case of Mott-Hubbard
insulators 4, 5, 6 or more broadly strongly correlated electrons.

It was recognized early on 7, 8 that a simplified, rotationally invariant
formulation of DFT$+U$ 9 provided a natural connection to the requirement of
piecewise linearity of the exact energy functional10, 11: the Hubbard
correction of Dudarev effectively removes the non-linear (almost quadratic)
behavior of the total energy with respect to the occupation of the Hubbard
manifold, and replaces it with a linear term. In this light, the strength of
the $U$ parameter can be determined fully from first principles, and obtained
from the second derivative of the energy (once the non-interacting terms are
removed) with respect to the occupations 7, 8, so that the quadratic curvature
is removed by DFT$+U$. This connection is heuristic, and relies on the (very
reasonable) assumption that the localized electrons in the $d$ or $f$ manifold
are only weakly interacting with the rest of the electron bath, so that this
manifold also follows the condition of piecewise linearity (PWL) that is in
principle valid only for the total energy functional10. The overstabilization
of fractional occupations in standard (semi)-local functionals is driven by
the incomplete cancellation from the xc functional to the Hartree term
(cancellation that is exact instead in Hartree-Fock), leading to a strong
(one-electron) self-interaction component. This is particularly severe for
localized electrons, that over-delocalize and over-hybridize with their
ligands. So, the DFT$+U$ correction should be seen as a self-interaction
correction, providing an approximate screened Fock-like contribution; this was
stated early on, highlighting how even for molecular systems containing only
one site (i.e., one TM atom) the electronic-structure description provided by
DFT$+U$ is greatly improved, both qualitatively and quantitatively 12. In
other words, while the original Hubbard model was concerned with the
correlation effects that appear on a lattice, where each site can only contain
0, 1 or 2 electrons, the functional form that the model inspires, once applied
to the continuum of the electron gas of approximate DFT, does not any more
account for correlations between sites, but counteracts the tendency of
strongly localized electrons to hybridize with their ligands, driven by
incomplete cancellation in the Hartree term of the one-electron self-
interaction. This can be easily argued not only by applying DFT$+U$ to
molecules 12, but also by the simple observation that changing the value of
$U$ changes the charge transfer to/from the ligands13, but does not transfer
charge to the Hubbard manifolds on the other sites, highlighting the
completely different physics of the Hubbard Hamiltonian on a lattice, and the
Hubbard $U$ correction in the continuum of Kohn-Sham electrons. Incidentally,
strongly correlated materials have almost invariably very localized $d$ or $f$
electrons, so do benefit from DFT$+U$ – but not because the latter improves
the treatment of correlations, but just because it decreases the self-
interaction errors (SIE). It is important to iterate that PWL is a property
that a system must obey on the global scale, i.e., with respect to the changes
in the total number of electrons14. DFT$+U$, on the other hand, seeks to
eliminate the curvature locally by restoring the piece-wise linear behavior in
a subspace known as the Hubbard manifold 8, 15. The relation between self-
interaction and piecewise linearity in many-electron systems is somewhat
confusing; in fact, Koopmans functionals 16, 17 impose (very accurately)
piecewise linearity on each orbital in a given system, but do not correct
self-interaction; the integer Koopmans functional (KI) provides the same exact
total energy of the base functional (LDA or GGA), but a completely different
spectrum.

A key challenge in DFT$+U$ arises from the (a priori) choice of the value of
the on-site Hubbard $U$ parameter. Often, $U$ is considered a tunable
quantity, adjusted to achieve agreement with experimental results for specific
properties of interest such as band gaps, lattice parameters, magnetic
moments, or formation enthalpies 18, 19, 20, 21. However, this empirical
procedure possesses limited predictive capabilities due to its reliance on
experimental data; not to mention that certain quantities, like band gaps, are
not even meant to be predicted by DFT. Alternatively, $U$ values can be
computed from first principles. The different methods that have been proposed
for this purpose can be classified into three groups: the constrained DFT
approach (cDFT) 22, 23, 24, 25, 26, 27, 28, 29, 30, the constrained random
phase approximation (cRPA) 31, 32, 33, 34, 35, and Hartree-Fock based
approaches 36, 37, 38, 39, 40, 41. The linear-response formulation of
constrained DFT (LR-cDFT) 8 has become a method of choice for many DFT+$U$
studies 42, 43, 44, 45, 46, 47, 48, 49, 50. A recent reformulation of this
method in terms of density-functional perturbation theory (DFPT) 51, 52 has
significantly enhanced its success. This reformulation enables the replacement
of computationally expensive supercells with primitive unit cells, utilizing
monochromatic perturbations. As a result, the computational burden of
determining Hubbard parameters is substantially reduced.

As evident from the large variety of available methods, considerable attention
has been devoted to the numerical evaluation of $U$ for a given manifold.
However, the question how this manifold must be defined, that is, which
electronic states require on-site Hubbard corrections, is a critical yet often
overlooked aspect in DFT$+U$. Recalling that the main motivation of Hubbard
$U$ corrections lies in the mitigation of local SIE through recovery of PWL of
the total energy, the Hubbard manifold should contain those and only those
states that substantially contribute to the former. Oftentimes, self-
interaction occurs in partially occupied $d$ and $f$ shells due to their high
electron count and localization; hence, these are the traditional targets of
Hubbard $U$ corrections. Nevertheless, self-interaction can also manifest
itself in $s$ and $p$ shells, and even in localized molecular orbitals (MO)53.
Moreover, the definition of the Hubbard manifold also involves choosing so-
called Hubbard-projector functions that suitably represent the states to be
corrected (details follow in Section 2.2). In fact, it has been shown that the
choice of the projector function can have a stronger impact on the calculation
than the numerical value of $U$, especially when orbital overlap is
appreciable54, 45, 55. Therefore, it is essential to choose a $U$ value that
is consistent with the projectors used (see appendix of Ref. 56).

The significance of judiciously selecting the Hubbard manifold becomes
particularly evident in the context of optical properties. Although DFT is a
total energy theory and KS orbitals have no direct physical meaning except for
the highest occupied molecular orbital (HOMO), KS eigenvalues are commonly
used as estimates for experimental quasiparticle energies. Unfortunately,
uncorrected (semi-)local DFT functionals tend to severely underestimate
properties such as band gaps and ionization energies; a fact that has been
linked to deviations from PWL 16, 57, 17. Thus, Koopmans-compliant
functionals16, 17, 58 systematically improve spectral properties by imposing
PWL on all orbitals in the system. For example, a recent study by Nguyen et
al. found an average band gap error of only 0.22 eV for a test set of 30
semiconductors59. Conversely, DFT$+U$ can only achieve such improvements if
the Hubbard manifold sufficiently overlaps with the frontier orbitals 60, 15,
45. This might explain why DFT$+U$ has been successfully applied to many Mott-
Hubbard insulators, while its predictive capabilities are less reliable when
applied to charge-transfer (CT) insulators. While in the former both the
valence band maximum (VBM) and the conduction band minimum (CBM) primarily
consist of localized TM $d$ states, the frontier states of the latter possess
significant ligand orbital character. In this case, the sole correction of TM
$d$ shells does not necessarily target all major sources of SIE. Additionally,
the representation of the hybridized orbitals through atomic-like Hubbard
projectors can be deceptive.

Some authors have explored including ligand orbitals in the on-site Hubbard
manifold, for example applying $U$ corrections to oxygen and S-$p$ shells of
TM oxides, perovskites and sulfides 48, 61, 62, 45, 63. The extended
DFT+$U$+$V$ framework53 was adopted for scenarios where hybridization plays an
eminent role. This method augments DFT$+U$ by an inter-site Hubbard $V$ term,
and so acts on combinations of projectors located on different sites,
enhancing accuracy and transferability 64, 65, 66, 67, 68, 54, 69, 70, 71, 72,
41, 73. Finally, many open-shell systems require unlike-spin terms such as
Hund’s $J$. In DFT+$U$+$J$, these terms introduce anisotropy in the Hubbard
correction with respect to the spin channels, while also reducing the
effective value of the $U$ parameter 60, 74, 75, 48, 62. On the other hand, it
is unlikely that multi-reference configurations, as common in open-shell
systems, might be described well with any $+U$ or $U$+$J$ correction.

What is striking about the aforementioned approaches is that the $d$ shell
Hubbard manifold is either extended by additional states ($U$ on ligands or
$V$ on molecular orbitals) or internally rebalanced (DFT+$U$+$J$). But what if
the practice of including the entire $d$ shell is, in itself, problematic?
Recent work by Mariano et al. suggests that this might indeed occur in some
cases. Their study revealed that applying $U$ corrections to the $d$ shells of
Fe in strong-field Fe(II) hexacomplexes leads to a spurious suppression of LS
states due to what appears to be an over-correction of hybridized $e_{g}$
orbitals 76, 77. Furthermore, the use of ab initio Hubbard $U$ parameters
derived from LR-cDFT led to a significant decrease in the overall accuracy of
the spin-state energies when compared to empirical $U$ values. These and other
observations suggest a critical analysis of the current practice of how
(shell-averaged) Hubbard $U$ corrections are applied to DFT. After all, the
correction of all magnetic quantum orbitals within a given shell using the
same scalar $U$ parameter is inherently a simplistic approximation 78.

Therefore, it is worthwhile to investigate whether departing from the shell-
averaged approximation can improve DFT$+U$ calculations concerning energetic,
structural, and magnetic properties. This possibility has been explored in a
few works. For example, Solovyev et al. showed that a selective application of
$U$ corrections to the $t_{2g}$ manifolds of LaMO3 perovskites (M=Ti-Cu)
significantly improves the agreement of band gaps and magnetic orderings with
experiments compared to both uncorrected LDA and shell-averaged LDA$+U$ 79.
The authors argue that $e_{g}$ electrons are reasonably described by the
uncorrected LDA functional because of their “itinerant“ behavior that arises
due to the strong $\sigma$ overlap between O-$2p$ and M-$e_{g}$ orbitals.
Pickett et al. adopted the cDFT approach to compute the Hubbard parameters of
Fe in bulk FeO as a matrix ($\mathbf{U}$) with $U$ values specific to the
$t_{2g}$ and $e_{g}$ manifolds, respectively 27. Using a different scheme,
orbital-resolved Hubbard parameters were obtained by mapping shell-averaged
$U$ and $J$ parameters onto orbital-dependent interaction and exchange
matrices $U_{mm^{\prime}}$ and $J_{mm^{\prime}}$ using atomic Slater integrals
and Gaunt’s numbers 80, 81, 82, 83. Other works focused on the spin-dependence
of $U$, finding that the spin-resolved on-site parameters can be pivotal for a
physical description of magnetic systems 50, 84, 49.

Despite these early efforts, orbital-resolved Hubbard $U$ parameters have not
gained widespread use in the DFT community. This is surprising in light of the
fact that orbital-specific Hubbard manifolds are quite common in the dynamical
mean-field theory (DMFT) community 85, 86, 87, 88. In the latter, Hubbard $U$
parameters are routinely computed for manifolds that range from groups of
orbitals (e.g., $t_{2g}$ models) to combinations of multiple shells localized
on different atoms (e.g., $d-p$ models) 86 using cRPA in conjunction with
various sets of Wannier projector functions that encompass the specific
manifold of interest 31, 32. To date, only a few publicly available DFT codes
incorporate this capability. For example, DFT codes that use Wannier function-
based Hubbard projectors (e.g., ONETEP89) indirectly facilitate orbital-
resolved DFT$+U$ calculations as the Wannier functions can be chosen to only
represent the desired subset of orbitals, e.g., $t_{2g}$55.

In this paper we present a user-friendly yet general implementation of
orbital-resolved DFT$+U$ that works with any kind of Hubbard projector. The
numerical values of the necessary parameters are extracted from first
principles using an orbital-resolved LR-cDFT approach 8. We benchmark the
orbital-resolved scheme by carrying out calculations of bulk pyrite (FeS2),
bulk pyrolusite (β-MnO2) and six Fe(II) molecular hexacomplexes of varying
ligand strength. For all of these strongly covalent compounds, the refined
approach leads to a substantial improvement in the prediction of structural
and energetic properties, aligning more closely with experimental observations
than conventional DFT$+U$ or even DFT${+U\\!+\\!V}$. Moreover, we observe that
the orbital-resolved $U$ parameters are considerably smaller, by up to 80%,
than their corresponding shell-averaged counterparts.

The remainder of the paper is organized as follows. In Section 2, we describe
the transition from the customary, shell-averaged implementations of the
DFT$+U$ energy functional and the LR-cDFT approach to their generalized,
orbital-resolved forms. Subsequently, Section 3 briefly lists the relevant
technical details of the calculations, whose results are presented and
discussed in Section 4. Section 5 is dedicated to the relationship between the
orbital-resolved (on-site) Hubbard $U$ presented in this work and the (inter-
site) $V$ terms of Ref. 53. Finally, we summarize the main conclusions in
Section 6.

##  2 DFT$+U$: orbital-resolved Hubbard parameters

###  2.1 The energy functional

As this work focuses on the on-site $U$ term of Hubbard-corrected DFT, we
start from the widespread shell-averaged and rotationally-invariant
formulation of DFT+$U$ by Dudarev et al.9. Written in a way that allows for
the simultaneous correction of multiple subshells on the same atom, the energy
functional reads90:

$E_{\mathrm{DFT}+U}=E_{\mathrm{DFT}}+E_{U}=E_{\mathrm{DFT}}+\sum_{I,\sigma}\sum_{nl}\frac{U^{I}_{nl}}{2}\,\mathrm{Tr}\left[\mathbf{n}^{I\sigma}_{nl}(1-\mathbf{n}^{I\sigma}_{nl})\right]\,,$
(1)

where $E_{\mathrm{DFT}}$ is the DFT total energy computed with standard
(semi-)local xc functionals, $U^{I}_{nl}$ is an effective on-site Hubbard
parameter, and $\mathbf{n}^{I\sigma}_{nl}$ is the orbital occupation matrix.
The summation over the principal ($n$) and orbital ($l$) quantum numbers
implies that Hubbard $U$ corrections can be applied simultaneously to multiple
subshells of the same atom53, 45, although it usually suffices to treat the
valence shell alone. For the sake of clarity, we henceforth omit the index
$nl$, assuming that the Hubbard manifolds consist of, at most, one subshell
per atom. The orbital occupation matrix defines how the Hubbard correction is
applied. For a given spin $\sigma$, its elements are computed by projecting
the valence KS wave functions onto the Hubbard manifold of an atom $I$:

${n}^{I\sigma}_{mm^{\prime}}=\sum_{\mathbf{k},\nu}f_{\mathbf{k},\nu}^{\sigma}\langle\psi_{\mathbf{k},\nu}^{\sigma}|\hat{P}^{I}_{m^{\prime}m}|\psi_{\mathbf{k},\nu}^{\sigma}\rangle\,,$
(2)

where $m$ and $m^{\prime}$ are the magnetic quantum numbers of the Hubbard
manifold, $\nu$ represents the band labels of the KS wave functions,
$\mathbf{k}$ indicates points in the first Brillouin zone, and
$f_{\mathbf{k},\nu}^{\sigma}$ are occupations of the KS wave functions
$\psi^{\sigma}_{\mathbf{k},\nu}$. If atomic-like orbitals ($\varphi^{I}_{m}$)
are used as projector functions (vide infra), one may define
$\hat{P}^{I}_{m^{\prime}m}=|\varphi^{I}_{m^{\prime}}\rangle\langle\varphi^{I}_{m}|$.
Note, however, that if plane-wave basis sets are used, the expression of
$\hat{P}$ also depends on the type of the pseudopotential 52. Applying Hubbard
$U$ corrections also modifies the KS potential of the target orbitals
according to

$\hat{V}_{U}^{\sigma}=\sum_{I}U^{I}\sum_{mm^{\prime}}\left(\frac{\delta_{mm^{\prime}}}{2}-n^{I\sigma}_{mm^{\prime}}\right)|\varphi^{I}_{m}\rangle\langle\varphi^{I}_{m^{\prime}}|,$
(3)

where $\delta_{mm^{\prime}}$ is the Kronecker delta. It is evident from Eq.
(3) that Hubbard $U$ corrections exert a stabilizing influence on fully
occupied orbitals, reducing their KS potential by up to $U/2$, while producing
the opposite effect on empty orbitals.

Provided a correct normalization of the Hubbard projector functions, the
eigenvalues of $n^{I\sigma}_{mm^{\prime}}$ express the occupation of the
$2l+1$ orbitals with numbers between 0 (fully empty) and 1 (fully occupied).
They are obtained by solving the eigenvalue problem 8

$\mathbf{n}^{I\sigma}\mathbf{v}_{i}^{I\sigma}=\lambda_{i}^{I\sigma}{\mathbf{v}}_{i}^{I\sigma},$
(4)

where $\lambda_{i}^{I\sigma}$ and $\mathbf{v}_{i}^{I\sigma}$ are the
eigenvalues and eigenvectors, respectively, and $i$ is a dummy index running
from 1 to $2l+1$. This allows us to re-write the $E_{U}$ contribution of Eq.
(1) more concisely in terms of these eigenvalues:

$E_{U}=\sum_{I,\sigma}\frac{U^{I}}{2}\sum_{i}^{2l+1}\lambda_{i}^{I\sigma}(1-\lambda_{i}^{I\sigma})\,.$
(5)

Within this diagonal representation, Eq. (5) can be readily generalized to
become an orbital-resolved correction:

$E_{U}=\sum_{I,\sigma}\sum_{i}^{2l+1}\frac{U_{i}^{I}}{2}\lambda_{i}^{I\sigma}(1-\lambda_{i}^{I\sigma})\,,$
(6)

Here, $U_{i}^{I}$ is now an effective on-site Hubbard parameter specific to
the $i$th orbital of the Hubbard manifold localized on atom $I$.

Similar generalizations of the shell-averaged DFT+$U$ functional were already
postulated by Pickett et al. and used by Solovyev et al., although these early
works adopted DFT$+U$ energy functionals different from Dudarev et al.’s
formulation.

In practical calculations, not every orbital requires a distinct $U_{i}^{I}$,
and it might suffice to distinguish between the irreducible representations
that follow from local point group symmetry. For instance, considering the $d$
shell of an octahedrally coordinated atom with local $O_{h}$ point group
symmetry, Eq. (6) can be written as

$E_{U}=\sum_{I,\sigma}\frac{U_{t_{2g}}^{I}}{2}\sum_{i\in\\{t_{2g}\\}}\lambda_{i}^{I\sigma}(1-\lambda_{i}^{I\sigma})+\sum_{I,\sigma}\frac{U_{e_{g}}^{I}}{2}\sum_{i\in\\{e_{g}\\}}\lambda_{i}^{I\sigma}(1-\lambda_{i}^{I\sigma})\,,$
(7)

where $U_{t_{2g}}^{I}$ and $U_{e_{g}}^{I}$ are the Hubbard parameters for the
$t_{2g}$ and $e_{g}$ orbitals, respectively, and $i\in\\{t_{2g}\\}$ means that
$i$ runs over the orbital indices of the $t_{2g}$ subshell, while
$i\in\\{e_{g}\\}$ means that $i$ runs over the orbital indices of the $e_{g}$
subshell. It is also possible to selectively exclude specific orbitals from
receiving Hubbard corrections by setting $U^{I}_{i}=0$.

The calculation of the Hubbard potential (Eq. (3)) is also performed in the
diagonal representation by rotating the atomic-like orbitals $\varphi^{I}_{m}$
from the global to the local coordinate system using the eigenvectors
$\mathbf{v}_{i}^{I\sigma}$. However, after computing the orbital-resolved
contributions to the Hubbard potential $\hat{V_{U}}$, we perform a
backrotation into the non-diagonal representation as this allows to use
existing implementations of density mixing as well as the calculation of
forces and stresses with no further adaption.

A crucial aspect for a successful application of Hubbard corrections lies in
finding a suitable projector function for the target manifold, as this
controls the occupation eigenvalues $\lambda_{i}^{I\sigma}$ (Eq. (2)) which
govern the corrective Hubbard energy. The next subsection offers a concise
introduction to the prevailing and commonly utilized approaches.

###  2.2 Hubbard projector functions

There are many ways to define Hubbard projector functions within the DFT+$U$
approach, and the reader is referred to Ref. 91 and references therein for a
more comprehensive overview. For electronic-structure codes employing a
localized basis set, natural choices for Hubbard occupations are either
Mulliken or Löwdin population matrices 92. In contrast, codes based on plane-
wave basis sets often use atomic orbitals as projectors 8, 93, whose localized
functions are parameterized with free-atom calculations and then stored in the
pseudopotentials. During the generation of pseudopotentials, the atomic
orbitals are chosen to be orthonormal to all other orbitals centered on the
same atom. However, this cannot guarantee orthogonality to orbitals localized
on other sites during practical calculations. Hence, these projector functions
are referred to as nonorthogonalized atomic orbitals (NAO) 54. While
straightforward to implement and use, NAO may display spatially extended
‘tails’, potentially resulting in the same domain being tackled twice by
Hubbard corrections if the atomic orbitals of two neighboring atoms overlap.

Truncation spheres 93 provide a means of ruling out such double counting by
cutting off possible tails. 111Note that this “double counting” must not be
confused with the double counting term in DFT+$U$ energy functionals derived
from the Hubbard model. However, the radius of these spheres represents an
additional parameter that affects the obtained occupation numbers 94. The
value of this parameter varies depending on the code utilized and may require
manual adjustment (e.g., muffin-tin radius in the linearized augmented plane-
wave approach) or can be embedded within the pseudopotential (typical of
projector-augmented-wave (PAW) pseudopotentials).

An alternative, parameter-free approach to circumvent truncation involves
orthogonalizing NAO across all atomic sites, e.g., using Löwdin’s scheme 95,
thus transforming them into orthogonalized atomic orbitals (OAO) 54, 91. The
resulting inter-site orthogonality clears the overlap of projectors on
different sites 91 and even accounts for possible hybridization between them,
albeit to a limited extent. Several benchmark studies have shown that OAO
projectors consistently outperform NAO projectors with respect to structural,
electronic and spectral properties54, 69, 45.

Finally, Wannier functions are also viable Hubbard projector functions 29, 96,
97. Specifically, maximally localized Wannier functions (MLWF) 98, 99 can
separate manifolds in a system-specific fashion 100, and can serve as
effective Hubbard projector functions within the generalized DFT+$U$ framework
presented in Eq. (6). Despite the extensive use of Wannier functions, their
current adoption in DFT+$U$ calculations remains limited 101, 102, 55,
possibly due to the non-trivial steps of finding an appropriate starting guess
and disentangling overlapping bands during wannierization. Moreover, features
relevant for practical studies such as the calculation of forces and stresses
are cumbersome to implement.

###  2.3 LR-cDFT to compute orbital-resolved $U$ parameters

LR-cDFT is based on using the DFT$+U$ energy functional to (heuristically)
restore piecewise-linearity for the Hubbard manifold in (semi-)local DFT
functionals suffering from electron self-interaction8, 12. An important
manifestation of the latter are so-called fractional charge errors (FCE),
which are spurious (usually convex) deviations from linearity of
$E_{\mathrm{DFT}}$ with respect to fractional addition or removal of charge
15, i.e.,

$\mathrm{FCE}=\frac{\partial^{2}E_{\mathrm{DFT}}}{\partial q^{2}},$ (8)

where $q$ represents the charge of the system under consideration. Inspection
of Eq. (1) shows that the Hubbard $U$ correction amounts to removing a
quadratic term and adding a linear one, scaled by the numerical value of $U$.
Note, however, that the curvature removed by Eq. (1) is not with respect to
the total charge of the system $q$ but instead with respect to the projected
local occupation of the Hubbard manifold $\mathbf{n}^{I\sigma}$. Thus, a
fundamental assumption is that the electrons in this localized Hubbard
manifold are the most affected by self-interaction and can be dealt with
separately. This interpretation makes it possible to define the value of
Hubbard $U$ as the one for which the second derivative of the total energy
functional becomes zero with respect to changes in the occupation of the shell
8,

$U^{I}=\frac{\partial^{2}E}{\partial(\Lambda^{I})^{2}}\rvert_{q}\,,$ (9)

where we now use
$\Lambda^{I}=\sum_{\sigma}\sum_{i}^{2l+1}\lambda^{I\sigma}_{i}$ to define the
local occupation of the Hubbard manifold. Because a direct control of orbital
occupations is not tractable in codes that obtain them as output quantities,
Lagrange multipliers $\alpha$ are introduced to linearly shift the potential
of the Hubbard manifold and thus indirectly control $\Lambda^{I}$ (see Refs.
27, 7, 8 for the derivation). Then, two response matrices are defined using
finite differences:

$(\chi)^{IJ}=\frac{\Delta\Lambda^{I}}{\Delta\alpha^{J}}\,,\quad(\chi_{0})^{IJ}=\frac{\Delta\Lambda^{I}_{0}}{\Delta\alpha^{J}}\,,$
(10)

where $\chi$ represents the self-consistent, screened response of the
manifold, whereas $\chi_{0}$ is due to the non-interacting, unscreened
response coming from the re-hybridization of the atomic orbital projectors
that results from the perturbation. In many plane-wave codes, these latter
non-interacting orbital occupations $\Lambda^{I}_{0}$ can be obtained from the
first iteration in the self-consistent cycle of the perturbative calculation.
Since the response $\chi_{0}$ is unrelated to electron self-interaction, it
must be subtracted from the interacting (i.e., screened) response when
computing $U$:

$U^{I}=\left(\chi_{0}^{-1}-\chi^{-1}\right)^{II}\,.$ (11)

In practice, the response functions are obtained by either applying multiple
small (positive and negative) perturbations to the shells of interest of a
converged ground state 8 (with periodic systems requiring a supercell approach
to avoid interactions between a perturbed Hubbard manifold and its periodic
images) or from DFPT 51, 103 using the response to monochromatic perturbations
in a primitive cell.

Orbital-resolved Hubbard $U$ parameters can be evaluated using the formalism
of Eqs. (10) and (11) by adapting the definition of the total occupation of
the Hubbard manifold. In the most general case, every magnetic quantum orbital
$i$ of the $nl$ subshell can acquire an individual Hubbard parameter. The
occupation of such a manifold is given by

$\Lambda^{I}_{i}=\sum_{\sigma}\lambda^{I\sigma}_{i}\,.$ (12)

With this, the elements of the response matrices (Eq. (10)) can be redefined
as

$(\chi)^{IJ}_{ij}=\frac{\Delta\Lambda^{I}_{i}}{\Delta\alpha^{J}_{j}}\,,\quad(\chi_{0})^{IJ}_{ij}=\frac{\Delta\Lambda^{I}_{i,0}}{\Delta\alpha^{J}_{j}}\,,$
(13)

while the expression for the orbital-resolved on-site Hubbard $U$ parameters
becomes

$U^{I}_{ij}=\left(\chi_{0}^{-1}-\chi^{-1}\right)^{II}_{ij}\,.$ (14)

The generalized formalism of Eq. (14) not only accommodates orbital-resolved
on-site $U$ parameters, but allows for the determination of inter-site
parameters ($U^{IJ}_{ij}$ with $I\neq J$, or, following the nomenclature of
DFT${+U\\!+\\!V}$, $V^{IJ}_{ij}$), as well as on-site inter-orbital parameters
($U^{II}_{ij}$ for $i\neq j$). Nevertheless, the primary emphasis of this
study lies in investigating the on-site intra-manifold parameters
$U^{II}_{ii}\equiv U^{I}_{i}$.

For the sake of illustration, let us consider a d shell of an atom exhibiting
local $O_{h}$ symmetry. Also, we assume that this atom be the only Hubbard
atom in the system, thereby enabling us to neglect inter-site responses and to
drop the superscript $I$. Then, one can define an orbital-resolved Hubbard
matrix of size $2\times 2$ that reads

$\begin{pmatrix}U_{t_{2g},t_{2g}}&U_{e_{g},t_{2g}}\\\
U_{t_{2g},e_{g}}&U_{e_{g},e_{g}}\end{pmatrix}\\\
=\begin{pmatrix}(\chi_{0})_{t_{2g},t_{2g}}&(\chi_{0})_{e_{g},t_{2g}}\\\
(\chi_{0})_{t_{2g},e_{g}}&(\chi_{0})_{e_{g},e_{g}}\end{pmatrix}^{-1}-\begin{pmatrix}(\chi)_{t_{2g},t_{2g}}&(\chi)_{e_{g},t_{2g}}\\\
(\chi)_{t_{2g},e_{g}}&(\chi)_{e_{g},e_{g}}\end{pmatrix}^{-1},$ (15)

where the indices represent the response in occupations of manifold $a$ to
perturbation of manifold $b$,
$(\chi)_{a,b}=\frac{\Delta\Lambda_{a}}{\Delta\alpha_{b}}$, for instance
$(\chi)_{e_{g},t_{2g}}=\frac{\Delta\Lambda_{e_{g}}}{\Delta\alpha_{t_{2g}}}$.
Pickett et al. showed that the orbital-resolved matrix elements of Eq. (15)
are related to the shell-averaged Hubbard parameter $U$ through a sum rule 27.
A slightly adapted form of this rule that accounts for the role of the non-
interacting response reads222Pickett et al. define the $d$ shell-averaged
Hubbard parameter (in their notation) as
$U_{dd}=\left[\sum_{a,b=t_{2g},eg}(\chi)_{a,b}\right]^{-1}$, whereas in this
work we use the definition of Eq. (16) due to the non-interacting (bare)
response.

$\displaystyle U$ $\displaystyle=$
$\displaystyle\left[\sum_{a,b=t_{2g},e_{g}}(\chi_{0})_{a,b}\right]^{-1}-\left[\sum_{a,b=t_{2g},e_{g}}(\chi)_{a,b}\right]^{-1}.$
(16)

The lesson to learn from Eq. (16) lies in the off-diagonal values of $\chi$,
whose physical implication consists in intra-shell screening (in this example
$e_{g}\leftrightarrow t_{2g}$). Since the sign of these off-diagonal values is
normally opposite to that of the diagonal ones 333This results from a physical
necessity: if a perturbation to the potential of a manifold $a$ causes the
occupation of $a$ to increase, the electrons must be borrowed from another
manifold, whose orbital occupation will naturally decrease, since the global
charge is kept constant in LR-cDFT., they diminish the contribution of the
term $\sum_{ab}(\chi)_{a,b}$, which may substantially increase $U$ values upon
computation of the inverse. Note that $\chi_{0}$ is not affected by this,
because the off-diagonal elements of the unscreened response are zero by
definition (except for numerical noise). In the extreme scenario where
perturbations are exclusively screened within the same shell, the sum over the
off-diagonal elements equals the trace of $\chi$, causing $\det(\chi)$ to
approach zero, leading to $U\to\infty$. Consequently, if the current
definition of the response matrices is used to compute orbital-resolved $U$
parameters, the determinant required for the full inversion of $\chi$
inadvertently reintroduces the shell-averaging of the response. This is
evident in the findings of the aforementioned study by Pickett et al., where
the computed $U_{t_{2g}}$ ($\equiv U_{t_{2g},t_{2g}}$) displayed only minimal
differences compared to $U_{e_{g}}$ ($\equiv U_{e_{g},e_{g}}$) and the shell-
averaged $U$ 27.

In order to derive on-site $U$ parameters that explicitly incorporate intra-
shell screening effects, it is necessary to set the off-diagonal matrix
elements of $\chi$ and $\chi_{0}$ to zero before computing their inverses.
This scheme was employed by Linscott et al., who refer to it as pointwise
(“$1\times 1$”) inversion, in order to compute screened and spin-resolved
Hubbard parameters for metal aquo complexes50. Solovyev et al. also applied
this strategy, albeit implicitly, to obtain screened $U$ parameters specific
to $t_{2g}$ in perovskites 79. Although the removal of off-diagonal elements
seems a drastic approximation, it allows for a more tailored definition of
Hubbard manifolds compared to the conventional DFT$+U$ approach. The latter
assumes (without proof) that ligand orbitals account for the majority of
screening while the role of intra-shell interactions is neglected. In fact,
several cRPA studies suggest that the opposite is true, showing that intra-
shell screening can be as significant, or even more so, as inter-shell
screening 86, 87. An important caveat to the LR-cDFT approach that also
affects the orbital-resolved form presented here is that unphysical $U$ values
may result when it is applied to fully occupied manifolds. This is a well-
known limitation 104, 105 that follows from the fact that the response of a
deep-lying state to a relatively small perturbation ($\approx 0.05\,$eV) is
often on the same order of magnitude as the numerical noise, leading to
instabilities during the inversion of the response matrices.

##  3 Computational details

All calculations are carried out using the Quantum ESPRESSO distribution 106,
107, 108. We have incorporated the capability to utilize and determine
orbital-resolved Hubbard $U$ parameters with the pw.x code and will make this
accessible in a future release. Structure and isosurface plots for pyrite
(FeS2) and pyrolusite (β-MnO2) are generated using VESTA109. Unless stated
otherwise, all systems are structurally optimized using the Broyden-Fletcher-
Goldfarb-Shanno (BFGS) algorithm 110 with convergence thresholds of
$10^{-4}\,$Ry, $10^{-3}\,$Ry/Bohr, and $0.5\,$kbar for the total energy,
forces, and pressure, respectively. KS wavefunctions (charge density) are
expanded in plane waves up to a kinetic-energy cutoff of 90 Ry (1080 Ry) using
PBE 111 pseudopotentials for pyrite and the Fe(II) molecular complexes, and
PBEsol pseudopotentials for β-MnO2 444 We used Fe.pbe-spn-kjpaw_psl.0.2.1.UPF
112, s_pbe_v1.4.uspp.F.UPF113, O.pbe-n-kjpaw_psl.0.1.UPF112, H.pbe-
rrkjus_psl.1.0.0.UPF114, P.pbe-n-rrkjus_psl.1.0.0.UPF114, C.pbe-n-
kjpaw_psl.1.0.0.UPF114, N.pbe-n-radius_5.UPF (THEOS pseudo) and
mn_pbesol_v1.5.uspp.F.UPF113 and O.pbesol-n-kjpaw_psl.0.1.UPF112 for β-MnO2.
taken from the SSSP Precision library v. 1.1.2 115, 116. The projected density
of states (PDOS) is obtained using a Gaussian smearing with a broadening
parameter of $0.02\,$Ry and employing the diag_basis feature of the projwfc.x
code, which projects the wavefunctions onto the eigenstates of the occupation
matrix rather than using unrotated atomic orbital projectors 54. This allows
for a clear distinction between $t_{g}$/$t_{2g}$ and $e_{g}$ states,
regardless of the orientation of the global coordinate system. Fe(II)
molecular complexes are simulated at a fixed $+2$ charge state in cubic boxes
with an edge length of $15\,$Å and using only the $\Gamma$ point to sample the
Brillouin zone. The total magnetization of the molecular complexes is always
fixed to either $4.0\,\mu_{\mathrm{B}}$ or $0.0\,\mu_{\mathrm{B}}$ in order to
compute high-spin (HS) and low-spin (LS) configurations, respectively. The
starting geometries of the Fe(II) molecular complexes are taken from the SI of
Ref. 76, whereas experimental structures are chosen as a starting points for
FeS2 and β-MnO2. The Brillouin zones of FeS2 and β-MnO2 are sampled with
uniform $\Gamma$-centered Monkhorst-Pack meshes of sizes $9\times 9\times 9$
and $4\times 4\times 6$, respectively.

We use Löwdin-orthogonalized atomic orbitals as Hubbard projectors (OAO) 95,
91 for all DFT+$U$ calculations including those with orbital-resolved $U$
parameters. Shell-averaged Hubbard parameters are evaluated using the DFPT
implementation 51, 52 of the HP code (hp.x) 67 included in Quantum ESPRESSO ,
employing $\mathbf{q}$ point meshes of size $2\times 2\times 2$ for FeS2 and
β-MnO2, and $1\times 1\times 1$ for the Fe(II) molecular complexes,
respectively. We compute orbital-resolved $U$ parameters according to the LR-
cDFT approach described in Section 2 by applying perturbations of
$\alpha=[-0.05,0.05]\,$eV to the manifold of interest and recording the non-
interacting and interacting responses of the orbital occupations. To avoid
interactions of perturbations with their periodic images, the calculations of
FeS2 and β-MnO2 are conducted in a $2\times 2\times 2$ supercell containing 96
atoms. We emphasize that DFPT and LR-cDFT are equivalent by construction, and
therefore yield the same Hubbard parameters when applied to identical systems
51. As the orbital response depends on the system’s electronic structure,
calculated Hubbard parameters may vary significantly when transitioning from a
PBE/PBEsol ground state to a PBE/PBEsol$+U$ one. Therefore, to achieve self-
consistency of the computed $U$ values, we employ an iterative procedure that
consists of structural optimizations and subsequent perturbative calculations
52, 30, 117. This procedure is repeated until the difference between the input
and output parameters falls below the predefined threshold of $\sim 0.1$ eV.
The resulting Hubbard parameters are reported in Section 4.

##  4 Results and discussion

###  4.1 Pyrite (FeS2)

####  4.1.1 Challenging theoretical description

Figure 1: The experimental crystal structure of pyrite (FeS2). Fe (brown) is
octrahedrally coordinated by S (yellow), which also forms characteristic S-S
dimers. The octahedra are slightly distorted and display angles different from
$90^{\circ}$.

Under normal conditions, pyrite is the stable polymorph of FeS2 and
crystallizes in the cubic space group Pa$\bar{3}$ with an experimental lattice
parameter $a=5.418\,$Å118. The crystal structure, shown in Figure 1, consists
of S2 dimers octahedrally coordinating Fe^2+ ions, which form a fcc sublattice
119. At 0 K, the compound is diamagnetic ($S=0$) due to the preferred LS
configuration of the Fe^2+ ions. Pyrite’s natural abundance, optical band gap
of $\sim$0.95 eV, and large optical absorption coefficient make it an
appealing material for photovoltaic applications. Nevertheless, despite
theoretical predictions suggesting an open circuit voltage of $\sim$ 0.71 V
based on the Shockley-Queisser equations, experimental results have
consistently fallen short, typically measuring values around 0.2 V 120. This
discrepancy is among several reasons motivating an accurate quantum-mechanical
description of pyrite’s ground state.

In early DFT studies, the LDA functional demonstrated exceptional accuracy in
predicting the equilibrium volume, band gap, and relative energy levels of the
S-$3p$ bands compared to the valence band maximum (VBM) 121, 119. PBE111 and
AM05122 provide a qualitatively similar picture but underestimate the band gap
by about 0.5 eV and 0.75 eV, respectively 120.

Figure 2: Stacked projected density of states (PDOS) of the nonmagnetic PBE
ground-state of pyrite FeS2. The small $t_{g}$ contributions555The point group
symmetry of Fe in pyrite is $T_{h}$ rather than $O_{h}$; thus, the triply
degenerate irreducible representation is correctly named $t_{g}$ instead of
$t_{2g}$123. in the conduction band and the $e_{g}$ contributions close to the
VBM are likely projection artifacts resulting from the slightly canted
geometry of the FeS6 octahedra.

At first glance, the dominant $e_{g}$ contributions in the conduction band
(Figure 5) indicate Mott-Hubbard insulation. However, a more detailed analysis
of the electronic band structure (Figure 3) reveals that the conduction band
minimum (CBM), located at the $\Gamma$ point, is composed of a $d$-$p$ σ^*
hybrid orbital with dominant S-$3p_{z}$ contributions.

Figure 3: (Fat-)Band structure of pyrite computed using the PBE functional.
The coloring indicates the most dominant contribution of the manifolds to the
individual bands. The uppermost valence bands are predominantly Fe-$t_{g}$,
whereas the conduction bands display a strong S-$3p_{z}$ character at $\Gamma$
but Fe-$e_{g}$ character elsewhere. Thus, the band gap of pyrite is not a
$d$-$d$ transition, but a $d$-$p$ one.

This observation holds the key to understanding why numerous electronic
structure methods, including more advanced approaches, face challenges in
improving the prediction of the band gap beyond the capabilities of PBE. For
instance, when employing $G_{0}W_{0}$, the small PBE band gap diminishes
further 124. In contrast, hybrid functionals such as PBE0125, HSE06126,
DSH127, and M06128 overestimate its value by more than $1\,$eV 129. Previous
PBE$+U$ studies focused on correcting the Fe-$d$ shell using shell-averaged
Hubbard parameters. For example, Sun et al. applied $U=2.0\,$eV (with PAW
projectors) 120, whereas Schena et al. used a combination of $U=3.0\,$eV with
$J=1.0\,$eV (Muffin-Tin sphere projectors)124. Notably, no set of ab initio
computed Hubbard parameters has been published for this material to date, to
the best of our knowledge. The rationale behind the use of empirical
parameters likely lies in a pronounced influence of shell-averaged $U$
corrections on the equilibrium properties, as we demonstrate hereinafter.

####  4.1.2 Double impact of shell-averaged U corrections

We perform PBE+$U$ calculations incorporating empirical on-site Hubbard $U$
corrections ranging from $1.0$ to $5.0\,$eV and analyze their impact on the
estimated band gap and on the crystal structure. To differentiate between a
band gap broadening caused by structural changes and one resulting directly
from the influence of $U$ on the electronic structure, we perform the
calculations under three different constrained conditions. Initially, we
retain the PBE ionic structure, keeping the cell vectors and ionic positions
fixed (setup 1). Subsequently, we perform another set of calculations with
ionic position relaxations while maintaining a constant cell volume (setup 2).
Lastly, we conduct full optimization, allowing for simultaneous adjustments of
both ionic positions and cell vectors (setup 3). While band gaps are
fundamentally outside the realm of DFT, DFT$+U$ should improve upon the
performance of the uncorrected functional as long as the compound’s frontier
states are well-represented by the Hubbard manifold57, 45. For this case, one
can argue that the Hubbard $U$ corrections (locally) act in the spirit of a
Koopmans-compliant functional17.

Figure 4 illustrates that even when no ionic relaxation is considered (setup
1), the band gap slowly but steadily expands with increasing $U$ values. The
shifts in band eigenvalues indicate that this expansion primarily stems from a
downshift of the $t_{g}$ orbitals’ KS potential in the valence region (see
Table SI 1). On the other hand, the CBM remains largely unaffected by the
Hubbard correction with the corresponding eigenvalues showing minimal changes.
Consistent with previous studies, we find that the experimental band gap is
accurately reproduced at $U\approx 2.0\,$eV, but is overestimated by
approximately 53% at $U=5.0\,$eV. This overestimation of the band gap at such
moderate values of $U$ is surprising since the CBM, consisting of S-$3p_{z}$
states, is not included in the Hubbard manifold, and its local deviation from
PWL remains uncorrected. Upon relaxing the ionic positions while keeping the
cell volume constant (setup 2), the band gap shows a more sensitive response
to higher values of $U$. The experimental band gap is already achieved at
$U\approx 1.5\,$eV and is overestimated by 115% at $U=5.0\,$eV. This trend
becomes even more pronounced when relaxation of the cell vectors is enabled
(setup 3), resulting in a band gap overestimation of 147% at $U=5.0\,$eV. As
shown in Figure 4(b), the additional expansion following ionic relaxations is
correlated with a significant contraction of the S-S bonds. This aligns well
with the findings of Eyert et al., who pointed out that the S-S bond length
governs the band gap, as it controls the dispersion of the lowest conduction
band around the $\Gamma$ point 121. Remarkably, in variable-cell calculations
(setup 3) the S-S bonds continue to contract with increasing values of $U$,
despite the simultaneous growth of the lattice parameter (inset of Figure
4(b)).

Figure 4: (a) calculated band gaps and (b) S-S bond lengths computed with
PBE$+U$ using three setups: at fixed cell and ionic positions (setup 1); at
constant cell volume relaxing the ionic positions (setup 2); and relaxing both
the ionic positions and the cell volume (setup 3). The inset in the bottom
panel shows the dependence of the lattice parameter on $U$ during volume
relaxations. Red dashed horizontal lines correspond to the experimental
values118, 130.

The differences observed between relaxed and fixed structure calculations
reveal that shell-averaged $U$ corrections modify the band gap through two
distinct mechanisms: first, through the expected downshift of the valence band
eigenvalues, and second by causing pronounced contractions of the S-S bonds.

We now shift our focus from conventional equilibrium properties to investigate
the underlying reasons behind the substantial influence of shell-averaged $U$
corrections on these bonds, which can be elucidated by examining the
occupations of the individual orbitals. Table 1 presents the eigenvalues of
the Fe-$d$ occupation matrix and the corresponding Hubbard energies $E_{U}$
for an illustrative case with $U=3.0\,$eV, both prior to and after complete
structural relaxation. As the eigenvalues of the $t_{g}$ orbitals ($d_{xy}$,
$d_{xz}$, $d_{yz}$) approach idempotency, this manifold is responsible for the
smaller share of the overall Hubbard energy, approximately 33%. In contrast,
the $e_{g}$ orbitals ($d_{x^{2}-y^{2}}$ and $d_{z^{2}}$) contribute to nearly
66% of the total Hubbard energy due to possessing occupation eigenvalues far
from either 0 or 1. Structural relaxation induces both intra-shell and inter-
shell charge transfers that allow for a 6% reduction of the total Hubbard
energy. At the intra-shell level, the occupation of $t_{g}$ orbitals grows at
the expense of $e_{g}$. Even more significantly, the overall $d$-occupancy
slightly drops from 6.996 to 6.970, indicating a migration of some $e_{g}$
electrons into adjacent S-$p$ orbitals. This transfer of charge effectively
enhances the S-S bond order and leads to the observed contraction of
$d_{\mathrm{S-S}}$. Thus, the unexpectedly potent influence of the value of
$U$ on the band gap is rooted in the correction of the $e_{g}$ manifold.

Table 1: Occupation eigenvalues $\lambda$ and corresponding Hubbard energies $E_{U}$ obtained for the Fe-$d$ orbitals in FeS2 using PBE$+U$ with $U=3.0\,$eV. Results are shown for the fixed PBE structure (setup 1, ${a=5.406}$ Å and a S-S distance of 2.195 Å, top) and after structural relaxation (setup 3, $a=5.469$ Å and a S-S distance of 2.132 Å, bottom). The last column shows the total occupations ($2\sum\lambda$) and Hubbard energies ($2\sum E_{U}$), which correspond to twice the sum (due to spin degeneracy) of the individual contributions. | Orbital | $d_{x^{2}-y^{2}}$ | $d_{z^{2}}$ | $d_{xy}$ | $d_{xz}$ | $d_{yz}$ | $2\Sigma$
---|---|---|---|---|---|---|---
setup 1 | $\lambda$ | 0.383 | 0.383 | 0.901 | 0.901 | 0.930 | 6.996
|  $E_{U}$ (meV) | 354 | 354 | 134 | 134 | 98 | 2148
setup 3 | $\lambda$ | 0.359 | 0.359 | 0.915 | 0.915 | 0.937 | 6.970
|  $E_{U}$ (meV) | 343 | 343 | 117 | 117 | 89 | 2018

We recall that the eigenvalues listed in Table 1 are not universal quantities
but the result of a projection of the KS wave functions onto atom-centered
OAO. The accuracy of on-site occupations provided by atomic-like projectors
relies on the similarity between a given orbital’s shape in the system under
inspection and that in a free atom, since the general shape and extension of
projector AOs is typically determined in free atom calculations. However, one-
center projections using atomic orbitals struggle to account for strong
orbital hybridization, where electrons localize ‘off-site’ between the bonded
atoms. In the case of FeS2, the $t_{g}$ manifold is only marginally bonding,
while the formally unoccupied $e_{g}$ orbitals form $\sigma$ MOs with
neighboring S-$3p_{z}$ orbitals. Consequently, it is misleading to interpret
the eigenvalues corresponding to $d_{x^{2}-y^{2}}$ and $d_{z^{2}}$ as
indicative of actual on-site orbital occupancies in the context of Hubbard $U$
corrections. Moreover, their numerical values significantly hinge on
computational factors like the chosen pseudopotentials and the specific charge
state for which they were parameterized 121, 56.

Within the setup of this study, the $U$ correction effectively penalizes the
hybridization between Fe and S; however, the use of different Hubbard
projector functions cause the opposite effect, for example, if NAO projectors
yield $e_{g}$ eigenvalues larger than $0.5$. In this case, the Hubbard
correction would draw electrons into the $e_{g}$ manifold rather than
expelling them (see Eq. (3)).

####  4.1.3 Absence of intra-shell screening in shell-averaged LR-cDFT

Having investigated how the shell-averaged $U$ parameter affects equilibrium
observables, our focus now shifts to understanding the implications when this
parameter is determined from first principles. Applying DFPT to the PBE ground
state of FeS2 yields a $U$ value of $7.37\,$eV. After undergoing four
iterations within the self-consistency loop detailed in Section 2.3, $U$
converges to $6.47\,$eV. Unfortunately, neither of these parameters can
reasonably reproduce the experimental characteristics of pyrite. In fact, both
values result in the stabilization of a spurious ferromagnetic ground state
($2\,\mu_{B}/$cell) in unrestricted open-shell calculations. For comparison,
all of the following $U$ parameters and observables are reported for the
nonmagnetic ground state that was enforced by fixing the total magnetization
to 0.0$\,\mu_{\mathrm{B}}$.

Figure 5: Unscreened/non-interacting (top) and screened/interacting (bottom)
responses of the occupations $\lambda$ of manifolds “$\mathbf{a}$” to
perturbations $\alpha$ applied to manifolds “$\mathbf{b}$” in FeS2, obtained
through the orbital-resolved LR-cDFT approach. The starting point is the fully
relaxed PBE ground state.

The reason for this significant overestimation of the shell-averaged Hubbard
$U$ parameter can be understood from Figure 5, which shows the response of
orbital occupancies to perturbations for both the entire $3d$-shell and for
its irreducible representations $t_{g}$ and $e_{g}$. This information cannot
be extracted from shell-averaged LR-cDFT, but can be recovered by adopting the
orbital-resolved approach. The substantial opposing responses of $e_{g}$
occupations to perturbations of $t_{g}$ and vice versa (depicted by dashed
lines in the bottom panel) suggest the presence of a robust intra-shell
$t_{g}\leftrightarrow e_{g}$ screening channel. Inserting the values of
$\chi_{0}$ and $\chi$ (represented by the line slopes in Figure 5) into Eq.
(16), one obtains

$\displaystyle U$ $\displaystyle=$
$\displaystyle\left[\sum\begin{pmatrix}-0.4608&0.0001\\\
0.0001&-0.4131\end{pmatrix}\right]^{-1}-\left[\sum\begin{pmatrix}-0.2266&0.1843\\\
0.1843&-0.2592\end{pmatrix}\right]^{-1}$ (17) $\displaystyle=$
$\displaystyle(-0.8737)^{-1}-(-0.1172)^{-1}=-1.144-(-8.532)=7.388\,\mathrm{eV},$

where the first and second matrix represent the unscreened and the screened
responses, respectively. Note the almost perfect agreement between this LR-
cDFT Hubbard parameter ($U=7.39\,$eV) and the one derived from DFPT
($U=7.37\,$eV); the minor discrepancies are due to numerical noise. The
individual matrix elements of Eq. (17) reveal that the apparent overestimation
of the shell-averaged Hubbard parameter primarily stems from the suppression
of intra-shell screening that occurs due to the simultaneous perturbation of
the $t_{g}$ and the $e_{g}$ manifolds. While the individual screened responses
of the $t_{g}$ and $e_{g}$ orbitals are substantial
($\chi_{t_{g},t_{g}}=-0.2266\,\mathrm{eV}^{-1}$ and
$\chi_{e_{g},e_{g}}=-0.2592\,\mathrm{eV}^{-1}$), the overall response of the
$d$ shell is massively reduced by the off-diagonal elements of $\chi$ (both
$0.1843\,\mathrm{eV}^{-1}$). Inversion of this apparently small screened
response of the $d$ shell ($-0.1172\,\mathrm{eV}^{-1}$) yields
$\chi^{-1}=-8.532\,$eV, which is what ultimately fuels the overestimation of
the shell-averaged Hubbard parameter.

It should be noted that this shell-averaged $U$ parameter also encompasses, to
a certain extent, the response of the S-$3p$ orbitals, because the $e_{g}$
eigenstates likely contain substantial S-$3p$ contributions. With these
insights in mind, it becomes pertinent to investigate whether the orbital-
resolved approach to Hubbard $U$ can offer improvements over the
unsatisfactory performance of shell-averaged PBE$+U$.

####  4.1.4 Application of the orbital-resolved $U$

Employing the orbital-resolved LR-cDFT methodology detailed in Section 2.3, we
proceed to compute individual Hubbard parameters for the $t_{g}$ and $e_{g}$
orbitals of FeS2. Prior to inverting the response matrix, all off-diagonal
intra-shell matrix elements are set to zero:

$(\chi)^{II}_{i\neq j}=0\textnormal{ , and }(\chi_{0})^{II}_{i\neq j}=0.$ (18)

This ensures that the resulting orbital-resolved Hubbard parameters
incorporate the effect of $t_{g}\leftrightarrow e_{g}$ screening. The
parameters $U_{t_{g}}$ and $U_{e_{g}}$ converge rapidly, reaching their self-
consistent values of $3.29\,$eV and $2.16\,$eV, respectively, within just
three iterations of self-consistency. The final response matrices read

$\displaystyle\begin{pmatrix}U_{t_{g}}&0\\\ 0&U_{e_{g}}\end{pmatrix}$
$\displaystyle=$ $\displaystyle\left[\begin{pmatrix}-0.3516&0\\\
0&-0.1948\end{pmatrix}\right]^{-1}-\left[\begin{pmatrix}-0.1998&0\\\
0&-0.1187\end{pmatrix}\right]^{-1}$ (19) $\displaystyle=$
$\displaystyle\begin{pmatrix}3.29&0\\\ 0&2.16\end{pmatrix}.$

Following the approach of Solovyev et al., one may also exclusively target the
$t_{g}$ manifold and not correct the $e_{g}$ states at all. This choice is
driven by the expectation that interactions of genuine on-site character
should manifest within the occupied $t_{g}$ orbitals rather than in the
hybridized and formally empty $e_{g}$ manifold. Moreover, as previously
stated, it is unlikely that one-center atomic orbital projectors are suited to
provide meaningful on-site occupation numbers for the $e_{g}$ orbitals. With
$U_{e_{g}}$ corrections absent, the converged value of $U_{t_{g}}$ slightly
decreases to $3.01\,$eV.

Table 2: Comparison of the equilibrium lattice parameter $a$, the S-S bond length $d_{\mathrm{S-S}}$ and the band gap $E_{\mathrm{g}}$ of pyrite derived from Hubbard-corrected PBE calculations targeting different manifolds. On-site Hubbard $U_{1}$ and $U_{2}$ refer to various parametrizations of the PBE+$U$ approach, while $V$ is the inter-site Hubbard parameter of the PBE+$U$+$V$ approach. | PBE | $+U$ | ${+U\\!+\\!V}$ | $+U^{\mathrm{emp}}$ | $+U_{t_{g}}+U_{e_{g}}$ | $+U_{t_{g}}$ | Expt.
---|---|---|---|---|---|---|---
$U_{1}\,$(eV) | — | 6.47 | 6.73 | 1.50 | 3.29 | 3.01 | —
$U_{2}$ or $V\,$(eV) | — | — | 0.62 | — | 2.16 | — | —
$a\,$(Å) | 5.41 | 5.63 | 5.51 | 5.43 | 5.46 | 5.44 | 5.42118
$d_{\mathrm{S-S}}\,$(Å) | 2.20 | 2.10 | 2.12 | 2.16 | 2.13 | 2.14 | 2.16118
$E_{\mathrm{g}}\,$(eV) | 0.41 | 2.63 | 2.06 | 1.00 | 1.70 | 1.45 | 0.95130

Having established the necessary parameters, we now focus on their impact in
practical calculations. Table 2 shows that the orbital-resolved approach
consistently outperforms the shell-averaged approach, except when $U$ is tuned
empirically ($U^{\mathrm{emp}}$). When utilizing only $U_{t_{g}}$, both the
lattice parameter and the S-S bond length exhibit deviations of less than 1%
from the experimental values. With orbital-resolved corrections applied to
${t_{g}}$ and $e_{g}$, this agreement with the experimental data is slightly
worsened but still surpasses that of shell-averaged PBE$+U$, where the lattice
parameter deviates by 4%. Furthermore, the significant overestimation of the
band gap amounting to 177% in the case of shell-averaged PBE$+U$ is
substantially mitigated to 53% when switching to $U_{t_{g}}$. The
DFT${+U\\!+\\!V}$ approach yields noticeable improvements compared to the
traditional $+U$ approach; however, it appears that the inter-site $V$ term
lacks the strength to fully restore the hybridization between Fe and S that is
suppressed by the shell-averaged on-site term. Therefore, the estimation of
the lattice parameters and the band gap still remains notably worse than for
the orbital-resolved PBE$+U$ approach.

These results underscore the critical importance of the choice of the Hubbard
manifold in DFT$+U$ calculations for the FeS2 polymorph pyrite. Under the
shell-averaged approximation, the ground-state properties are extremely
sensitive to the numerical value of $U$. This sensitivity primarily arises due
to the correction of the hybridized $e_{g}$ part of the $d$ shell. Moreover,
the shell-averaged Hubbard parameter derived from LR-cDFT lacks full
screening, as the intra-shell screening channel ($t_{g}\leftrightarrow e_{g}$)
is deactivated during the concurrent perturbation of both manifolds. On the
other hand, the orbital-resolved approach accounts for intra-shell screening
during perturbative calculations, resulting in significantly smaller Hubbard
parameters. We have shown that correcting the localized Fe-$t_{g}$ manifold
alone suffices to obtain structural properties in very good agreement with
experimental data, all without requiring empirical adjustments to the
parameters. A viable option for further improvements involves supplementing
the orbital-resolved $U$ corrections with inter-site $V$ terms, as discussed
in Section 5. Additionally, the PWL condition should also be enforced for the
S-$3p_{z}$ orbital composing the CBM. Given the strong $sp^{3}$ hybridization
of S, this would, however, require an adaption of the atomic-like projectors
used in this study.

###  4.2 Fe(II) molecular complexes

####  4.2.1 Adiabatic spin energy differences

A crucial quantity for many TM compounds is the total energy difference
between their HS and LS states, denoted as:

$\Delta E_{H-L}=E^{\mathrm{HS}}_{tot}-E^{\mathrm{LS}}_{tot}\textnormal{ .}$
(20)

Based on calculations of six Fe(II)-hexacomplexes, Mariano et al. showed that
shell-averaged Hubbard $U$ corrections to PBE and LDA introduce an unphysical
bias against LS configurations, primarily due to substantial discrepancies in
the Hubbard energies of LS and HS complexes 76. This arises from the electron
occupation patterns: LS-Fe(II) complexes exhibit a $t_{2g}^{6}/e_{g}^{0}$
electron configuration, while HS complexes a $t_{2g}^{4}/e_{g}^{2}$
configuration. In octahedral coordination environments, unoccupied $e_{g}$
orbitals display substantial $\sigma$ overlap with neighboring ligand
orbitals. Consequently, in LS-Fe(II) all four $e_{g}$ spin-orbitals can adopt
fractional occupation eigenvalues, whereas in HS-Fe(II) this only applies to
two spin-orbitals, as the other two are fully occupied. Considering that the
Hubbard energy functional (Eq. (1)) induces corrections of up to $U/2\,$eV for
fractional occupations, but of $0\,$eV for idempotent occupations, Hubbard
energies of LS-Fe(II) tend to be larger than those of HS-Fe(II). Furthermore,
the LR-cDFT approach was shown to generate larger shell-averaged $U$
parameters for LS-Fe(II) compared to HS-Fe(II) 76, thus enhancing this
disparity.

In the following, we investigate whether the orbital-resolved DFT$+U$ approach
can rectify the unphysical bias against LS states and deliver accurate
adiabatic spin energy differences. For this purpose, we compute $\Delta
E_{H-L}$ for all six octahedrally coordinated Fe(II)-hexacomplexes studied in
Ref. 76 and benchmark the results against coupled-cluster corrected CASPT2
(CASPT2/CC)131 values, which are also taken from Ref. 76. We compare the
performance of two distinct Hubbard manifolds. In manifold I, Hubbard
corrections are applied either only to $e_{g}$ states ([Fe(H2O)_6]^2+ and
[Fe(NH3)_6]^2+) or to only $t_{2g}$ states ([Fe(NCH)_6]^2+, [Fe(PH3)_6]^2+,
[Fe(CO)_6]^2+, [Fe(CNH)_6]^2+). While it would be desirable to obtain
$t_{2g}$-specific Hubbard parameters for the former two compounds, this proves
elusive due to the fact that the $t_{2g}$ orbitals are fully occupied in the
numerical sense ($\lambda\geq 0.997$) and thus respond non-linearly to
perturbations 104. For the same reason, no $U_{t_{2g}}^{\mathrm{HS}}$
parameters could be computed for the HS complexes either. Therefore, all
Hubbard corrections to $t_{2g}$ states utilize the $U_{t_{2g}}^{\mathrm{LS}}$
value. In manifold II, both $t_{2g}$ and $e_{g}$ orbitals are corrected
simultaneously. When comparing the present work with Ref. 76, it is important
to note that we use OAO projector functions and self-consistent procedures for
geometry relaxations and determination of Hubbard terms. A comprehensive list
of all complexes, the Hubbard manifolds considered and the corresponding
orbital-resolved and shell-averaged $U$ parameters is provided in Table 3.

Table 3: PBE+$U$ manifolds and their calculated on-site parameters in eV. For $t_{2g}$ the same parameter was applied to LS and HS complexes. The complexes are sorted in ascending order according to the ligand’s field strength. Complex | manifold I | manifold II | $U_{t_{2g}}^{\mathrm{LS}}$ | $U^{\mathrm{LS}}$ | $U^{\mathrm{HS}}_{e_{g}}$ | $U^{\mathrm{LS}}$ | $U^{\mathrm{HS}}$
---|---|---|---|---|---|---|---
[Fe(H2O)_6]^2+ | $e_{g}$ | — | — | 4.28 | 2.69 | 5.71 | 4.13
[Fe(NH3)_6]^2+ | $e_{g}$ | — | — | 3.64 | 3.07 | 6.14 | 4.43
[Fe(NCH)_6]^2+ | $t_{2g}$ |  $t_{2g}$+$e_{g}$ | 7.05 | 2.12 | 1.87 | 6.76 | 5.53
[Fe(PH3)_6]^2+ | $t_{2g}$ |  $t_{2g}$+$e_{g}$ | 4.13 | 2.29 | 2.12 | 6.88 | 4.78
[Fe(CO)_6]^2+ | $t_{2g}$ |  $t_{2g}$+$e_{g}$ | 4.63 | 1.92 | 1.89 | 7.16 | 5.43
[Fe(CNH)_6]^2+ | $t_{2g}$ |  $t_{2g}$+$e_{g}$ | 4.64 | 2.13 | 1.72 | 7.43 | 5.79

The resulting adiabatic spin energy differences are visualized in Figure 6.
Additionally, Table 4 lists the mean absolute error (MAE) of the different
corrections against the CASPT2/CC131 reference values. In line with the
findings of Mariano et al., shell-averaged PBE$+U$ drastically overstabilizes
HS configurations and fails to reproduce the trend of increasing $\Delta
E_{H-L}$. Conversely, uncorrected PBE displays its known bias towards LS
configurations. The introduction of orbital-resolved $U$ corrections
substantially enhances the predictive accuracy across all complexes analyzed,
outperforming both PBE and the shell-averaged PBE$+U$ by a large margin.
Correcting either $t_{2g}$ or $e_{g}$ (manifold I) closely aligns adiabatic
spin energy differences with the CASPT2/CC values, resulting in the lowest MAE
of $0.27\,$eV. The simultaneous correction of $t_{2g}$ and $e_{g}$ (manifold
II) yields a slightly higher MAE of $0.57\,$eV, with larger deviation
particularly for [Fe(PH3)6]^2+ and [Fe(CNH)6]^2+, but a better agreement for
[Fe(NCH)_6]^2+. In general, $U_{e_{g}}$ corrections are smaller and less
significant for the overall quantitative agreement than $U_{t_{2g}}$
corrections, since the $t_{2g}$ orbitals host more – and more localized –
electrons. The good overall agreement (qualitative and quantitative) of
orbital-resolved PBE$+U$ with the CASPT2/CC data is surprising since the
latter is a multiconfigurational approach, while the former can only improve
the accuracy of a single reference state.

Figure 6: Adiabatic spin energies calculated for six Fe(II)-hexacomplexes using different DFT corrections and CASPT2/CC as the reference method. Data for PBE, PBE+$U$, and CASPT2/CC was taken from the SI of Ref. 76. Table 4: Mean absolute error (MAE) of $\Delta E_{H-L}$ for various Hubbard correction manifolds against the CASPT2/CC reference value from Ref. 76. | PBE76 | +$U$ | +$U_{t_{2g}}$ or +$U_{e_{g}}$ | +$U_{t_{2g}}$+$U_{e_{g}}$
---|---|---|---|---
MAE (eV) | 1.07 | 1.94 | 0.27 | 0.57

The good performance of the orbital-resolved $U$ compared to the shell-
averaged approach stems from two main factors. First, excluding the $e_{g}$
orbitals from the Hubbard manifold of strong-field complexes (manifold I)
eliminates the primary cause of spuriously large Hubbard energies. Again, the
exclusion of $e_{g}$ orbitals can be justified with the inadequate
representation of their strongly hybridized character by atomic orbital
projectors. Second, even when the $e_{g}$ orbitals are included (manifold II),
independently computed $U_{e_{g}}$ parameters are significantly smaller than
$U_{t_{2g}}$ or shell-averaged $U$ parameters owing to the explicit
incorporation of intra-shell screening. For instance, in [Fe(CO)_6]^2+
$U_{e_{g}}^{\mathrm{LS}}=1.92\,$eV compared to $U^{\mathrm{LS}}=7.16\,$eV.
Hence, the orbital-resolved approach effectively diminishes the impact of non-
ideal Hubbard projectors.

####  4.2.2 Piece-wise linearity of the total energy

As mentioned earlier, the use of on-site Hubbard $U$ corrections in
(semi-)local DFT is motivated by the mitigation of SIE, which have been linked
to the spurious global deviation from PWL of the total energy with respect to
fractional addition or removal of electronic charge $q$ to the entire system
12, 15, 132. However, a study by Zhao et al. indicates that (shell-averaged)
$U$ corrections to TM d shells might be unfit for this purpose for a wide
range of TM/ligand combinations when the Hubbard $U$ parameters are derived
from LR-cDFT15.

Thus, to directly benchmark the performance of orbital-resolved DFT+$U$, we
explicitly determine the global curvature following the approach presented in
Ref. 15. For our analysis, we chose the strongest ligand complex
[Fe(CNH)_6]^q+ and perform several fixed-charge calculations where the total
charge is varied between $q=2$ and $q=3$ in increments of $0.1$e-.

Figure 7: Deviation from piecewise linearity of the total energy upon
fractional addition/removal of charge of [Fe(CNH)6]^(q+) in low-spin (left)
and high-spin (right) configurations. Note that the data points of
PBE$+U_{t_{2g}}$ and PBE$+U_{t_{2g}}+U_{e_{g}}$ overlap in the LS case.

A useful metric to assess a functional’s deviation from PWL is given by
$E_{\mathrm{dev}}$, which is computed by subtracting the DFT total energy of a
fractional charge calculation from the linear interpolation between the
energies of the integer-charge endpoints $E_{tot}(q=2)$ and $E_{tot}(q=3)$.
Note that this explicit approach differs from the method employed by Mariano
et al., who approximate $E_{\mathrm{dev}}$ using a cubic interpolation
parameterized by the energy difference $\Delta
E_{q}=E_{tot}(q=2)-E_{tot}(q=3)$ and the HOMO and LUMO eigenvalues of the
integer-charge system 76.

Figure 7 shows that while shell-averaged Hubbard $U$ corrections markedly
reduce the deviation from PWL of HS [Fe(CNH)_6]^q+, the global curvature of
the LS compound is not eased. In fact, in the LS case the typical convex-
shaped FCE displayed by bare PBE is transformed into a concave one of similar
magnitude. In contrast, the orbital-resolved corrections reduce the deviation
from PWL by an order of magnitude for both the LS and the HS configurations.
It is worth noting that the remaining curvature barely differs between
manifold I (+$U_{t_{2g}}$ only) and manifold II (+$U_{t_{2g}}$  +
$U_{e_{g}}$). Again, this observation is likely related to the low numerical
value of the $U_{e_{g}}$ parameters, which further supports the previous
assumption that the primary contributor to the FCE is the rather localized
$t_{2g}$ manifold. For the same reason, applying $U$ corrections to ligand
orbitals like C-$p$ is unlikely to significantly affect the FCE.

###  4.3 The importance of correcting ligand states in β-MnO2

Unlike in TM sulfides, partially filled d shells of TM oxides are a frequent
target for on-site corrections (see the discussion e.g. in Ref. 73). While
Mott-Hubbard insulators such as the prototypical wüstite (FeO) present an
energy gap separating different $d$ states, so-called charge-transfer
insulators, often involving highly oxidized and/or heavy TM species, exhibit
O-$2p\rightarrow d$ band gaps 133. The significant fraction of electrons
localized on O-$2p$ orbitals adds complexity to the definition of the Hubbard
manifold for DFT$+U$ calculations, as studied below for the case of β-MnO2.

Crystallizing in the rutile structure (space group P$4_{2}$/mnm), β-MnO2
displays a complex helical (screw-type spiral) antiferromagnetic (AFM)
ordering below $T_{N}=92\,$K 134. As done here, this ordering can be
approximated by a collinear arrangement of the Mn^4+ ions along [001],
referred to as A1-AFM69. Considerable uncertainty surrounds the electronic
properties of β-MnO2. Sato et al. measured a large electrical resistivity at
$0.3\,$K, suggesting insulating behavior 134. However, the magnitude of the
band gap has not been yet determined with universally accepted accuracy. Mid-
twentieth-century works reported narrow band gaps around $0.26-0.28\,$eV 135,
136, while more recent optical absorption studies reported much larger values
ranging from $1.5\,$eV (nanocacti and nanorods) 137 to $2.0\,$eV for β-MnO2
nanostructures grown on fluorine-doped tin oxide 138. These results align with
a reported hybrid functional DFT (PBE0) value of $1.5\,$eV, but not with a
prediction of the closely related HSE03 functional ($0.6\,$eV) 139. Previous
works employing the shell-averaged DFT+$U$ approach suggest that a band gap is
not opened unless the Mn-$d$ states are corrected using $U$ values larger than
6 eV139, 54. In general, the calculated properties are sensitive to the choice
of projector functions. For instance, NAO projectors stabilize a ferromagnetic
ordering, whereas OAO projectors favor the expected A1-AFM ordering 54. Only
by adding inter-site $V$ between Mn-$d$ and O-$p$ states 54 or Hund’s $J$
corrections 140, small gaps ranging from $0.25$ to $0.32\,$eV emerge. A band
gap of $0.8\,$eV was also obtained 83 with the anisotropic DFT+$U$+$J$
approach of Czyżyk and Sawatzky. In this method, shell-averaged Hubbard $U$
and $J$ corrections are augmented with orbital-resolved $U_{mm^{\prime}}$ and
$J_{mm^{\prime}}$ matrices parameterized from summation relations involving
Slater integrals of atomic Hartree-Fock calculations 80. This approach thus
differs from the here-presented scheme, which does not rely on model systems
like the free atom and includes intra-manifold screening effects, such as
$t_{2g}\leftrightarrow e_{g}$.

The wide range of band gap values reported in both the theoretical and
experimental literature underscores the pressing need for a more profound
understanding of the insulating behavior of β-MnO2. Considering the charge-
transfer nature of the band gap and the coexistence of relatively localized
and strongly hybridized states in both the Mn-$d$ shell and the O-$p$ shell,
orbital-resolved DFT+$U$ calculations offer a promising starting point. In the
interest of comparability, our calculations are carried out using a setup
equivalent to that of Ref. 54, imposing the A1-AFM ordering in a $2\times
2\times 2$ supercell containing 48 atoms and considering once again different
Hubbard manifolds.

Table 5: Comparison of equilibrium properties of A1-AFM β-MnO2 obtained from PBEsol without and with various Hubbard corrections. $U_{1}$ and $U_{2}$ refer to various parameterizations of the PBEsol+$U$ approach, while $V$ is the inter-site Hubbard parameter of the PBEsol+$U$+$V$ approach. $a$ and $c$ are the lattice parameters, $V_{0}$ is the unit cell volume, and $E_{g}$ is the band gap value. All the presented results (including those from Ref. 54) were obtained using OAO projectors.  | PBEsol | +$U^{\mathrm{Mn}}_{t_{2g}}$ |  $+U^{\mathrm{Mn}}_{t_{2g}}+U^{\mathrm{O}}_{p_{z}}$ |  $+U^{\mathrm{Mn}}$ (Ref. 54) |  ${+U^{\mathrm{Mn}}\\!+\\!V^{\mathrm{Mn-O}}}$ (Ref. 54) | Expt.
---|---|---|---|---|---|---
$U_{1}\,$ (eV) | — | 1.59 | 1.64 | 6.34 | 6.76 | —
$U_{2}/V\,$ (eV) | — | — | 4.62 | — | (0.99, 1.10) | —
$a\,$ (Å) | 4.37 | 4.37 | 4.38 | 4.40 | 4.39 | 4.40 141
$c\,$ (Å) | 2.83 | 2.86 | 2.88 | 2.94 | 2.92 | 2.88 141
$V_{0}\,$ (Å3) | 54.68 | 54.64 | 55.26 | 57.07 | 56.35 | 55.79 141
$E_{\mathrm{g}}\,$ (eV) | 0.00 | 0.00 | 1.01 | 0.02 | 0.32 |  0.26135, 1.50 137, 2.0138

The first manifold solely considers the ${t_{2g}}$ orbitals in order to single
out the impact of the localized Mn-$d$ states. 666Here, we use the $t_{2g}$
irreducible representation in spite of the actual $D_{4h}$ site symmetry of
Mn^4+ and therefore neglect the small energetic splitting between the
$A_{1g}$, $B_{1g}$ and $B_{2g}$ representations. As shown in Table 5, this
correction only marginally influences the band gap or the structural
properties in comparison to bare PBEsol. This observation aligns with the
modest value of $U^{\mathrm{Mn}}_{t_{2g}}=1.59\,$eV, indicative of minimal
local curvature due to the ${t_{2g}}$ states. A much more significant
improvement is achieved by considering the localized O-${p_{z}}$ orbitals
(Figure 8). Correcting these states with a self-consistent Hubbard parameter
$U^{\mathrm{O}}_{p_{z}}=4.62\,$eV (in addition to the
$U^{\mathrm{Mn}}_{t_{2g}}$ correction) leads to an insulating gap of
$1.01\,$eV. Moreover, the lattice parameters approach their experimental
values, resulting in an error on the predicted cell volume smaller than 1%.
Conversely, a shell-averaged $U$ correction as large as $6.34\,$eV for Mn-$d$
states only results in less accurate cell parameters and a still negligible
band gap. This situation is only partly mended by the application of an
additional inter-site $V$ parameter 54.

Figure 8: Isosurface plot of the highest fully occupied KS band of A1-AFM
β-MnO2 obtained using PBEsol$+U^{\mathrm{Mn}}_{t_{2g}}$, showing $\pi$
interactions between O-${p_{z}}$ and Mn-${t_{2g}}$ orbitals. Mn atoms of
opposite spin-polarization are colored dark and light brown, O atoms are red.
For clarity, the interaction is only shown for one MnO6 octahedron.

The importance of including O-${p_{z}}$ in the Hubbard manifold is evidenced
by the PDOS depicted in Figure 9. The states near the Fermi level exhibit
notable O-${p}$ contributions, regardless of the correction applied. These
contributions primarily stem from localized ${p_{z}}$ orbitals, which engage
in $\pi$ interactions with the Mn-${t_{2g}}$ manifold (Figure 8). The other
two $p$ orbitals are subject to $sp^{2}$ hybridization and contribute to the
formation of the lower-lying $\sigma$ MOs in conjunction with Mn-${e_{g}}$
states. This observation underscores that applying shell-averaged $U$
corrections to the entire O-$p$ shell would weaken the $\sigma$ bonds and lead
to Mn-O underbinding. For the same reason, applying corrections to the
Mn-${e_{g}}$ states is neither required nor useful.

Figure 9: Stacked spin-resolved PDOS of the A1-AFM collinear ordering of
β-MnO2 computed with PBEsol and two types of orbital-resolved $U$ corrections.
Metallic behavior is observed when using PBEsol and
PBEsol$+U^{\mathrm{Mn}}_{t_{2g}}$, whereas an insulating character is obtained
with PBEsol$+U^{\mathrm{Mn}}_{t_{2g}}+U^{\mathrm{O}}_{p_{z}}$. On each panel,
the upper and lower subpannels correspond to the spin-up and spin-down
channels, respectively.

In conclusion, the orbital-resolved approach pinpoints deviations from PWL in
the localized O-$p_{z}$ orbitals as the potential root-cause of the absence of
a band gap in β-MnO2 when using DFT with semi-local functionals. This outcome
underscores the frequently overlooked significance of ligand orbitals in
electronic localization and establishes a foundation for future
investigations. Again, achieving more comprehensive descriptions of this
intricate charge-transfer insulator could entail incorporating the orbital-
resolved $U$ approach alongside inter-site $V$ and unlike-spin Hund’s $J$
terms.

##  5 Outlook: orbital-resolved inter-site Hubbard interactions?

Hubbard ${U}$ corrections are primarily designed to mitigate the FCE in KS
states exhibiting strong atomic-like character, that is, at the on-site level.
Orbitals involved in strong inter-atomic hybridization should not be corrected
using $U$ terms. This limitation is a consequence of the one-center projectors
that are typically used to determine the occupation numbers for the Hubbard
$U$ energy functional, which inherently lack the ability to adequately
represent localization on hybrid (i.e., molecular) orbitals. To counteract the
FCE within covalent environments, the extended DFT+$U\\!+\\!V$ framework was
introduced in Ref. 53, where conventional shell-averaged $U$ corrections are
augmented with an inter-site term scaled by a Hubbard parameter $V$ in a way
that restores ligand hybridization through two-center (dual) occupation
numbers. However, although the +$V$ correction often improves the predictions
of electronic and structural properties compared to plain DFT+$U$, the results
presented for FeS2 and MnO2 exemplify that it only partially mitigates the
impact of shell-averaged $U$ corrections. Specifically, inter-site terms
cannot rectify a substantial bias in the total energy arising from on-site
terms since the values of $U$ (for first-row TM elements) are typically
between 4 and 10 eV, which is several times greater than typical values of
$V$, which amount to $\approx 1\,$eV. For instance, the Hubbard $U$ and $V$
parameters of the strong-field molecular complex [Fe(CNH)_6]^2+ from Sec. 4
are $8.42$ (for Fe-$3d$ states) and $1.12\,$eV (between Fe-$3d$ and C-$2p$
states) for the LS configuration, and $6.57$ and $1.11\,$eV for the HS
configuration, respectively. Applying these parameters, one obtains an
adiabatic spin energy difference $\Delta E_{H-L}=-0.42\,$eV that erroneously
suggests a HS ground state for [Fe(CNH)_6]^2+ (the CASPT2/CC reference is
$+2.82\,$eV, see Figure 6).

Given that the majority of spurious contributions to the total energy arise
from the shell-averaged form of $U$, the implementation of orbital-resolved
DFT+$U\\!+\\!V$ is an intriguing prospect, as already hinted in the concluding
remarks of Ref. 53. This extension could exert precise control over electron
localization on molecular orbitals, which could be achieved through individual
Hubbard $V^{IJ}_{ij}$ parameters. Moreover, such highly tailored but still
fully first-principles corrections should allow to eliminate potential
conflicts arising from the simultaneous treatment of orbitals by both $U$ and
$V$ terms. With regard to the compounds discussed in the present work, the
orbital-resolved $U$ corrections of the localized $t_{g}$/$t_{2g}$ orbitals
could be augmented by $V$ terms specific to the interaction of $e_{g}$ with
neighboring ligand orbitals. Note that orbital-resolved inter-site $V$
parameters are readily obtained as a by-product of the calculation of $U$
parameters within the generalized LR-cDFT approach. However, practical
calculations with this fully-resolved approach would require the a priori
assignment of two-center occupation numbers to specific MOs, which involves
additional preparatory effort. Desirable will be the conception and
implementation of appropriate automation workflows in supporting tools such as
AiiDA142, 143. Alternatively, a simplified approach could involve using the
Hubbard $V$ correction in its current (averaged) form while retaining the
orbital-dependence of $U$. This would serve as a pragmatic ad-hoc solution to
the issues of $U$ at the expense of consistency.

Another avenue for extension pertains to Hund’s $J$ corrections, whose
combination with the orbital-resolved $U$ should be straightforward. Such
‘unlike-spin’ terms play a pivotal role in addressing the fractional spin
error that occurs in systems characterized by significant magnetic coupling
144, 62, 50, 145, 132, 84, 75.

##  6 Conclusions

In this study, we have introduced an orbital-resolved generalization of the
DFT$+U$ functional originally formulated by Dudarev et al.. Our implementation
is agnostic about the specific nature of the Hubbard projectors employed and
maintains full invariance against basis vector rotations. As a result,
calculations involving forces and stresses require no adjustments stemming
from this expansion. The true potential of the scheme emerges when the
required Hubbard parameters are derived from first principles. For this
purpose, we have employed an adapted version of the LR-cDFT approach, enabling
perturbative calculations with orbital resolution. This approach intrinsically
incorporates intra-shell screening, which often causes orbital-resolved $U$
parameters to be significantly smaller than their shell-averaged counterparts.
Provided a proper selection of the target manifold, the orbital resolution
therefore enables a more surgical use of DFT+$U$ that avoids overcorrections,
as comparative calculations of six Fe(II) molecular hexacomplexes as well as
of the charge-transfer insulators pyrite and pyrolusite underscore.

For instance, the orbital-resolved approach effectively addresses the bias in
adiabatic spin energies towards HS states observed in shell-averaged DFT$+U$
76. Particularly noteworthy is its success in accurately predicting spin
energies across a diverse spectrum of Fe(II) hexacomplexes, achieved through
selective corrections of the highly localized $t_{2g}$ orbitals. Explicit
fractional charge calculations on Fe[CNH]6^2+ suggest that these improvements
are not coincidental: the orbital-resolved approach reliably counteracts the
spurious global curvature of DFT with (semi-)local functionals with respect to
the fractional addition or removal of electrons, reducing the FCE by an order
of magnitude. In contrast, shell-averaged Hubbard $U$ corrections exhibit a
mixed performance: while they effectively diminish the FCE in the HS states,
their efficacy wanes when addressing the LS states. Here, instead of
rectifying the FCE, the convex error characteristic of (semi-)local DFT is
converted into a comparably pronounced concave error. The results also suggest
that a sole correction of the localized $t_{2g}$ orbitals performs slightly
better than a joint correction of $t_{2g}$ and the rather hybridized $e_{g}$
manifold, even when using orbital-resolved $U$ parameters. The superiority of
the orbital-resolved formulation over the shell-averaged approximation also
extends to the charge-transfer insulators FeS2 and β-MnO2. The correct non-
magnetic ground state of the former is only stabilized when orbital-resolved
Hubbard parameters are used. The magnitude of the large experimental band gap
of the latter can only be achieved by applying a pinpointed Hubbard correction
to the frontier O-$p_{z}$ orbitals, whereas standard Hubbard $U$ corrections
to the $d$ shell of Mn fail to open a significant band gap.

As these examples illustrate, the success of the orbital-resolved formulation
predominantly originates from the exclusion of hybridized orbitals from the
Hubbard manifold. From a theoretical perspective, the necessity to exclude
hybridized orbitals is rooted in the very definition of $U$ as an on-site
term. However, practically implementing this definition poses challenges,
particularly given that many DFT$+U$ investigations rely on atomic-like
orbitals as Hubbard projectors. This approach can lead to occupancy
eigenvalues far from 0 or 1, especially in compounds featuring covalent bonds.
These fractional values, however, are not related to the electron self-
interaction, as they arise from a one-center projector being applied to a two-
center phenomenon, namely a molecular orbital. In such scenarios, the orbital-
resolved approach provides an ad-hoc solution that allows to circumvent the
potential conflict of on-site corrections with the intricacies of covalently
bonded systems.

In a broader context, it is crucial to emphasize that the refined approach
presented in this work constitutes just one among various potent (Hubbard)
corrections available to DFT. As such, its effectiveness and ability to
achieve consistent enhancements relies on a thoughtful and technically well-
executed application. In practical terms, this is achieved by adopting a more
nuanced strategy for the determination of Hubbard manifolds, tailored to the
specific Hubbard projectors employed. Such a strategy should consider
localized ligand orbitals for Hubbard corrections in situations where these
act as frontier states. In addition, it might be desirable to address the FCE
in hybrid (i.e., molecular) orbitals. For this purpose, the use of Wannier
functions as Hubbard projectors or the adoption of an extended orbital-
resolved DFT$+U\\!+\\!V$ approach should be contemplated. A pivotal aspect of
our future investigations will focus on developing a protocol to automate the
choice of Hubbard manifolds based on measurable criteria. This step aims to
streamline and standardize the approach’s application, ensuring its systematic
and effective use across diverse systems and scenarios.

We thank Louis Ponet for fruitful discussions. EM was funded by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) – Project No.
286518848 – RTG 2247 $\mathrm{QM}^{3}$. IT and NM acknowledge support from the
Swiss National Science Foundation (SNSF) through its National Centre of
Competence in Research (NCCR) MARVEL (Grant No. 205602). NM acknowledges
funding from the U Bremen Excellence Chair programme. The authors also
acknowledge the computing time granted by the Resource Allocation Board and
provided on the supercomputers Lise and Emmy at NHR@ZIB and NHR@Göttingen as
part of the NHR infrastructure (project hbc00053) and the Swiss National
Supercomputing Centre (CSCS) under project No. s1073.

Shifts of VBM and CBM KS eigenvalues in pyrite due to Hubbard corrections.

## References

  * Hohenberg and Kohn 1964 Hohenberg, P.; Kohn, W. Inhomogeneous Electron Gas. _Phys. Rev._ 1964, _136_ , B864–B871
  * Kohn and Sham 1965 Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. _Phys. Rev._ 1965, _140_ , A1133–A1138
  * Hubbard 1963 Hubbard, J. Electron correlations in narrow energy bands. _Proc. R. Soc. A: Math. Phys._ 1963, _276_ , 238–257
  * Anisimov et al. 1991 Anisimov, V.; Zaanen, J.; Andersen, O. Band theory and Mott insulators: Hubbard $U$ instead of Stoner $I$. _Phys. Rev. B_ 1991, _44_ , 943
  * Anisimov et al. 1993 Anisimov, V. I.; Solovyev, I. V.; Korotin, M. A.; Czyżyk, M. T.; Sawatzky, G. A. Density-functional theory and NiO photoemission spectra. _Phys. Rev. B_ 1993, _48_ , 16929–16934
  * Liechtenstein et al. 1995 Liechtenstein, A. I.; Anisimov, V. I.; Zaanen, J. Density-functional theory and strong interactions: Orbital ordering in Mott-Hubbard insulators. _Phys. Rev. B_ 1995, _52_ , R5467–R5470
  * Cococcioni 2002 Cococcioni, M. A LDA+U Study of Selected Iron Compounds. Ph.D. thesis, International School for Advanced Studies (SISSA), Trieste, 2002
  * Cococcioni and de Gironcoli 2005 Cococcioni, M.; de Gironcoli, S. Linear response approach to the calculation of the effective interaction parameters in the LDA + U method. _Phys. Rev. B_ 2005, _71_ , 035105
  * Dudarev et al. 1998 Dudarev, S. L.; Botton, G. A.; Savrasov, S. Y.; Humphreys, C. J.; Sutton, A. P. Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+U study. _Phys. Rev. B_ 1998, _57_ , 1505–1509
  * Perdew et al. 1982 Perdew, J. P.; Parr, R. G.; Levy, M.; Balduz, J. L. Density-Functional Theory for Fractional Particle Number: Derivative Discontinuities of the Energy. _Phys. Rev. Lett._ 1982, _49_ , 1691–1694
  * Yang et al. 2000 Yang, W.; Zhang, Y.; Ayers, P. W. Degenerate Ground States and a Fractional Number of Electrons in Density and Reduced Density Matrix Functional Theory. _Phys. Rev. Lett._ 2000, _84_ , 5172–5175
  * Kulik et al. 2006 Kulik, H. J.; Cococcioni, M.; Scherlis, D. A.; Marzari, N. Density Functional Theory in Transition-Metal Chemistry: A Self-Consistent Hubbard U Approach. _Phys. Rev. Lett._ 2006, _97_ , 103001
  * Zhao and Kulik 2018 Zhao, Q.; Kulik, H. J. Where Does the Density Localize in the Solid State? Divergent Behavior for Hybrids and DFT+U. _J. Chem. Theory Comput._ 2018, _14_ , 670–683
  * Levy 1982 Levy, M. Electron densities in search of Hamiltonians. _Phys. Rev. A_ 1982, _26_ , 1200–1208
  * Zhao et al. 2016 Zhao, Q.; Ioannidis, E. I.; Kulik, H. J. Global and local curvature in density functional theory. _J. Chem. Phys._ 2016, _145_ , 054109
  * Dabo et al. 2010 Dabo, I.; Ferretti, A.; Poilvert, N.; Li, Y.; Marzari, N.; Cococcioni, M. Koopmans’ Condition for Density-Functional Theory. _Phys. Rev. B_ 2010, _82_ , 115121
  * Borghi et al. 2014 Borghi, G.; Ferretti, A.; Nguyen, N. L.; Dabo, I.; Marzari, N. Koopmans-Compliant Functionals and Their Performance against Reference Molecular Data. _Phys. Rev. B_ 2014, _90_ , 075135
  * Wang et al. 2006 Wang, L.; Maxisch, T.; Ceder, G. Oxidation energies of transition metal oxides within the GGA+U framework. _Phys. Rev. B_ 2006, _73_ , 195107
  * Le Bacq et al. 2004 Le Bacq, O.; Pasturel, A.; Bengone, O. Impact on electronic correlations on the structural stability, magnetism, and voltage of LiCoPO4 battery. _Phys. Rev. B_ 2004, _69_ , 245107
  * Aykol and Wolverton 2014 Aykol, M.; Wolverton, C. Local environment dependent GGA+$U$ method for accurate thermochemistry of transition metal compounds. _Phys. Rev. B_ 2014, _90_ , 115105
  * Isaacs and Marianetti 2017 Isaacs, E.; Marianetti, C. Compositional phase stability of strongly correlated electron materials within DFT+$U$. _Phys. Rev. B_ 2017, _95_ , 045141
  * Dederichs et al. 1984 Dederichs, P. H.; Blügel, S.; Zeller, R.; Akai, H. Ground States of Constrained Systems: Application to Cerium Impurities. _Phys. Rev. Lett._ 1984, _53_ , 2512–2515
  * McMahan et al. 1988 McMahan, A.; Martin, R.; Satpathy, S. Calculated effective Hamiltonian for La2CuO4 and solution in the impurity Anderson approximation. _Phys. Rev. B_ 1988, _38_ , 6650
  * Gunnarsson et al. 1989 Gunnarsson, O.; Andersen, O.; Jepsen, O.; Zaanen, J. Density-functional calculation of the parameters in the Anderson model: Application to Mn in CdTe. _Phys. Rev. B_ 1989, _39_ , 1708
  * Hybertsen et al. 1989 Hybertsen, M.; Schlüter, M.; Christensen, N. Calculation of Coulomb-interaction parameters for La2CuO4 using a constrained-density-functional approach. _Phys. Rev. B_ 1989, _39_ , 9028
  * Gunnarsson 1990 Gunnarsson, O. Calculation of parameters in model Hamiltonians. _Phys. Rev. B_ 1990, _41_ , 514
  * Pickett et al. 1998 Pickett, W. E.; Erwin, S. C.; Ethridge, E. C. Reformulation of the LDA + U method for a local-orbital basis. _Phys. Rev. B_ 1998, _58_ , 1201–1209
  * Solovyev and Imada 2005 Solovyev, I.; Imada, M. Screening of Coulomb interactions in transition metals. _Phys. Rev. B_ 2005, _71_ , 045103
  * Nakamura et al. 2006 Nakamura, K.; Arita, R.; Yoshimoto, Y.; Tsuneyuki, S. First-principles calculation of effective onsite Coulomb interactions of $3d$ transition metals: Constrained local density functional approach with maximally localized Wannier functions. _Phys. Rev. B_ 2006, _74_ , 235113
  * Shishkin and Sato 2016 Shishkin, M.; Sato, H. Self-consistent parametrization of DFT + U framework using linear response approach: Application to evaluation of redox potentials of battery cathodes. _Phys. Rev. B_ 2016, _93_ , 085135
  * Springer and Aryasetiawan 1998 Springer, M.; Aryasetiawan, F. Frequency-dependent screened interaction in Ni within the random-phase approximation. _Phys. Rev. B_ 1998, _57_ , 4364
  * Kotani 2000 Kotani, T. Ab initio random-phase-approximation calculation of the frequency-dependent effective interaction between 3d electrons: Ni, Fe, and MnO. _J. Phys. Condens. Matter_ 2000, _12_ , 2413
  * Aryasetiawan et al. 2004 Aryasetiawan, F.; Imada, M.; Georges, A.; Kotliar, G.; Biermann, S.; Lichtenstein, A. I. Frequency-dependent local interactions and low-energy effective models from electronic structure calculations. _Phys. Rev. B_ 2004, _70_ , 195104
  * Aryasetiawan et al. 2006 Aryasetiawan, F.; Karlsson, K.; Jepsen, O.; Scönberger, U. Calculations of Hubbard $U$ from first-principles. _Phys. Rev. B_ 2006, _74_ , 125106
  * Karlsson et al. 2010 Karlsson, K.; Aryasetiawan, F.; Jepsen, O. Method for calculating the electronic structure of correlated materials from a truly first-principles LDA + U scheme. _Phys. Rev. B_ 2010, _81_ , 245113
  * Mosey and Carter 2007 Mosey, N. J.; Carter, E. A. Ab initio evaluation of Coulomb and exchange parameters for DFT + U calculations. _Phys. Rev. B_ 2007, _76_ , 155123
  * Mosey et al. 2008 Mosey, N. J.; Liao, P.; Carter, E. A. Rotationally invariant ab initio evaluation of Coulomb and exchange parameters for DFT+U calculations. _J. Chem. Phys._ 2008, _129_ , 014103
  * Andriotis et al. 2010 Andriotis, A.; Sheetz, R.; Menon, M. LSDA+$U$ method: A calculation of the $U$ values at the Hartree-Fock level of approximation. _Phys. Rev. B_ 2010, _81_ , 245103
  * Agapito et al. 2015 Agapito, L.; Curtarolo, S.; Buongiorno Nardelli, M. Reformulation of DFT+$U$ as a Pseudohybrid Hubbard Density Functional for Accelerated Materials Discovery. _Phys. Rev. X_ 2015, _5_ , 011006
  * Tancogne-Dejean and Rubio 2020 Tancogne-Dejean, N.; Rubio, A. Parameter-free hybridlike functional based on an extended Hubbard model: DFT+U+V. _Phys. Rev. B_ 2020, _102_ , 155117
  * Lee and Son 2020 Lee, S.-H.; Son, Y.-W. First-principles approach with a pseudohybrid density functional for extended Hubbard interactions. _Phys. Rev. Res._ 2020, _2_ , 043410
  * Cococcioni and Marzari 2019 Cococcioni, M.; Marzari, N. Energetics and cathode voltages of Li M PO 4 olivines ( M = Fe , Mn) from extended Hubbard functionals. _Phys. Rev. Materials_ 2019, _3_ , 033801
  * Ricca et al. 2019 Ricca, C.; Timrov, I.; Cococcioni, M.; Marzari, N.; Aschauer, U. Self-consistent site-dependent DFT+$U$ study of stoichiometric and defective SrMnO3. _Phys. Rev. B_ 2019, _99_ , 094102
  * Sun et al. 2020 Sun, Y.; Cococcioni, M.; Wentzcovitch, R. _Phys. Rev. Mater._ 2020, _4_ , 063605
  * Kirchner-Hall et al. 2021 Kirchner-Hall, N. E.; Zhao, W.; Xiong, Y.; Timrov, I.; Dabo, I. Extensive Benchmarking of DFT+U Calculations for Predicting Band Gaps. _Appl. Sci._ 2021, _11_ , 2395
  * Xiong et al. 2021 Xiong, Y.; Campbell, Q.; Fanghanel, J.; Badding, C.; Wang, H.; Kirchner-Hall, N.; Theibault, M.; Timrov, I.; Mondschein, J.; Seth, K.; Katz, R.; Molina Villarino, A.; Pamuk, B.; Penrod, M.; Khan, M.; Rivera, T.; Smith, N.; Quintana, X.; Orbe, P.; Fennie, C.; Asem-Hiablie, S.; Young, J.; Deutsch, T.; Cococcioni, M.; Gopalan, V.; Abruña, H.; Schaak, R.; Dabo, I. Optimizing accuracy and efficacy in data-driven materials discovery for the solar production of hydrogen. _Energy Environ. Sci._ 2021, _14_ , 2335
  * Zhou et al. 2021 Zhou, J.-J.; Park, J.; Timrov, I.; Floris, A.; Cococcioni, M.; Marzari, N.; Bernardi, M. Ab Initio Electron-Phonon Interactions in Correlated Electron Systems. _Phys. Rev. Lett._ 2021, _127_ , 126404
  * Moore et al. 2022 Moore, G. C.; Horton, M. K.; Ganose, A. M.; Siron, M.; Linscott, E.; O’Regan, D. D.; Persson, K. A. High-throughput determination of Hubbard U and Hund J values for transition metal oxides via linear response formalism. 2022; http://arxiv.org/abs/2201.04213, arXiv:2201.04213
  * Haddadi et al. 2023 Haddadi, F.; Linscott, E.; Timrov, I.; Marzari, N.; Gibertini, M. On-Site and Inter-Site Hubbard Corrections in Magnetic Monolayers: The Case of FePS3 and CrI3. 2023; http://arxiv.org/abs/2306.06286, arXiv:2201.04213
  * Linscott et al. 2018 Linscott, E. B.; Cole, D. J.; Payne, M. C.; O’Regan, D. D. Role of spin in the calculation of Hubbard U and Hund’s J parameters from first principles. _Phys. Rev. B_ 2018, _98_ , 235157
  * Timrov et al. 2018 Timrov, I.; Marzari, N.; Cococcioni, M. Hubbard parameters from density-functional perturbation theory. _Phys. Rev. B_ 2018, _98_ , 085127
  * Timrov et al. 2021 Timrov, I.; Marzari, N.; Cococcioni, M. Self-consistent Hubbard parameters from density-functional perturbation theory in the ultrasoft and projector-augmented wave formulations. _Phys. Rev. B_ 2021, _103_ , 045141
  * Leiria Campo Jr and Cococcioni 2010 Leiria Campo Jr, V.; Cococcioni, M. Extended DFT + U \+ V method with on-site and inter-site electronic interactions. _J. Phys. Condens. Matter_ 2010, _22_ , 055602
  * Mahajan et al. 2021 Mahajan, R.; Timrov, I.; Marzari, N.; Kashyap, A. Importance of intersite Hubbard interactions in beta - MnO2 : A first-principles DFT+U+V study. _Phys. Rev. Mater._ 2021, _5_ , 104402
  * O’Regan et al. 2010 O’Regan, D. D.; Hine, N. D. M.; Payne, M. C.; Mostofi, A. A. Projector self-consistent DFT + U using nonorthogonal generalized Wannier functions. _Phys. Rev. B_ 2010, _82_ , 081102
  * Kulik and Marzari 2008 Kulik, H. J.; Marzari, N. A self-consistent Hubbard U density-functional theory approach to the addition-elimination reactions of hydrocarbons on bare FeO+. _J. Chem. Phys._ 2008, _129_ , 134314
  * Kraisler and Kronik 2013 Kraisler, E.; Kronik, L. Piecewise Linearity of Approximate Density Functionals Revisited: Implications for Frontier Orbital Energies. _Phys. Rev. Lett._ 2013, _110_ , 126403
  * Linscott et al. 2023 Linscott, E. B.; Colonna, N.; De Gennaro, R.; Nguyen, N. L.; Borghi, G.; Ferretti, A.; Dabo, I.; Marzari, N. Koopmans : An Open-Source Package for Accurately and Efficiently Predicting Spectral Properties with Koopmans Functionals. _J. Chem. Theory Comput._ 2023, _19_ , 7097–7111
  * Nguyen et al. 2018 Nguyen, N. L.; Colonna, N.; Ferretti, A.; Marzari, N. Koopmans-Compliant Spectral Functionals for Extended Systems. _Phys. Rev. X_ 2018, _8_ , 021051
  * Himmetoglu et al. 2014 Himmetoglu, B.; Floris, A.; de Gironcoli, S.; Cococcioni, M. Hubbard-corrected DFT energy functionals: The LDA+U description of correlated systems. _Int. J. Quantum Chem._ 2014, _114_ , 14–49
  * May and Kolpak 2020 May, K. J.; Kolpak, A. M. Improved description of perovskite oxide crystal structure and electronic properties using self-consistent Hubbard U corrections from ACBN0. _Phys. Rev. B_ 2020, _101_ , 165117
  * Orhan and O’Regan 2020 Orhan, O. K.; O’Regan, D. D. First-principles Hubbard U and Hund’s J corrected approximate density functional theory predicts an accurate fundamental gap in rutile and anatase TiO 2. _Phys. Rev. B_ 2020, _101_ , 245137
  * 63 Gelin, S.; Kirchner-Hall, N.; Katzbaer, R.; Theibault, M.; Xiong, Y.; Zhao, W.; Khan, M.; Andrewlavage, E.; Orbe, P.; Baksa, S.; Cococcioni, M.; Timrov, I.; Campbell, Q.; Abruna, H.; Schaak, R.; Dabo, I. Ternary oxides of $s-$ and $p-$block metals for photocatalytic solar-to-hydrogen conversion. _arXiv:2303.03332_
  * Timrov et al. 2020 Timrov, I.; Agrawal, P.; Zhang, X.; Erat, S.; Liu, R.; Braun, A.; Cococcioni, M.; Calandra, M.; Marzari, N.; Passerone, D. Electronic structure of Ni-substituted LaFeO3 from near edge x-ray absorption fine structure experiments and first-principles simulations. _Phys. Rev. Res._ 2020, _2_ , 033265
  * Ricca et al. 2020 Ricca, C.; Timrov, I.; Cococcioni, M.; Marzari, N.; Aschauer, U. Self-consistent DFT+$U$+$V$ study of oxygen vacancies in SrTiO3. _Phys. Rev. Res._ 2020, _2_ , 023313
  * Floris et al. 2020 Floris, A.; Timrov, I.; Himmetoglu, B.; Marzari, N.; de Gironcoli, S.; Cococcioni, M. Hubbard-corrected density functional perturbation theory with ultrasoft pseudopotentials. _Phys. Rev. B_ 2020, _101_ , 064305
  * Timrov et al. 2022 Timrov, I.; Aquilante, F.; Cococcioni, M.; Marzari, N. Accurate Electronic Properties and Intercalation Voltages of Olivine-Type Li-Ion Cathode Materials from Extended Hubbard Functionals. _PRX Energy_ 2022, _1_ , 033003
  * Kulik and Marzari 2011 Kulik, H. J.; Marzari, N. Transition-metal dioxides: A case for the intersite term in Hubbard-model functionals. _J. Chem. Phys._ 2011, _134_ , 094103
  * Mahajan et al. 2022 Mahajan, R.; Kashyap, A.; Timrov, I. Pivotal Role of Intersite Hubbard Interactions in Fe-Doped $\alpha$-MnO${}_{\textrm{2}}$. _J. Phys. Chem. C_ 2022, _126_ , 14353–14365
  * Timrov et al. 2023 Timrov, I.; Kotiuga, M.; Marzari, N. Unraveling the effects of inter-site Hubbard interactions in spinel Li-ion cathode materials. _Phys. Chem. Chem. Phys._ 2023, _25_ , 9061
  * Binci et al. 2023 Binci, L.; Kotiuga, M.; Timrov, I.; Marzari, N. Hybridization driving distortions and multiferroicity in rare-earth nickelates. _Phys. Rev. Research_ 2023, _5_ , 033146
  * 72 Bonfá, P.; Onuorah, I.; Lang, F.; Timrov, I.; Monacelli, L.; Wang, C.; Sun, X.; Petracic, O.; Pizzi, G.; Marzari, N.; Blundell, S.; De Renzi, R. Magnetostriction-driven muon localisation in an antiferromagnetic oxide. _arXiv:2305.12237_
  * 73 Gebreyesus, G.; Bastonero, L.; Kotiuga, M.; Marzari, N.; Timrov, I. Understanding the role of Hubbard corrections in the rhombohedral phase of BaTiO3. _arXiv:2309.04348_
  * Shishkin and Sato 2017 Shishkin, M.; Sato, H. Challenges in computational evaluation of redox and magnetic properties of Fe-based sulfate cathode materials of Li- and Na-ion batteries. _J. Phys. Condens. Matter_ 2017, _29_ , 215701
  * Shishkin and Sato 2021 Shishkin, M.; Sato, H. Evaluation of redox potentials of cathode materials of alkali-ion batteries using extended DFT+U+U up/down method: The role of interactions between the electrons with opposite spins. _J. Chem. Phys._ 2021, _154_ , 114709
  * Mariano et al. 2020 Mariano, L. A.; Vlaisavljevich, B.; Poloni, R. Biased Spin-State Energetics of Fe(II) Molecular Complexes within Density-Functional Theory and the Linear-Response Hubbard U Correction. _J. Chem. Theory Comput._ 2020, _16_ , 6755–6762
  * Mariano et al. 2021 Mariano, L. A.; Vlaisavljevich, B.; Poloni, R. Improved Spin-State Energy Differences of Fe(II) Molecular and Crystalline Complexes via the Hubbard U -Corrected Density. _J. Chem. Theory Comput._ 2021, _17_ , 2807–2816
  * O’Regan 2012 O’Regan, D. D. _Optimised Projections for the Ab Initio Simulation of Large and Strongly Correlated Systems_ ; Springer Berlin Heidelberg: Berlin, Heidelberg, 2012; pp 89–123
  * Solovyev et al. 1996 Solovyev, I.; Hamada, N.; Terakura, K. t23 versus all 3d localization in LaMO3 perovskites (M=Ti–Cu): First-principles study. _Phys. Rev. B_ 1996, _53_ , 7158–7170
  * Czyżyk and Sawatzky 1994 Czyżyk, M. T.; Sawatzky, G. A. Local-density functional and on-site correlations: The electronic structure of La2CuO4 and LaCuO3. _Phys. Rev. B_ 1994, _49_ , 14211–14228
  * Ylvisaker et al. 2009 Ylvisaker, E. R.; Pickett, W. E.; Koepernik, K. Anisotropy and magnetism in the LSDA + U method. _Phys. Rev. B_ 2009, _79_ , 035103
  * Kasinathan et al. 2007 Kasinathan, D.; Koepernik, K.; Pickett, W. E. Pressure-driven magnetic moment collapse in the ground state of MnO. _New Journal of Physics_ 2007, _9_ , 235–235
  * Tompsett et al. 2012 Tompsett, D. A.; Middlemiss, D. S.; Islam, M. S. Importance of Anisotropic Coulomb Interactions and Exchange to the Band Gap and Antiferromagnetism of beta-MnO2 from DFT+U. _Phys. Rev. B_ 2012, _86_ , 205126
  * Shishkin and Sato 2019 Shishkin, M.; Sato, H. DFT+ U in Dudarev’s formulation with corrected interactions between the electrons with opposite spins: The form of Hamiltonian, calculation of forces, and bandgap adjustments. _J. Chem. Phys._ 2019, _151_ , 024102
  * Georges et al. 1996 Georges, A.; Kotliar, G.; Krauth, W.; Rozenberg, M. J. Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions. _Rev. Mod. Phys._ 1996, _68_ , 13–125
  * Vaugier et al. 2012 Vaugier, L.; Jiang, H.; Biermann, S. Hubbard U and Hund exchange J in transition metal oxides: Screening versus localization trends from constrained random phase approximation. _Phys. Rev. B_ 2012, _86_ , 165105
  * Hampel et al. 2021 Hampel, A.; Lee-Hand, J.; Georges, A.; Dreyer, C. E. Correlation-induced octahedral rotations in SrMoO3. _Phys. Rev. B_ 2021, _104_ , 035102
  * Merkel and Ederer 2023 Merkel, M. E.; Ederer, C. Calculation of Screened Coulomb Interaction Parameters for the Charge-Disproportionated Insulator CaFeO3. 2023; arXiv:2310.16541
  * Prentice et al. 2020 Prentice, J. C. A.; Aarons, J.; Womack, J. C.; Allen, A. E. A.; Andrinopoulos, L.; Anton, L.; Bell, R. A.; Bhandari, A.; Bramley, G. A.; Charlton, R. J.; Clements, R. J.; Cole, D. J.; Constantinescu, G.; Corsetti, F.; Dubois, S. M.-M.; Duff, K. K. B.; Escartín, J. M.; Greco, A.; Hill, Q.; Lee, L. P.; Linscott, E.; O’Regan, D. D.; Phipps, M. J. S.; Ratcliff, L. E.; Serrano, A. R.; Tait, E. W.; Teobaldi, G.; Vitale, V.; Yeung, N.; Zuehlsdorff, T. J.; Dziedzic, J.; Haynes, P. D.; Hine, N. D. M.; Mostofi, A. A.; Payne, M. C.; Skylaris, C.-K. The ONETEP linear-scaling density functional theory program. _J. Chem. Phys._ 2020, _152_ , 174111
  * Kulik 2015 Kulik, H. J. Perspective: Treating electron over-delocalization with the DFT+U method. _J. Chem. Phys._ 2015, _142_ , 240901
  * Timrov et al. 2020 Timrov, I.; Aquilante, F.; Binci, L.; Cococcioni, M.; Marzari, N. Pulay forces in density-functional theory with extended Hubbard functionals: From nonorthogonalized to orthogonalized manifolds. _Phys. Rev. B_ 2020, _102_ , 235159
  * Kulik et al. 2016 Kulik, H. J.; Seelam, N.; Mar, B. D.; Martínez, T. J. Adapting DFT+ U for the Chemically Motivated Correction of Minimal Basis Set Incompleteness. _J. Phys. Chem. A_ 2016, _120_ , 5939–5949
  * Rohrbach et al. 2004 Rohrbach, A.; Hafner, J.; Kresse, G. Molecular adsorption on the surface of strongly correlated transition-metal oxides: A case study for CO/NiO(100). _Phys. Rev. B_ 2004, _69_ , 075413
  * Nawa et al. 2018 Nawa, K.; Akiyama, T.; Ito, T.; Nakamura, K.; Oguchi, T.; Weinert, M. Scaled effective on-site Coulomb interaction in the DFT+ U method for correlated materials. _Phys. Rev. B_ 2018, _97_ , 035117
  * Löwdin 1950 Löwdin, P. On the Non‐Orthogonality Problem Connected with the Use of Atomic Wave Functions in the Theory of Molecules and Crystals. _J. Chem. Phys._ 1950, _18_ , 365–375
  * Korotin et al. 2012 Korotin, D.; Kukolev, V.; Kozhevnikov, A.; Novoselov, D.; Anisimov, V. Electronic correlations and crystal structure distortions in BaBiO3. _J. Phys. Condens. Matter_ 2012, _24_ , 415603
  * Novoselov et al. 2015 Novoselov, D.; Korotin, D.; Anisimov, V. Hellmann-Feynman forces within the DFT+$U$ in Wannier functions basis. _J. Phys. Condens. Matter_ 2015, _27_ , 325602
  * Marzari and Vanderbilt 1997 Marzari, N.; Vanderbilt, D. Maximally localized generalized Wannier functions for composite energy bands. _Phys. Rev. B_ 1997, _56_ , 12847–12865
  * Marzari et al. 2012 Marzari, N.; Mostofi, A.; Yates, J.; Souza, I.; Vanderbilt, D. Maximally localized Wannier functions: Theory and applications. _Rev. Mod. Phys._ 2012, _84_ , 1419–1475
  * Qiao et al. 2023 Qiao, J.; Pizzi, G.; Marzari, N. Automated mixing of maximally localized Wannier functions into target manifolds. _npj Comp. Mater._ 2023, _9_
  * Tesch and Kowalski 2022 Tesch, R.; Kowalski, P. M. Hubbard U parameters for transition metals from first principles. _Phys. Rev. B_ 2022, _105_ , 195153
  * Sakuma and Aryasetiawan 2013 Sakuma, R.; Aryasetiawan, F. First-principles calculations of dynamical screened interactions for the transition metal oxides MO (M=Mn, Fe, Co, Ni). _Phys. Rev. B_ 2013, _87_ , 165118
  * Timrov et al. 2022 Timrov, I.; Marzari, N.; Cococcioni, M. HP – A code for the calculation of Hubbard parameters using density-functional perturbation theory. _Comput. Phys. Commun._ 2022, _279_ , 108455
  * Kulik and Marzari 2010 Kulik, H. J.; Marzari, N. Systematic study of first-row transition-metal diatomic molecules: A self-consistent DFT+U approach. _J. Chem. Phys._ 2010, _133_ , 114103
  * Yu and Carter 2014 Yu, K.; Carter, E. Communication: Comparing ab initio methods of obtaining effective U parameters for closed-shell materials. _J. Chem. Phys._ 2014, _140_ , 121105
  * Giannozzi et al. 2009 Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G.; Cococcioni, M.; Dabo, I.; Dal Corso, A.; De Gironcoli, S.; Fabris, S.; Fratesi, G.; Gebauer, R.; Gerstmann, U.; Gougoussis, C.; Kokalj, A.; Lazzeri, M.; Martin-Samos, L.; Marzari, N.; Mauri, F.; Mazzarello, R.; Paolini, S.; Pasquarello, A.; Paulatto, L.; Sbraccia, C.; Scandolo, S.; Sclauzero, G.; Seitsonen, A.; Smogunov, A.; Umari, P.; Wentzcovitch, R. Quantum ESPRESSO: A modular and open-source software project for quantum simulations of materials. _J. Phys. Condens. Matter_ 2009, _21_ , 395502
  * Giannozzi et al. 2017 Giannozzi, P.; Andreussi, O.; Brumme, T.; Bunau, O.; Buongiorno Nardelli, M.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Cococcioni, M.; Colonna, N.; Carnimeo, I.; Dal Corso, A.; de Gironcoli, S.; Delugas, P.; DiStasio Jr., R. A.; Ferretti, A.; Floris, A.; Fratesi, G.; Fugallo, G.; Gebauer, R.; Gerstmann, U.; Giustino, F.; Gorni, T.; Jia, J.; Kawamura, M.; Ko, H.-Y.; Kokalj, A.; Küçükbenli, E.; Lazzeri, M.; Marsili, M.; Marzari, N.; Mauri, F.; Nguyen, N. L.; Nguyen, H.-V.; Otero-de-la Rosa, A.; Paulatto, L.; Poncé, S.; Rocca, D.; Sabatini, R.; Santra, B.; Schlipf, M.; Seitsonen, A.; Smogunov, A.; Timrov, I.; Thonhauser, T.; Umari, P.; Vast, N.; Baroni, S. Advanced capabilities for materials modelling with Quantum ESPRESSO. _J. Phys. Condens. Matter_ 2017, _29_ , 465901
  * Giannozzi et al. 2020 Giannozzi, P.; Baseggio, O.; Bonfà, P.; Brunato, D.; Car, R.; Carnimeo, I.; Cavazzoni, C.; de Gironcoli, S.; Delugas, P.; Ferrari Ruffino, F.; Ferretti, A.; Marzari, N.; Timrov, I.; Urru, A.; Baroni, S. Quantum ESPRESSO toward the exascale. _J. Chem. Phys._ 2020, _152_ , 154105
  * Momma and Izumi 2011 Momma, K.; Izumi, F. VESTA for Three-Dimensional Visualization of Crystal, Volumetric and Morphology Data. _J. Appl. Crystallogr._ 2011, _44_ , 1272–1276
  * Fletcher 1987 Fletcher, R. _Practical Methods of Optimization_ , 2nd ed.; Wiley: Chichester, 1987
  * Perdew et al. 1996 Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. _Phys. Rev. Lett._ 1996, _77_ , 3865–3868
  * Kucukbenli et al. 2014 Kucukbenli, E.; Monni, M.; Adetunji, B. I.; Ge, X.; Adebayo, G. A.; Marzari, N.; de Gironcoli, S.; Corso, A. D. Projector augmented-wave and all-electron calculations across the periodic table: a comparison of structural and energetic properties. 2014; http://arxiv.org/abs/1404.3015, arXiv:1404.3015
  * Garrity et al. 2014 Garrity, K. F.; Bennett, J. W.; Rabe, K. M.; Vanderbilt, D. Pseudopotentials for high-throughput DFT calculations. _Comput. Mater. Sci._ 2014, _81_ , 446–452
  * Dal Corso 2014 Dal Corso, A. Pseudopotentials periodic table: From H to Pu. _Comput. Mater. Sci._ 2014, _95_ , 337–350
  * Lejaeghere et al. 2016 Lejaeghere, K.; Bihlmayer, G.; Björkman, T.; Blaha, P.; Blügel, S.; Blum, V.; Caliste, D.; Castelli, I. E.; Clark, S. J.; Dal Corso, A.; de Gironcoli, S.; Deutsch, T.; Dewhurst, J. K.; Di Marco, I.; Draxl, C.; Dułak, M.; Eriksson, O.; Flores-Livas, J. A.; Garrity, K. F.; Genovese, L.; Giannozzi, P.; Giantomassi, M.; Goedecker, S.; Gonze, X.; Grånäs, O.; Gross, E. K. U.; Gulans, A.; Gygi, F.; Hamann, D. R.; Hasnip, P. J.; Holzwarth, N. A. W.; Iuşan, D.; Jochym, D. B.; Jollet, F.; Jones, D.; Kresse, G.; Koepernik, K.; Küçükbenli, E.; Kvashnin, Y. O.; Locht, I. L. M.; Lubeck, S.; Marsman, M.; Marzari, N.; Nitzsche, U.; Nordström, L.; Ozaki, T.; Paulatto, L.; Pickard, C. J.; Poelmans, W.; Probert, M. I. J.; Refson, K.; Richter, M.; Rignanese, G.-M.; Saha, S.; Scheffler, M.; Schlipf, M.; Schwarz, K.; Sharma, S.; Tavazza, F.; Thunström, P.; Tkatchenko, A.; Torrent, M.; Vanderbilt, D.; van Setten, M. J.; Van Speybroeck, V.; Wills, J. M.; Yates, J. R.; Zhang, G.-X.; Cottenier, S. Reproducibility in density functional theory calculations of solids. _Science_ 2016, _351_
  * Prandini et al. 2018 Prandini, G.; Marrazzo, A.; Castelli, I. E.; Mounet, N.; Marzari, N. Precision and efficiency in solid-state pseudopotential calculations. _Npj Comput. Mater._ 2018, _4_ , 72
  * Moynihan et al. 2017 Moynihan, G.; Teobaldi, G.; O’Regan, D. D. A self-consistent ground-state formulation of the first-principles Hubbard U parameter validated on one-electron self-interaction error. 2017; http://arxiv.org/abs/1704.08076, arXiv:1704.08076
  * Zuñiga-Puelles et al. 2019 Zuñiga-Puelles, E.; Cardoso-Gil, R.; Bobnar, M.; Veremchuk, I.; Himcinschi, C.; Hennig, C.; Kortus, J.; Heide, G.; Gumeniuk, R. Structural stability and thermoelectric performance of high quality synthetic and natural pyrites (FeS2). _Dalton Trans._ 2019, _48_ , 10703–10713
  * Banjara et al. 2018 Banjara, D.; Malozovsky, Y.; Franklin, L.; Bagayoko, D. First-principles studies of electronic, transport and bulk properties of pyrite FeS2. _AIP Advances_ 2018, _8_ , 025212
  * Sun et al. 2011 Sun, R.; Chan, M. K. Y.; Ceder, G. First-principles electronic structure and relative stability of pyrite and marcasite: Implications for photovoltaic performance. _Phys. Rev. B_ 2011, _83_ , 235311
  * Eyert et al. 1998 Eyert, V.; Höck, K.-H.; Fiechter, S.; Tributsch, H. Electronic structure of FeS2: The crucial role of electron-lattice interaction. _Phys. Rev. B_ 1998, _57_ , 6350–6359
  * Armiento and Mattsson 2005 Armiento, R.; Mattsson, A. E. Functional designed to include surface effects in self-consistent density functional theory. _Phys. Rev. B_ 2005, _72_ , 085108
  * Burns 1977 Burns, G. _Introduction to group theory with applications_ ; Materials science and technology; Academic Press: New York, 1977
  * Schena et al. 2013 Schena, T.; Bihlmayer, G.; Blügel, S. First-principles studies of FeS2 using many-body perturbation theory in the G0W0 approximation. _Phys. Rev. B_ 2013, _88_ , 235203
  * Perdew et al. 1996 Perdew, J. P.; Ernzerhof, M.; Burke, K. Rationale for mixing exact exchange with density functional approximations. _J. Chem. Phys._ 1996, _105_ , 9982–9985
  * Heyd et al. 2003 Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Hybrid functionals based on a screened Coulomb potential. _J. Chem. Phys._ 2003, _118_ , 8207–8215
  * Cui et al. 2018 Cui, Z.-H.; Wang, Y.-C.; Zhang, M.-Y.; Xu, X.; Jiang, H. Doubly Screened Hybrid Functional: An Accurate First-Principles Approach for Both Narrow- and Wide-Gap Semiconductors. _J. Phys. Chem. Lett._ 2018, _9_ , 2338–2345
  * Zhao and Truhlar 2008 Zhao, Y.; Truhlar, D. G. The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals. _Theor. Chem. Acc._ 2008, _120_ , 215–241
  * Zhang and Jiang 2021 Zhang, M.-Y.; Jiang, H. Accurate Prediction of Band Structure of FeS2: A Hard Quest of Advanced First-Principles Approaches. _Front. Chem._ 2021, _9_ , 747972
  * Ennaoui and Tributsch 1984 Ennaoui, A.; Tributsch, H. Iron sulphide solar cells. _Sol. Cells_ 1984, _13_ , 197–200
  * Phung et al. 2018 Phung, Q. M.; Feldt, M.; Harvey, J. N.; Pierloot, K. Toward Highly Accurate Spin State Energetics in First-Row Transition Metal Complexes: A Combined CASPT2/CC Approach. _J. Chem. Theory Comput._ 2018, _14_ , 2446–2455
  * Bajaj et al. 2017 Bajaj, A.; Janet, J. P.; Kulik, H. J. Communication: Recovering the flat-plane condition in electronic structure theory at semi-local DFT cost. _J. Chem. Phys._ 2017, _147_ , 191101
  * Zaanen et al. 1985 Zaanen, J.; Sawatzky, G. A.; Allen, J. W. Band gaps and electronic structure of transition-metal compounds. _Phys. Rev. Lett._ 1985, _55_ , 418–421
  * Sato et al. 2000 Sato, H.; Enoki, T.; Isobe, M.; Ueda, Y. Transport properties and magnetism of a helically Hund-coupled conductor: beta-MnO2. _Phys. Rev. B_ 2000, _61_ , 3563–3569
  * Chevillot and Brenet 1959 Chevillot, J.-P.; Brenet, J. Influence d’ions étrangers sur la semi-conductibilité des bioxydes de manganèse. _C. R. Acad. Sci. Paris_ 1959, _278_ , 776–778
  * Druilhe 1966 Druilhe, R. Nouveaux Semi-Conducteurs Magnetiques, Cristal-Lisant Dans La Structure Du Rutile. _C. R. Hebd. Séances Acad. Sci., Ser. B_ 1966, 653
  * Kumar et al. 2017 Kumar, N.; Sen, A.; Rajendran, K.; Rameshbabu, R.; Ragupathi, J.; Therese, H. A.; Maiyalagan, T. Morphology and phase tuning of alpha- and beta-MnO ${}_{\textrm{2}}$ nanocacti evolved at varying modes of acid count for their well-coordinated energy storage and visible-light-driven photocatalytic behaviour. _RSC Adv._ 2017, _7_ , 25041–25053
  * Bigiani et al. 2020 Bigiani, L.; Gasparotto, A.; Maccato, C.; Sada, C.; Verbeeck, J.; Andreu, T.; Morante, J. R.; Barreca, D. Dual Improvement of beta‐MnO${}_{\textrm{2}}$ Oxygen Evolution Electrocatalysts via Combined Substrate Control and Surface Engineering. _ChemCatChem_ 2020, _12_ , 5984–5992
  * Franchini et al. 2007 Franchini, C.; Podloucky, R.; Paier, J.; Marsman, M.; Kresse, G. Ground-state properties of multivalent manganese oxides: Density functional and hybrid density functional calculations. _Phys. Rev. B_ 2007, _75_ , 195128
  * Lim et al. 2016 Lim, J. S.; Saldana-Greco, D.; Rappe, A. M. Improved pseudopotential transferability for magnetic and electronic properties of binary manganese oxides from DFT+U+J calculations. _Phys. Rev. B_ 2016, _94_ , 165151
  * Bolzan et al. 1993 Bolzan, A.; Fong, C.; Kennedy, B.; Howard, C. Powder Neutron Diffraction Study of Pyrolusite, beta-MnO2. _Aust. J. Chem._ 1993, _46_ , 939
  * Huber et al. 2020 Huber, S. P.; Zoupanos, S.; Uhrin, M.; Talirz, L.; Kahle, L.; Häuselmann, R.; Gresch, D.; Müller, T.; Yakutovich, A. V.; Andersen, C. W.; Ramirez, F. F.; Adorf, C. S.; Gargiulo, F.; Kumbhar, S.; Passaro, E.; Johnston, C.; Merkys, A.; Cepellotti, A.; Mounet, N.; Marzari, N.; Kozinsky, B.; Pizzi, G. AiiDA 1.0, a Scalable Computational Infrastructure for Automated Reproducible Workflows and Data Provenance. _Sci. Data_ 2020, _7_ , 300
  * Uhrin et al. 2021 Uhrin, M.; Huber, S. P.; Yu, J.; Marzari, N.; Pizzi, G. Workflows in AiiDA: Engineering a high-throughput, event-based engine for robust and modular computational workflows. _Comput. Mater. Sci._ 2021, _187_ , 110086
  * Himmetoglu et al. 2011 Himmetoglu, B.; Wentzcovitch, R. M.; Cococcioni, M. First-principles study of electronic and structural properties of CuO. _Phys. Rev. B_ 2011, _84_ , 115108
  * Mori-Sánchez et al. 2009 Mori-Sánchez, P.; Cohen, A. J.; Yang, W. Discontinuous Nature of the Exchange-Correlation Functional in Strongly Correlated Systems. _Phys. Rev. Lett._ 2009, _102_ , 066403