# Machine learning Hubbard parameters with equivariant neural networks

Martin Uhrin  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 (EPFL), CH-1015 Lausanne, Switzerland
Université Grenoble Alpes, 1130 Rue de la Piscine, BP 75, 38402 St Martin
D’Heres, France  Email<EMAIL_ADDRESS>Austin Zadoks  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 (EPFL), CH-1015 Lausanne, Switzerland  Luca Binci  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 (EPFL), CH-1015 Lausanne, Switzerland  Department of Materials
Science and Engineering, University of California at Berkeley, Berkeley,
California 94720, United States  Materials Sciences Division, Lawrence
Berkeley National Laboratory, Berkeley, CA, 94720, USA  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 (EPFL), CH-1015 Lausanne, Switzerland  Laboratory for Materials
Simulations (LMS), Paul Scherrer Institut (PSI), CH-5232 Villigen PSI,
Switzerland  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 (EPFL), CH-1015 Lausanne,
Switzerland  Laboratory for Materials Simulations (LMS), Paul Scherrer
Institut (PSI), CH-5232 Villigen PSI, Switzerland  Email<EMAIL_ADDRESS>
###### Abstract

Density-functional theory with extended Hubbard functionals (DFT+$U$+$V$)
provides a robust framework to accurately describe complex materials
containing transition-metal or rare-earth elements. It does so by mitigating
self-interaction errors inherent to semi-local functionals which are
particularly pronounced in systems with partially-filled $d$ and $f$
electronic states. However, achieving accuracy in this approach hinges upon
the accurate determination of the on-site $U$ and inter-site $V$ Hubbard
parameters. In practice, these are obtained either by semi-empirical tuning,
requiring prior knowledge, or, more correctly, by using predictive but
expensive first-principles calculations. Here, we present a machine learning
model based on equivariant neural networks which uses atomic occupation
matrices as descriptors, directly capturing the electronic structure, local
chemical environment, and oxidation states of the system at hand. We target
here the prediction of Hubbard parameters computed self-consistently with
iterative linear-response calculations, as implemented in density-functional
perturbation theory (DFPT), and structural relaxations. Remarkably, when
trained on data from 11 materials spanning various crystal structures and
compositions, our model achieves mean absolute relative errors of 3% and 5%
for Hubbard $U$ and $V$ parameters, respectively. By circumventing
computationally expensive DFT or DFPT self-consistent protocols, our model
significantly expedites the prediction of Hubbard parameters with negligible
computational overhead, while approaching the accuracy of DFPT. Moreover,
owing to its robust transferability, the model facilitates accelerated
materials discovery and design via high-throughput calculations, with
relevance for various technological applications.

## INTRODUCTION

A fundamental tool in investigating compounds involving transition-metal (TM)
and rare-earth (RE) compounds is density-functional theory (DFT), [1, 2] a
cornerstone for first-principles simulations in physics, chemistry, and
materials science. In practical applications, DFT necessitates approximations
to the exchange-correlation (xc) functional, with the local spin-density
approximation (LSDA) and spin-polarized generalized-gradient approximation
($\sigma$-GGA) being the most prevalent choices. However, these approximations
yield unsatisfactory outcomes for various properties of TM and RE compounds,
primarily due to significant self-interaction errors (SIEs) [3, 4, 5] that are
particularly pronounced for localized $d$ and $f$ electrons. To address these
challenges, more accurate approaches surpassing the limitations of “standard
DFT” have been devised. Noteworthy among these are Hubbard-corrected DFT (so-
called DFT+$U$ [6, 7, 8] and its extension DFT+$U$+$V$, [9, 10, 11] whose role
in addressing SIEs, rather than correlation errors, was first pointed out in
Ref. 4), meta-GGA functionals, [12, 13, 14] and hybrid functionals. [15, 16,
17] While these methods offer valuable insights, each comes with inherent
limitations and challenges, as discussed in, e.g., Ref. 18.

Initial structure, $U_{\text{in}},V_{\text{in}}$
DFT+$U_{\text{in}}$+$V_{\text{in}}$ geom. optimisation 1. DFPT calculation 2.
$|U_{\text{out}}-U_{\text{in}}|<\Delta$
$|V_{\text{out}}-V_{\text{in}}|<\Delta$ 3. $U_{\text{in}}=U_{\text{out}}$
$V_{\text{in}}=V_{\text{out}}$ Relaxed structure Self-consistent structure,
$U_{\text{SC}},V_{\text{SC}}$ Relaxed
structure$U_{\text{out}},V_{\text{out}}$NoYes

(a)

(b)

Figure 1: (a) Protocol for the self-consistent calculation of Hubbard
parameters using DFPT. [19] $U_{\mathrm{in}}$ and $V_{\mathrm{in}}$ represent
the input Hubbard parameters, while $U_{\mathrm{out}}$ and $V_{\mathrm{out}}$
denote the output parameters, with $\Delta$ representing the convergence
threshold. $U_{\mathrm{SC}}$ and $V_{\mathrm{SC}}$ are the final SC Hubbard
parameters. (b) Convergence of the Hubbard $U$ parameter for Mn-3d states in
LiMnPO4 using the self-consistent protocol.[18] The inset displays the crystal
structure of the material, where arrows indicate the spin direction, and Li
atoms are depicted in grey, O in red, Mn in violet, and P in yellow.

Hubbard-corrected DFT in particular stands out for its greater accuracy with
only a marginal increase in computational cost over standard DFT functionals.
[20] In the DFT+$U$+$V$ scheme, a corrective Hubbard energy $E_{U+V}$ is
introduced alongside the approximate DFT energy $E_{\mathrm{DFT}}$:

$E_{\mathrm{DFT}+U+V}=E_{\mathrm{DFT}}+E_{U+V}.$ (1)

In the simplified rotationally-invariant formulation, [8] the extended Hubbard
correction energy for a manifold with angular momentum $\ell$ takes the form:
[21]

$E_{U+V}=\frac{1}{2}\sum_{I\sigma}U^{I}\mathrm{Tr}\left[\mathbf{n}_{\ell}^{I}\left(\mathbf{1}-\mathbf{n}_{\ell}^{I\sigma}\right)\right]\\\
-\frac{1}{2}\sum_{I}\sum_{J(J\neq
I)}^{*}\sum_{\sigma}V^{IJ}\mathrm{Tr}\left[\mathbf{n}^{IJ\sigma}_{\ell}\mathbf{n}^{JI\sigma}_{\ell}\right]\,,$
(2)

where $I$ and $J$ are atomic site indices, while $m$ and $m^{\prime}$ denote
magnetic quantum numbers associated with a specific angular momentum. $U^{I}$
and $V^{IJ}$ are effective on-site and inter-site Hubbard parameters,
respectively. The asterisk in the sum signifies that, for each atom $I$, the
index $J$ encompasses all its neighbors up to a given distance. The
generalized occupation matrices $\mathbf{n}^{IJ\sigma}_{\ell}$ are derived
from the projection of the Kohn-Sham (KS) states onto localized atom-centered
orbitals $\phi^{I}_{m}(\mathbf{r})$ (Hubbard projector functions) of
neighboring atoms:

$\mathbf{n}_{\ell}^{IJ\sigma}\equiv
n^{IJ\sigma}_{mm^{\prime}}=\sum_{v,\mathbf{k}}f^{\sigma}_{v,\mathbf{k}}\braket{\psi^{\sigma}_{v,\mathbf{k}}}{\phi^{J}_{m^{\prime}}}\braket{\phi^{I}_{m}}{\psi^{\sigma}_{v,\mathbf{k}}}\,,$
(3)

where $v$ and $\sigma$ represent the band and spin labels of the KS
wavefunctions $\psi^{\sigma}_{v,\mathbf{k}}(\mathbf{r})$, respectively,
$\mathbf{k}$ denotes points in the first Brillouin zone (BZ), and
$f^{\sigma}_{v,\mathbf{k}}$ are the occupations of the KS states. The two
terms in eq. 2 — proportional to on-site $U^{I}$ and inter-site $V^{IJ}$ —
counteract one another. The on-site term promotes localization on atomic
sites, suppressing hybridization with neighbors, while the inter-site term
favors hybridized states with components on neighboring atoms. Consequently,
the values of $U^{I}$ and $V^{IJ}$ are critical to optimizing the extent of
the localization and hybridization within Hubbard-corrected DFT. However,
these parameters are not known a priori and must be determined in some way. In
passing we note that for the sake of simplicity, hereafter we drop the
superscripts $I$ and $J$ in the notations of Hubbard parameters and occupation
matrices, unless required for clarity.

Hubbard $U$ can be fit semi-empirically to reproduce a target property from
experimental data [22, 23, 24, 25, 26] or from other advanced first-principles
methods (e.g. $GW$ [27] or hybrids [15, 16, 17]); however, this approach has
many limitations. Fitted parameters are not guaranteed to, and often do not,
accurately predict properties other than that used in fitting, and calibrating
a single Hubbard $U$ to multiple properties is non-trivial, requiring advanced
algorithms like Bayesian optimization (BO). [28, 29] This method also
precludes materials discovery, where properties are by definition unknown, and
targeting results of advanced first-principles methods, like DFT with hybrid
functionals, inherits the limitations of those methods. [30, 31] Moreover, the
inter-site $V$ parameters, which are necessary to properly describe materials
with pronounced covalent interactions, are difficult to fit semi-empirically
due to the high-dimensional regression procedure they would require. An
attractive alternative is computing Hubbard parameters using first-principles
methods such as constrained DFT (cDFT), [32, 33, 34, 35, 36, 37, 38, 39, 40]
Hartree-Fock-based approaches, [41, 42, 43, 44, 10, 11] and the constrained
random-phase approximation (cRPA). [45, 46, 47, 48] These methods do not rely
on data from experiments or advanced simulation methods and, moreover, are
able to provide values not only for $U$ but also for $V$. However, first-
principles approaches are considerably more computationally expensive than DFT
+$U$(+$V$) ground-state calculations themselves, making their application
feasible but demanding for large systems or high-throughput studies.

The linear-response (LR) formulation of cDFT (LR-cDFT) [21] has witnessed
popularity due to its simplicity and accuracy; however, it demands
computationally expensive supercell calculations. A recent reformulation of
LR-cDFT in terms of DFPT [49, 19] significantly reduces the computational
burden for determining Hubbard parameters by replacing cumbersome supercell
calculations with faster unit-cell calculations and by leveraging symmetries
that further diminishes the computational cost by reducing the number of
perturbations in reciprocal space (see Sec. S1 in the supplemental information
(SI)). Its physical rationale relies on heuristically imposing piecewise
linearity of the energy of the system as a function of the population of the
Hubbard manifold. [21] Despite the numerous successful applications of DFPT in
computing Hubbard parameters, [50, 51, 18, 52, 53, 54, 55, 56, 57] this
approach introduces a significant overhead compared to DFT +$U$+$V$ ground-
state calculations. Moreover, it has been shown that jointly optimizing
Hubbard parameters and the crystal structure, rather than relying on the
equilibrium geometry obtained using (semi-)local functionals, can
significantly improve the accuracy of the final properties of interest. [58]
This is due in part to taking into account the geometry dependence of
$U\\{\mathbf{R}\\}$ and thus the contribution $dU/d\mathbf{R}$ to the
Hellmann-Feynman forces. [59] To do so, a self-consistent procedure combining
DFPT and structural optimizations (fig. 1(a)) can be used (see Sec. S2 in the
SI). [60, 19] Figure 1(b) shows a typical self-consistent determination of
Hubbard parameters for LiMnPO4 following this procedure, where each iteration
costs approximately one order of magnitude more computational time than a
Hubbard-corrected DFT electronic ground-state calculation. Consequently, the
adoption of acceleration techniques is highly sought after to expedite the
determination of Hubbard parameters while providing a level of accuracy very
close to that of DFPT and, prospectively, a simple estimate of
$dU/d\mathbf{R}$ and $dV/d\mathbf{R}$.

The application of machine learning (ML) algorithms to predict Hubbard
parameters is gaining popularity, as exemplified by a recent work, Ref. 61,
which employs a random forest regression model. Ref. 61 uses crystal structure
parameters, such as bond lengths and angles (among others), as descriptors for
an ML model trained on a database of Mn oxides. In this study, the training
data consist of $U$ values fitted using Bayesian optimization to replicate
band gaps and band structures obtained from HSE calculations — a concept
initially proposed in the earlier work of Ref. 28 within the context of
DFT+$U$ and subsequently extended to the DFT+$U$+$V$ framework. [62] The
resulting ML model predicts $U$ values such that subsequent DFT+$U$
calculations using these values yield band gaps that are $2-3$ times smaller
than the reference HSE values, while the band structures appear qualitatively
similar. Despite the significance of these efforts in applying ML to the
prediction of Hubbard parameters, numerous critical issues persist, such as
the accuracy of the reference training data (HSE often provides unreliable
band gaps for solids [30, 31]), the choice of the target properties for the
training data (fitting $U$ to reproduce band gaps is questionable), the
importance of self-consistency between Hubbard parameters and the crystal
structure, as well as the inclusion of inter-site $V$ which has been neglected
thus far, to name a few. Addressing these points is highly relevant and
necessitates further investigation, which serves as the motivation for this
work.

Here, we introduce a novel ML approach based on equivariant neural-networks
(ENNs) which aims to replace computationally demanding first-principles DFPT
calculations of Hubbard $U$ and $V$ parameters while providing a negligible
loss of accuracy for the vast majority of practical applications. Crucially,
the model employs: 1) atomic occupation matrices within the DFT+$U$+$V$
framework as descriptors of the geometry, electronic structure (e.g. oxidation
states (OSs)) and local chemical environments in materials, as well as 2)
DFPT-based Hubbard parameters, and 3) interatomic distances. This model is
general and can be applied to materials with ionic, covalent, and mixed ionic-
covalent interactions. It is trained on all intermediate Hubbard parameters
and occupation matrices obtained during the self-consistent cycle and directly
provides the final self-consistent values of $U$ and $V$, thus also bypassing
intermediate structural optimizations. The utilization of ENNs facilitates the
exploitation of the inherent O(3) group structure of the occupation matrices,
ensuring excellent model performance even with scarce training data.
Equivariant models have demonstrated state-of-the-art accuracy and
transferability in ML interaction potentials, [63, 64] while our work is the
first to incorporate electronic-structure degrees of freedom as explicit
features in solids.

## RESULTS

### Occupation matrices as the model inputs

The goal of our ML approach is to replace the self-consistent procedure
represented in fig. 1(a) and provide the final Hubbard parameters using as
input results from an initial DFT(+$U$+$V$) calculation. Conventional ML
methods for atomistic systems primarily use the ionic structures as inputs
(see e.g. Refs. 65, 66, 67, 68), side-stepping the need to explicitly
calculate or even consider the electronic structure. However, as a consequence
such models may not be particularly sensitive to changes in electronic
structure that only lead to subtle changes in local atomic geometry, as can
happen during a change of OS. To overcome this limitation, our models take as
input the on-site occupation matrices
$\mathbf{n}^{\sigma}_{\ell}\equiv\mathbf{n}_{\ell}^{II\sigma}$ [see eq. 3,
$I=J$], which describe the local electronic structure around an atom,
including the OS which is reflected in the occupation matrix eigenvalues. [69]
This is a particularly compelling choice given the significant variations in
Hubbard parameters seen for different OSs of TM elements. [18, 52]

As with any learning task, the ML model should respect the way in which the
target physical properties transform under global rotations, translations,
reflections and permutations of labels. In the present case, the output $U$
and $V$ values are left unchanged by these transformations, however the
entries of the occupation matrix _do_ change under rotation, making
$n^{\sigma}_{mm^{\prime}}$ unsuitable for use as inputs to non-symmetry aware
models. For this reason, we make use of an _equivariant_ learning model, [70,
71, 72, 73] which has the property that the learned function $f:X\to Y$ obeys
the following relation:

$D_{Y}[g]f(x)=f(D_{X}[g]x),\quad\forall g\in G,\,\,\forall x\in X,$ (4)

where the equivariance is with respect to a group $G$ (in our case SE(3)), and
$D_{Y}[g]$ and $D_{X}[g]$ are representations that act on the vector spaces
$Y$ and $X$, respectively. For example, applying a rotation to the inputs and
then applying $f$ must produce the same outputs as applying the rotation to
the outputs of $f$ given the original (unrotated) input. While our model
outputs scalars (which are symmetry invariant), by using an equivariant model
we ensure that the inputs and all intermediate (hidden) features transform
together under group actions. This property has been shown to give state-of-
the-art accuracy and transferability, particularly within the domain of ML
interaction potentials. [63, 74, 68]

### Equivariant descriptors

By virtue of being atom-centred, the on-site occupation matrix
$\mathbf{n}^{\sigma}_{\ell}$ is naturally invariant to global translations.
However, to build a rotationally invariant model, it is convenient to re-
express the occupation matrix in terms of the irreducible representations
(irreps) of the O(3) group. If we let $D^{(\ell,P)}$ be the irrep with degree
$\ell$ and parity $P$, we can express the occupation matrices as the tensor
products $D^{(\ell,P)}\otimes D^{(\ell,P)}$, where e.g. $\ell=1$ for $p$
orbitals, and $\ell=2$ for $d$ orbitals. These can then be decomposed into a
direct sum of irreps, e.g. for $\ell=1$ and $P=-1$ the occupation matrix is a
$3\times 3$, rank-2, tensor $D^{(1,-1)}\otimes D^{(1,-1)}$, that is symmetric,
as
$\mathbf{n}^{\sigma}_{\ell}=\left(\mathbf{n}^{\sigma}_{\ell}\right)^{\intercal}$.
In practice, this decomposition is achieved by applying a change of basis for
$\mathbf{n}^{\sigma}_{\ell}$ that transforms it into the irrep basis (see SI
S3 for details). At this stage, the remaining symmetry to be addressed is that
of label permutation. In the case of the atom labels for the inter-site term,
the model should be invariant to permutation of $I$ and $J$, however in our
training data the Hubbard $V$ correction is always applied to a
$d$-block/$p$-block atom pairs which have occupation matrices of different
dimension ($5\times 5$ for $d$ and $3\times 3$ for $p$), thereby making them
distinguishable.

Attribute | Irreducible representation
---|---
$U^{I}$, $V^{IJ}$, $r_{IJ}$ | $D^{(0,1)}$
Atomic species | $\bigoplus^{N_{s}}D^{(0,1)}$
$\mathbf{n}^{\sigma}_{\mathrm{p}}$ | $D^{(0,1)}\oplus D^{(2,1)}$
$\mathbf{n}^{\sigma}_{\mathrm{d}}$ | $D^{(0,1)}\oplus D^{(2,1)}\oplus D^{(4,1)}$
Table 1: Irreducible representations of various attributes (inputs and
outputs) of the ML model. $U^{I}$ and $V^{IJ}$ are the on-site and inter-site
Hubbard parameters, respectively, $r_{IJ}$ is the interatomic distance between
species $I$ and $J$, $N_{s}$ is the number of atomic species (one-hot
vectors), $\mathbf{n}^{\sigma}_{\mathrm{p}}$ and
$\mathbf{n}^{\sigma}_{\mathrm{d}}$ are the occupation matrices for the p and d
orbitals, respectively, in the special representation for the ML model [see
eqs. 5 and 6]. The irreducible representations of the occupation matrices can
be readily obtained using the e3nn [73, 75] library.

$N$ layers$\bullet$ $U_{\text{in}}$$\bullet$ Species$\bullet$
$\mathbf{x}_{\text{d}}^{1}$$\bullet$ $\mathbf{x}_{\text{d}}^{2}$ $\bigotimes$
Gated non-linearity $\bigotimes$$U_{\text{out}}$= Input/output tensor= Hidden
state$\otimes$= Tensor product (with learnable weights) (a)

$N$ layers$N$ layers$\bullet$ Species$\bullet$
$\mathbf{x}_{\text{d}}^{1}$$\bullet$ $\mathbf{x}_{\text{d}}^{2}$ $\bullet$
Species$\bullet$ $\mathbf{x}_{\text{p}}^{1}$$\bullet$
$\mathbf{x}_{\text{p}}^{2}$ $\bullet$ $V_{\text{in}}$$\bullet$ $r_{IJ}$
$\bigotimes$$\bigotimes$ Gated non-linearity $\bigotimes$$V_{\text{out}}$ (b)

Figure 2: Schematic illustration of the equivariant neural-network ML model
for predicting on-site Hubbard $U$ (a) and inter-site Hubbard $V$ (b)
parameters. $U_{\mathrm{in}}$ and $V_{\mathrm{in}}$ are the input Hubbard
parameters, while $U_{\mathrm{out}}$ and $V_{\mathrm{out}}$ are the outputs.
The atomic species enter as one-hot tensors, and $r_{IJ}$ is the interatomic
distance between sites $I$ and $J$. $\mathbf{x}^{1}_{\mathrm{d}}$ and
$\mathbf{x}^{2}_{\mathrm{d}}$ are the occupation matrices for the $d$
orbitals, while $\mathbf{x}^{1}_{\mathrm{p}}$ and
$\mathbf{x}^{2}_{\mathrm{p}}$ for the $p$ orbitals, all in the special
representation for the ML model [see eqs. 5 and 6]. Tensors are represented as
squares: open for inputs and outputs, and filled for intermediate features.

There remains a permutational invariance to be imposed on the spin labels
($\sigma=\uparrow$ or $\downarrow$). We achieve this using permutationally
invariant polynomials:

$\displaystyle\mathbf{x}^{1}_{\ell}$
$\displaystyle=\mathbf{n}_{\ell}^{\uparrow}+\mathbf{n}_{\ell}^{\downarrow}\,,$
(5) $\displaystyle\mathbf{x}^{2}_{\ell}$
$\displaystyle=\mathbf{n}_{\ell}^{\uparrow}\otimes\mathbf{n}_{\ell}^{\downarrow}\,.$
(6)

In the following, for the sake of convenience, we label these tensors in their
irrep form as $\mathbf{x}^{i}_{\mathrm{p}}$ for $p$ orbitals and
$\mathbf{x}^{i}_{\mathrm{d}}$ for $d$ orbitals, where $i$ can be 1 or 2. We
note in passing that in the non-spin-polarized case, one can set
$\mathbf{n}_{\ell}^{\uparrow}=\mathbf{n}_{\ell}^{\downarrow}=\mathbf{n}_{\ell}/2$,
where $\mathbf{n}_{\ell}$ is divided by two to account for spin degeneracy,
and then proceed as above. In addition, atomic species (i.e. atomic types) are
encoded as one-hot vectors in our ML model. For example, for a model that
supports $N_{s}$ distinct atomic species, we choose an arbitrary labelling of
a $N_{s}$ long vector, $[\text{Mn},\text{Fe},\text{Ni},\text{O},\ldots]$, and
place a 1 at the entry corresponding to the site element and zeros everywhere
else (Mn = $[1,0,0,0,0,\ldots]$, Fe = $[0,1,0,0,0,\ldots]$, etc). A summary of
the irreps of all the model inputs is presented in Table 1. In the case of the
ML model for Hubbard $V$, the interatomic distances $r_{JI}$ are also included
as an additional input, giving some basic information of the ionic structure.

### Equivariant neural-network model

To construct our ML model, we use the e3nn library [73, 75] and PyTorch; [76]
our codebase is open-source and freely available (see code availability
below). As shown in fig. 2, we define two separate models, one for predicting
the on-site $U$ values and the other for the inter-site $V$ ones. By using
separate models, we provide flexibility for the model to be used in
calculations where only Hubbard $U$ is applied (i.e. DFT+$U$). In each case,
the model starts with one or more nodes that carry attributes, expressed as a
direct sum of irreps, which represent the inputs to the learned function. The
inputs pass through a series of repeated layers made up of a tensor product
(analogous to all-to-all connected layers in a traditional neural network) and
a gated non-linearity (an equivariant version of a traditional activation
function [75]) followed by a final tensor product before readout. All of the
tensor products have learnable weights, meaning that every pair of input
irreps that contribute to one output irrep has a learnable scalar parameter
that is optimised during training. We use the AdamW [77] optimiser to minimise
the loss, which is the mean squared error between predicted and training
Hubbard parameters.

### Training and validation datasets

To train and validate our ENN model, we curated a dataset comprising materials
with diverse crystal structures. This dataset is constructed based on an
investigation of various Li-ion cathode materials, covering olivine-type, [18]
spinel-type, [52] and layered-type [78] structures at different Li
concentrations, employing the self-consistent DFT+$U$ and DFT+$U$+$V$
protocols outlined above. Additionally, we include materials such as tunnel-
and rutile-type MnO2 [79, 51] and perovskite-type rare-earth nickelates [53]
in our dataset. The crystal and electronic structure details, magnetic
ordering, and various other properties of these materials are extensively
documented in their respective publications. These are the compounds
containing TM elements Ni, Mn, and Fe, and a list of all these materials is
provided in table 2. It is crucial to note that, while our current dataset is
relatively small in terms of material count, ongoing efforts are directed
towards a high-throughput exploration involving hundreds of materials, which
will serve as an expanded training and validation set for our ENN model in the
future. As a proof of concept, this initial dataset demonstrates the efficacy
of our ML model, showcasing its ability to predict Hubbard parameters even
with a limited dataset. Despite the size of the dataset, the model exhibits
accurate predictions, setting the stage for further refinement and expansion
with the upcoming comprehensive investigation.

Crystal structure type | Chemical composition | #
---|---|---
Olivine | LixFePO4 | 5
LixMnPO4 | 5
LixFe0.5Mn0.5PO4 | 5
Spinel | LixMn2O4 | 2
LixMn1.5Ni0.5O4 | 2
Layered | LixNiO2 | 2
LixMnO2 | 2
Tunnel |  $\alpha$-MnO2 | 1
Rutile |  $\beta$-MnO2 | 1
Perovskite | YNiO3 | 1
PrNiO3 | 1
Total |  | 27
Table 2: A list of the materials that are used to train and validate the ML
model. The last column shows the number of materials. For the olivine-type
materials $x=0,0.25,0.50,0.75,0$, while for the spinel- and layered-type
materials $x=0,1$.

For each material in our dataset, SC Hubbard $U$ and $V$ parameters are
computed using the self-consistent protocol illustrated in fig. 1(a). We
include all $V$ parameters for Hubbard-active atom pairs whose DFPT values are
greater than $0.3$ eV, while in Ref. 62 only $V$ for the nearest neighbor
couples were included in the training of their ML model. Throughout the
iterative process, all intermediate converged occupation matrices and Hubbard
parameters are systematically saved. This dataset is then used for training
our ENN model. The training data, as summarized in table 3, is quite
substantial when considering the multitude of these first-principles
calculations. The dataset for training $V$ comprises approximately seven times
more data points than that for training $U$ due to the greater abundance of
atom pairs involved in the computation of $V$, in contrast to the local and
consequently sparser nature of $U$. For each material and chemical element,
this data is randomly split, holding back 20% from each material for
validation.

Chemical element | Number of data points
---|---
| $U$ | $V$
Ni | 396 (124) | 67,802 (63,332)
Mn | 856 (284) | 162,272 (153,232)
Fe | 138 (120) | 22,511 (21,483)
Total | 1,390 (528) | 252,585 (238,047)
Table 3: Chemical elements and the total number of data points for both on-
site $U$ and inter-site $V$ Hubbard parameters determined using the SC
protocol illustrated in fig. 1(a). These data points are aggregated across all
the materials listed in table 2 that contain the respective chemical element.
Numbers in brackets indicate the unique data points after de-duplication (see
text).

In the first-principles calculations, certain input Hubbard parameters
($U_{\mathrm{in}}$ and $V_{\mathrm{in}}$) yield a specific ground state and
occupation matrices, and the response of the system due to a perturbation
through DFPT provides corresponding output Hubbard parameters
($U_{\mathrm{out}}$ and $V_{\mathrm{out}}$). This one-to-one correspondence
allows our ML model to discern the mapping between input and output Hubbard
parameters and the associated intermediate occupation matrices. In addition,
it is worth noting that many of the structures encompass symmetry-equivalent
sites, and correlations between results from successive self-consistent steps
are prevalent. To address this issue, we implement a de-duplication procedure
aimed at mitigating potential biases in our training data. This involves
calculating a symmetry-invariant distance between all pairs of ENN inputs for
each atomic species in each material, which is then used to cluster duplicates
based on a specified distance threshold (see Sec. S4 in the SI). In random
train/validate splits, we sample from these clusters rather than individual
data points, ensuring that only one example from each cluster is included in
the validation set.

### Accuracy of the ML model in predicting the SC Hubbard parameters

First, we aim to learn the final SC Hubbard parameters for the TM elements
across all the materials in our dataset. This is achieved using all attributes
listed in table 1 as inputs for our ML model. Figure 3 shows the parity plots
obtained for the on-site $U$ and inter-site $V$ Hubbard parameters. The mean
absolute relative error (MARE) over all species consistently remains below 3%
and 5% for the $U$ and $V$ parameters, respectively (per-species distributions
are reported in Sec. S5 in the SI). The result for Hubbard $U$ is particularly
promising, considering the relatively small amount of available training data.
The higher MARE obtained for $V$ is likely due to the significantly increased
number of degrees of freedom in its model, owing to the involvement of pairs
of atoms. It is worth noting that in order to attain this accuracy in
predicting the $U$ and $V$ values, utilizing a relatively small batch size of
$8-16$ was beneficial in preventing overfitting.

(a)

(b)

Figure 3: Parity plots showing the prediction accuracy on an unseen validation
dataset, where the energies in the legend are the RMSE categorized by
element(s) and the overall RMSE across all elements. All attributes listed in
table 1 are used as inputs for the ML model. (a) Hubbard $U$ for 3d states of
TM elements, (b) Hubbard $V$ between the TM-3d and O-2p states.

Second, we examine the extent to which our ML model depends upon the input
Hubbard parameters for its prediction accuracy. To investigate this, we repeat
the same numerical experiment as above, eliminating the input Hubbard
parameters ($U_{\text{in}}$ and $V_{\text{in}}$) from our model. As detailed
in Sec. S6 in the SI, the parity plots closely resemble those in fig. 3. The
reduction in accuracy of our ML model resulting from this simplification is
relatively small, namely the overall root-mean-square error (RMSE) is
increases by 29 and 18 meV for the on-site $U$ and inter-site $V$,
respectively. This finding suggests that the input Hubbard parameters may not
be the most critical input attributes for our ML model, and that the input
occupation matrices contain sufficient information to make low-error
predictions.

$x$ | Property | LixMnPO4 | LixFePO4
---|---|---|---
$0-1$ |  $\Phi$ (DFPT) | 4.205 | 3.544
$\Phi$ (ML) | 4.194 | 3.544
$\Delta\Phi$ |  $-0.26$% | 0.00%
0 |  $m$ (DFPT) | 3.9738 | 4.1828
$m$ (ML) | 3.9720 | 4.1832
$\Delta m$ |  $-0.05$% |  $0.01$%
1 |  $m$ (DFPT) | 4.7482 | 3.7388
$m$ (ML) | 4.7486 | 3.7391
$\Delta m$ |  $0.01$% |  $0.01$%
Table 4: Comparison of the open-circuit voltages $\Phi$ (in V) and magnetic
moments for TM elements (in $\mu_{\mathrm{B}}$) for LixMnPO4 and LixFePO4
computed within DFT+$U$+$V$, with $U$ and $V$ obtained from first principles
using DFPT and predicted using the ML model. The voltages are computed using
the total energy differences for the Li concentrations $x=0$ and $x=1$, while
the magnetic moments are computed as the trace of the difference between the
spin-up and spin-down occupation matrices. [18] $\Delta\Phi$ and $\Delta m$
are the relative differences between the voltages and magnetic moments based
on DFPT and ML, respectively.

Although the model’s relative error is low, it is important to evaluate its
impact on downstream properties, such as e.g. voltages in Li-ion battery
cathode materials. Table 4 shows a comparison of calculated open-circuit
voltage (OCV) [18] within the DFT+$U$+$V$ framework, using first-principles SC
Hubbard parameters from DFPT and those predicted by our ML model. We find that
the differences between the computed and predicted $U$ and $V$ values are less
than $0.12$ eV for LixMnPO4 and LixFePO4 ($x=0$ and $x=1$). Generally, such a
small variation in the values of Hubbard parameters has a negligible impact on
the vast majority of various physical and chemical properties of materials.
The differences observed in the OCV in Table 4 are indeed negligible,
indicating the accuracy and reliability of the ML-predicted Hubbard
parameters.

### ML model’s performance using a reduced number of iterations in the SC
protocol

In the preceding section, we assessed the ML model’s performance in predicting
final self-consistent Hubbard parameters using training data that consists of
a subset of DFPT calculations from all iterations of the SC protocol (see fig.
1(a)). However, this approach relies on conducting numerous computationally
intensive DFPT calculations to generate these training data (often $2-5$ but
occasionally up to $10$ per material to reach self-consistency). In this
experiment, we assess the performance of the model by training on only the
first $N_{\mathrm{iter}}-1$ DFPT result, and asking it to predict the outcome
of the next DFPT calculation. This gives an indication of whether the model
can be useful in a scenario where a limited number of self-consistent steps
are carried out, such as during the early stages of a high-throughput
screening study.

Figure 4 shows the average RMSE over all elements as the number of SC
iterations used for training $N_{\mathrm{iter}}$ increases. The reference data
corresponds to RMSE evaluated using DFPT calculations, i.e. excluding the ML
component. This RMSE is computed by determining the difference between the
Hubbard parameters at the current $N_{\mathrm{iter}}$th and the previous
$(N_{\mathrm{iter}}-1)$th iteration. It is evident that the reference RMSE
decreases non-monotonically, which was observed in previous studies, [19] and
eventually diminishes at large $N_{\mathrm{iter}}$. It is noteworthy that the
optimal $N_{\mathrm{iter}}$ for achieving self-consistency in the SC protocol
varies for different materials listed in table 2. Therefore, for larger
$N_{\mathrm{iter}}$ values in fig. 4, fewer data points are available compared
to smaller $N_{\mathrm{iter}}$ values.

Figure 4: RMSE as a function of the number of iterations $N_{\mathrm{iter}}$
in the SC protocol (see fig. 1(a)). The reference data represents the RMSE
obtained from successive DFPT calculations (without ML involvement), whereas
the ML data denotes the RMSE resulting from ML predictions for the
$N_{\mathrm{iter}}$th iteration based on training the model on all preceding
$N_{\mathrm{iter}}-1$ iterations, with validation conducted using the DFPT
data for the $N_{\mathrm{iter}}$th iteration as a reference.

The ML-based RMSE in fig. 4 at the $N_{\mathrm{iter}}$th iteration is computed
by training the model on data from all previous $N_{\mathrm{iter}}-1$
iterations and predicting the Hubbard parameters for the $N_{\mathrm{iter}}$th
iteration. Section S7 in the SI contains parity plots for various values of
$N_{\mathrm{iter}}$, ranging from 2 to 7, alongside the reference results
broken down by element as well as the results for the training set. It is
evident from fig. 4 that the ML-based RMSE steadily decreases, and for all
$N_{\mathrm{iter}}$ values except $N_{\mathrm{iter}}=6$, it is smaller than
the RMSE of the reference dataset. This result is shows that our ML model can
improve upon a DFPT-based SC result that was terminated early, particularly
during the first few iterations, thereby facilitating faster convergence of
Hubbard parameters. For instance, the reference value for
$N_{\mathrm{iter}}=2$ indicates that on average, a single-shot DFPT
calculation yields Hubbard parameters with an RMSE of $\sim 1$ eV, while the
ML model trained on single-shot data can reduce this error to under $0.8$ eV.
Incorporating additional $N_{\mathrm{iter}}$ iterations in the ML training
data leads to steady improvements, although the relative performance compared
to the reference gradually diminishes until reaching the floor of the model’s
accuracy. This floor is partly determined by the convergence of the underlying
DFPT-based Hubbard parameters in the SC protocol, conducted with a tolerance
of approximately $\Delta=0.01-0.1$ eV (see fig. 1(a)). Therefore, the most
significant computational savings in predicting the SC Hubbard parameters can
be achieved using our ML model, which acts as a surrogate, after only a few
iterations in the DFPT-based SC protocol. This substantially reduces the
overall computational cost and renders it a feasible approach, particularly
for high-throughput screening scenarios where a priori training data is
unavailable for all materials.

### Transferability of the ML model

To test the transferability of the model, we isolate training data from the
olivines (Lix$M$PO4 with $M$ = Fe, Mn, or
Fe${}_{\text{0.5}}$Mn${}_{\text{0.5}}$) for which we have calculations at Li
concentrations $x=0,0.25,0.5,0.75,1$. Within this class of materials, the OS
of TM elements Mn and Fe changes from $+2$ to $+3$ upon delithiation. [80] As
mentioned earlier, this change of OS is directly reflected in changes to the
occupation matrices. [18] Consequently, the Hubbard parameters for TM elements
at different OS also exhibit variation. For instance, the SC Hubbard $U$
parameter for Mn changes from 4.56 to 6.26 eV when transitioning from $+2$ to
$+3$, while for Fe this change is from 5.29 to 5.43 eV. [18] In addition, in
Sec. S8 of the SI we present the distribution of Hubbard $U$ parameters for Fe
and Mn ions in the olivines extracted from the DFPT-based SC protocol for
various Li concentrations. Our ML model effectively captures these changes in
the Hubbard parameters and accurately predicts their values based on the
occupation matrices for each TM ion at different Li concentrations.

The outcomes of our numerical experiments are shown in fig. 5. As a reference,
we compute the RMSE individually for each material by training the ML model on
80% of the data from all five concentrations and validating it on the
remaining 20% of the data after de-duplication. This yields RMSE values of
126, 182, and 420 meV for the Fe, Mn, and mixed Fe-Mn olivines, respectively
(see the horizontal lines in fig. 5). Subsequently, we investigate the RMSE
values for scenarios in which the ML model is trained on fewer concentrations,
denoted as $N_{\mathrm{c}}$, to assess the sensitivity of the ML model to the
amount of the training data and its transferability for predicting Hubbard
parameters at other concentrations. To accomplish this, we train the ML model
on $N_{\mathrm{c}}$ concentrations and validate it on the remaining
$5-N_{\mathrm{c}}$ concentrations, and then we average out over all possible
permutations of the concentrations between the training and validation
datasets. The average RMSE values are represented by dots, while the error
bars in fig. 5 indicate the maximum and minimum RMSE values resulting from the
various permutations. It can be seen that as $N_{\mathrm{c}}$ increases, the
average RMSE for each material decreases, reaching values of approximately
106, 301, and 652 meV at $N_{\mathrm{c}}=4$ for the Fe, Mn, and mixed Fe-Mn
olivines, respectively.

This experiment demonstrates that, depending on the desired accuracy, it is
not necessary to include training data from the exact compositions that the
model is being evaluated on, and even using examples at one or two
concentrations may be sufficient to extrapolate to other parts of the
composition range. Given the computational cost of SC calculations, this can
lead to a significant speedup, particularly in cases requiring larger
supercells and various configurations of Li ion distributions within the
structure. Overall, these findings demonstrate the good transferability of our
ML model, indicating its potential utility and reliability for predicting
Hubbard parameters in materials not included in the model’s training dataset.

(a) LixFePO4 (Ref. $=0.126$ eV)

(b) LixMnPO4 (Ref. $=0.182$ eV)

(c) LixFe0.5Mn0.5PO4 (Ref. $=0.420$ eV)

Figure 5: RMSE as a function of $N_{\mathrm{c}}$ (the number of Li
concentrations $x$) for the three olivines: (a) LixFePO4, (b) LixMnPO4, and
(c) LixFe0.5Mn0.5PO4. The ML model is trained and validated using
$N_{\mathrm{c}}$ and the remaining $5-N_{\mathrm{c}}$ concentrations,
respectively. Error bars indicate the range of RMSE values obtained by
considering various permutations of concentrations in the training and
validation datasets, the dots represent the average. Reference data is the
RMSE computed by training the ML model on 80% of all five concentrations, with
validation performed on the remaining 20% after de-duplication. The line
represents the mean results over three runs with different random
initializations, the confidence interval show the standard deviation.

## DISCUSSION

We have introduced a novel equivariant ML model designed for predicting the
self-consistent on-site $U$ and inter-site $V$ Hubbard parameters, thereby
circumventing the computationally intensive DFPT-based protocols. The model
incorporates three input descriptors: Hubbard parameters, inter-atomic
distances, and, notably, atomic occupation matrices. The latter play a pivotal
role in encoding essential information about the electronic structure and
local chemical environment within materials. Such an ML model holds
significant promise, particularly for high-throughput investigations and
large-scale systems, scenarios where DFPT-based approaches become too
computationally expensive. Furthermore, the model demonstrates good
transferability, rendering it reliable for predicting Hubbard parameters in
materials not included in its training dataset, or as a very accurate first
guess for further DFPT refinements.

The usage of our ML model is straightforward and entails two DFT-based
calculations for a given material. Initially, a ground-state calculation
employing DFT+$U$+$V$ with initial guesses for $U$ and $V$ (which can be set
to zero) is conducted to determine the atomic occupation matrices required as
input for the model. Subsequently, a final structural optimization using the
model-predicted SC Hubbard parameters yields a self-consistent structural-
electronic ground state. The computational cost of these two calculations and
a model evaluation is negligible compared to the DFPT-based SC protocol they
replace. Furthermore, we have demonstrated that the full SC protocol can be
reduced to just a few iterations, which prove adequate for achieving a RMSE
within a few percent of the real value.

Unlike other ML models designed for predicting Hubbard parameters, [28, 61,
29, 62] our model does not rely on experimental data or information from other
state-of-the-art computational methods such as $GW$ and hybrid functionals,
which possess inherent limitations and specific ranges of applicability.
Instead, our model exclusively relies on linear-response theory through DFPT,
which provides material-specific Hubbard parameters directly reflecting the
local chemistry and OS of TM elements. Notably, our model not only predicts
on-site Hubbard $U$ parameters but also inter-site Hubbard $V$ parameters,
crucial for materials characterized by significant covalent interactions, an
aspect that was disregarded in previous studies. [81] Furthermore, the
architecture of our ML model is highly versatile, permitting easy integration
of additional inputs and outputs, facilitating exploration of diverse learning
tasks beyond learning Hubbard parameters. We expect that similar hybrid ML-
accelerated electronic structure methodologies, maintaining accuracy and
transferability, will become prevalent, potentially yielding a comparable
impact on the field as observed with ML interaction potentials.

Lastly, it is essential to highlight the limitations of the trained models we
have presented. These were trained on data generated using a specific
computational setup (see below), which must remain exactly the same when the
model is applied to other systems. In other words, the $U$ and $V$ values
predicted by the model are not transferable across electronic-structure codes
and even across different pseudopotentials within the same code. [82]
Furthermore, the training data only encompass Fe, Mn, and Ni with Hubbard $U$
corrections, while $V$ parameters are available for pairs involving these
elements and O. Consequently, these models can effectively predict Hubbard
parameters for other materials with Hubbard corrections on these atoms or atom
pairs, provided that the model input is generated using identical
pseudopotentials, Hubbard projectors, and functional (further elaborated in
the methods section). Nevertheless, with a more diverse dataset, the models
could be easily extended to accommodate a broader range of compositions.
Indeed, ongoing efforts aim to establish a comprehensive database of Hubbard
parameters for various TM-containing materials, akin to Ref. 83, thereby
significantly broadening the scope of our model to encompass numerous TM
elements across diverse OS and chemical environments. Moreover, our model can
integrate into automated AiiDA workflows, [84, 85] enabling non-experts to
harness it effortlessly and access SC Hubbard parameters with minimal
intervention. Consequently, we believe that our ML model represents a
significant advancement in expediting materials discovery, design, and
understanding based on the DFT+$U$+$V$ approach, thereby unlocking new avenues
for technological progress and breakthroughs.

## METHODS

All calculations are performed using the plane-wave pseudopotential method as
implemented in the Quantum ESPRESSO distribution. [86, 87, 88] We use the
exchange-correlation functional constructed using $\sigma$-GGA with the PBEsol
prescription, [89] and the pseudopotentials are taken from the SSSP library
v1.1 (efficiency). [90, 91] For metallic ground states, we use Gaussian
smearing. To construct the Hubbard projectors, we use atomic orbitals which
are orthonormalized using Löwdin’s method. [92, 93] Structural optimizations
are performed using DFT+$U$+$V$ [58] with the Broyden-Fletcher-Goldfarb-Shanno
(BFGS) algorithm [94] and convergence thresholds for the total energy of
$10^{-6}$ Ry, for forces of $10^{-5}$ Ry/Bohr, and for pressure of $0.5$ Kbar.
The DFPT calculations of Hubbard parameters are performed using the HP code,
[95] with an accuracy of $0.01-0.1$ eV for the computed values of $U$ and $V$.
The information about the kinetic-energy cutoff, $\mathbf{k}$ and $\mathbf{q}$
points sampling of the Brillouin zone for each system are detailed in Sec. S9
in the SI. More technical details can be found directly in the source files
publicly available through the Materials Cloud Archive (see below).

## DATA AVAILABILITY

The data used to produce the results of this work will be available in the
Materials Cloud Archive prior to the publication at
https://archive.materialscloud.org.

## CODE AVAILABILITY

The code for the equivariant neural networks is available at
https://github.com/camml-lab/hubbardml/. Results were generated using the code
at git commit d67007e04676741818f5bcebb4373746351b4b3d.

## ACKNOWLEDGEMENTS

We thank Francesco Aquilante, Mario Geiger, and Jigyasa Nigam for fruitful
discussions. We acknowledge support by the NCCR MARVEL, a National Centre of
Competence in Research, funded by the Swiss National Science Foundation (Grant
number 205602). This work was supported by a grant from the Swiss National
Supercomputing Centre (CSCS) under project ID s1073 (Piz Daint) and ID
465000416 (LUMI-G). This work has been partially supported by MIAI@Grenoble
Alpes, (ANR-19-P3IA-0003).

## Author contributions statement

M.U. developed and trained ML model architecture; A.Z. wrote data processing
code; L.B. and I.T. performed self-consistent calculations of Hubbard
parameters; M.U., A.Z., L.B., N.M. and I.T. analysed the results and wrote the
manuscript.

## References

  * [1] Hohenberg, P. & Kohn, W.  Inhomogeneous electron gas.  _Phys. Rev._ 136, B864, DOI: 10.1103/PhysRev.136.B864 (1964).
  * [2] Kohn, W. & Sham, L.  Self-consistent equations including exchange and correlation effects.  _Phys. Rev._ 140, A1133, DOI: 10.1103/PhysRev.140.A1133 (1965).
  * [3] Perdew, J. & Zunger, A.  Self-interaction correction to density-functional approximations for many-electron systems.  _Phys. Rev. B_ 23, 5048, DOI: 10.1103/PhysRevB.23.5048 (1981).
  * [4] Kulik, H., Cococcioni, M., Scherlis, D. & Marzari, N.  Density Functional Theory in Transition-Metal Chemistry: A Self-Consistent Hubbard U Approach.  _Phys. Rev. Lett._ 97, 103001, DOI: 10.1103/PhysRevLett.97.103001 (2006).
  * [5] Mori-Sánchez, P., Cohen, A. & Yang, W.  Many-electron self-interaction error in approximate density functionals.  _J. Chem. Phys._ 125, 201102, DOI: 10.1063/1.2403848 (2006).
  * [6] Anisimov, V., Zaanen, J. & Andersen, O.  Band theory and Mott insulators: Hubbard $U$ instead of Stoner $I$.  _Phys. Rev. B_ 44, 943, DOI: 10.1103/PhysRevB.44.943 (1991).
  * [7] Liechtenstein, A., Anisimov, V. & Zaanen, J.  Density-functional theory and strong interactions: Orbital ordering in Mott-Hubbard insulators.  _Phys. Rev. B_ 52, R5467, DOI: 10.1103/PhysRevB.52.R5467 (1995).
  * [8] Dudarev, S., Botton, G., Savrasov, S., Humphreys, C. & Sutton, A.  Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+$U$ study.  _Phys. Rev. B_ 57, 1505, DOI: 10.1103/PhysRevB.57.1505 (1998).
  * [9] Campo Jr, V. L. & Cococcioni, M.  Extended DFT+$U$+$V$ method with on-site and inter-site electronic interactions.  _J. Phys.: Condens. Matter_ 22, 055602, DOI: 10.1088/0953-8984/22/5/055602 (2010).
  * [10] Tancogne-Dejean, N. & Rubio, A.  Parameter-free hybridlike functional based on an extended Hubbard model: DFT+U+V.  _Phys. Rev. B_ 102, 155117, DOI: 10.1103/PhysRevB.102.155117 (2020).
  * [11] Lee, S.-H. & Son, Y.-W.  First-principles approach with a pseudohybrid density functional for extended Hubbard interactions.  _Phys. Rev. Research_ 2, 043410, DOI: 10.1103/PhysRevResearch.2.043410 (2020).
  * [12] Sun, J., Ruzsinszky, A. & Perdew, J.  Strongly Constrained and Appropriately Normed Semilocal Density Functional.  _Phys. Rev. Lett._ 115, 036402, DOI: 10.1103/PhysRevLett.115.036402 (2015).
  * [13] Bartók, A. & Yates, J.  Regularized SCAN functional.  _J. Chem. Phys._ 150, 161101, DOI: 10.1063/1.5094646 (2019).
  * [14] Furness, J., Kaplan, A., Ning, J., Perdew, J. & Sun, J.  Accurate and Numerically Efficient r2SCAN Meta-Generalized Gradient Approximation.  _J. Phys. Chem. Lett._ 11, 8208, DOI: 10.1021/acs.jpclett.0c02405 (2020).
  * [15] Adamo, C. & Barone, V.  Toward reliable density functional methods without adjustable parameters: The PBE0 model.  _J. Chem. Phys._ 110, 6158, DOI: 10.1063/1.478522 (1999).
  * [16] Heyd, J., Scuseria, G. & Ernzerhof, M.  Hybrid functionals based on a screened Coulomb potential.  _J. Chem. Phys._ 118, 8207, DOI: 10.1063/1.1564060 (2003).
  * [17] Heyd, J., Scuseria, G. & Ernzerhof, M.  Erratum: “Hybrid functionals based on a screened Coulomb potential” [J. Chem. Phys. 118, 8207 (2003)].  _J. Chem. Phys._ 124, 219906, DOI: 10.1063/1.2204597 (2006).
  * [18] 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_ 1, 033003, DOI: 10.1103/PRXEnergy.1.033003 (2022).
  * [19] 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_ 103, 045141, DOI: 10.1103/PhysRevB.103.045141 (2021).
  * [20] Himmetoglu, B., Floris, A., de Gironcoli, S. & Cococcioni, M.  Hubbard-corrected DFT energy functionals: The LDA+U description of correlated systems.  _Int. J. Quant. Chem._ 114, 14, DOI: 10.1002/qua.24521 (2014).
  * [21] Cococcioni, M. & de Gironcoli, S.  Linear response approach to the calculation of the effective interaction parameters in the LDA+U method.  _Phys. Rev. B_ 71, 035105, DOI: 10.1103/PhysRevB.71.035105 (2005).
  * [22] Wang, L., Maxisch, T. & Ceder, G.  Oxidation energies of transition metal oxides within the GGA+U framework.  _Phys. Rev. B_ 73, 195107, DOI: 10.1103/PhysRevB.73.195107 (2006).
  * [23] Kubaschewski, O., Alcock, C. & Spencer, P.  _Materials thermochemistry_ (Pergamon Press: Elmsford, New York, 1993), 6th edn.
  * [24] Le Bacq, O., Pasturel, A. & Bengone, O.  Impact on electronic correlations on the structural stability, magnetism, and voltage of LiCoPO4 battery.  _Phys. Rev. B_ 69, 245107, DOI: 10.1103/PhysRevB.69.245107 (2004).
  * [25] Aykol, M. & Wolverton, C.  Local environment dependent GGA+$U$ method for accurate thermochemistry of transition metal compounds.  _Phys. Rev. B_ 90, 115105, DOI: 10.1103/PhysRevB.90.115105 (2014).
  * [26] Isaacs, E. & Marianetti, C.  Compositional phase stability of strongly correlated electron materials within DFT+$U$.  _Phys. Rev. B_ 95, 045141, DOI: 10.1103/PhysRevB.95.045141 (2017).
  * [27] Hedin, L.  New Method for Calculating the One-Particle Green’s Function with Application to the Electron-Gas Problem.  _Phys. Rev._ 139, A796, DOI: 10.1103/PhysRev.139.A796 (1965).
  * [28] Yu, M., Yang, S., Wu, C. & Marom, N.  Machine learning the Hubbard $U$ parameter in DFT+$U$ using Bayesian optimization.  _npj Comput. Mater._ 6, 180, DOI: 10.1038/s41524-020-00446-9 (2020).
  * [29] Tavadze, P. _et al._ Exploring DFT+$U$ parameter space with a Bayesian calibration assisted by Markov chain Monte Carlo sampling.  _npj Comput. Mater._ 7, 182, DOI: 10.1038/s41524-021-00651-0 (2021).
  * [30] Skone, J., Govoni, M. & Galli, G.  Self-consistent hybrid functional for condensed systems.  _Phys. Rev. B_ 89, 195112, DOI: 10.1103/PhysRevB.89.195112 (2014).
  * [31] Skone, J., Govoni, M. & Galli, G.  Nonempirical range-separated hybrid functionals for solids and molecules.  _Phys. Rev. B_ 93, 235106, DOI: 10.1103/PhysRevB.93.235106 (2016).
  * [32] Dederichs, P., Blügel, S., Zeller, R. & Akai, H.  Ground States of Constrained Systems: Application to Cerium Impurities.  _Phys. Rev. Lett._ 53, 2512, DOI: 10.1103/PhysRevLett.53.2512 (1984).
  * [33] McMahan, A., Martin, R. & Satpathy, S.  Calculated effective Hamiltonian for La2CuO4 and solution in the impurity Anderson approximation.  _Phys. Rev. B_ 38, 6650, DOI: 10.1103/PhysRevB.38.6650 (1988).
  * [34] 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_ 39, 1708, DOI: 10.1103/PhysRevB.39.1708 (1989).
  * [35] Hybertsen, M., Schlüter, M. & Christensen, N.  Calculation of Coulomb-interaction parameters for La2CuO4 using a constrained-density-functional approach.  _Phys. Rev. B_ 39, 9028, DOI: 10.1103/PhysRevB.39.9028 (1989).
  * [36] Gunnarsson, O.  Calculation of parameters in model Hamiltonians.  _Phys. Rev. B_ 41, 514, DOI: 10.1103/PhysRevB.41.514 (1990).
  * [37] Pickett, W., Erwin, S. & Ethridge, E.  Reformulation of the LDA+$U$ method for a local-orbital basis.  _Phys. Rev. B_ 58, 1201, DOI: 10.1103/PhysRevB.58.1201 (1998).
  * [38] Solovyev, I. & Imada, M.  Screening of Coulomb interactions in transition metals.  _Phys. Rev. B_ 71, 045103, DOI: 10.1103/PhysRevB.71.045103 (2005).
  * [39] 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_ 74, 235113, DOI: 10.1103/PhysRevB.74.235113 (2006).
  * [40] 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_ 93, 085135, DOI: 10.1103/PhysRevB.93.085135 (2016).
  * [41] Mosey, N. & Carter, E.  Ab initio evaluation of Coulomb and exchange parameters for DFT+U calculations.  _Phys. Rev. B_ 76, 155123, DOI: 10.1103/PhysRevB.76.155123 (2007).
  * [42] Mosey, N., Liao, P. & Carter, E.  Rotationally invariant ab initio evaluation of Coulomb and exchange parameters for DFT+U calculations.  _J. Chem. Phys._ 129, 014103, DOI: 10.1063/1.2943142 (2008).
  * [43] 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_ 81, 245103, DOI: 10.1103/PhysRevB.81.245103 (2010).
  * [44] Agapito, L., Curtarolo, S. & Buongiorno Nardelli, M.  Reformulation of DFT+$U$ as a Pseudohybrid Hubbard Density Functional for Accelerated Materials Discovery.  _Phys. Rev. X_ 5, 011006, DOI: 10.1103/PhysRevX.5.011006 (2015).
  * [45] Springer, M. & Aryasetiawan, F.  Frequency-dependent screened interaction in Ni within the random-phase approximation.  _Phys. Rev. B_ 57, 4364, DOI: 10.1103/PhysRevB.57.4364 (1998).
  * [46] 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_ 12, 2413, DOI: 10.1088/0953-8984/12/11/307 (2000).
  * [47] Aryasetiawan, F. _et al._ Frequency-dependent local interactions and low-energy effective models from electronic structure calculations.  _Phys. Rev. B_ 70, 195104, DOI: 10.1103/PhysRevB.70.195104 (2004).
  * [48] Aryasetiawan, F., Karlsson, K., Jepsen, O. & Scönberger, U.  Calculations of Hubbard $U$ from first-principles.  _Phys. Rev. B_ 74, 125106, DOI: 10.1103/PhysRevB.74.125106 (2006).
  * [49] Timrov, I., Marzari, N. & Cococcioni, M.  Hubbard parameters from density-functional perturbation theory.  _Phys. Rev. B_ 98, 085127, DOI: 10.1103/PhysRevB.98.085127 (2018).
  * [50] Timrov, I. _et al._ Electronic structure of Ni-substituted LaFeO3 from near edge x-ray absorption fine structure experiments and first-principles simulations.  _Phys. Rev. Research_ 2, 033265, DOI: 10.1103/PhysRevResearch.2.033265 (2020).
  * [51] Mahajan, R., Kashyap, A. & Timrov, I.  Pivotal Role of Intersite Hubbard Interactions in Fe-Doped $\alpha$-MnO2.  _J. Phys. Chem. C_ 126, 14353, DOI: https://doi.org/10.1021/acs.jpcc.2c04767 (2022).
  * [52] Timrov, I., Kotiuga, M. & Marzari, N.  Unraveling the effects of inter-site Hubbard interactions in spinel Li-ion cathode materials.  _Phys. Chem. Chem. Phys._ 25, 9061, DOI: 10.1039/d3cp00419h (2023).
  * [53] Binci, L., Kotiuga, M., Timrov, I. & Marzari, N.  Hybridization driving distortions and multiferroicity in rare-earth nickelates.  _Phys. Rev. Research_ 5, 033146, DOI: 10.1103/PhysRevResearch.5.033146 (2023).
  * [54] 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.  _Phys. Rev. Materials_ 8, 014007, DOI: 10.1103/PhysRevMaterials.8.014007 (2024).
  * [55] Bonfà, P. _et al._ Magnetostriction-driven muon localisation in an antiferromagnetic oxide.  _Phys. Rev. Lett._ 132, 046701, DOI: 10.1103/PhysRevLett.132.046701 (2024).
  * [56] Gelin, S. _et al._ Ternary oxides of $s-$ and $p-$block metals for photocatalytic solar-to-hydrogen conversion.  _PRX Energy_ 3, 013007, DOI: 10.1103/PRXEnergy.3.013007 (2024).
  * [57] Macke, E., Timrov, I., Marzari, N. & Ciacchi, L. C.  Orbital-resolved DFT+$U$ for molecules and solids.  _arXiv:2312.13580_ .
  * [58] 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_ 102, 235159, DOI: 10.1103/PhysRevB.102.235159 (2020).
  * [59] Kulik, H. & Marzari, N.  Accurate potential energy surfaces with a DFT+$U$(R) approach.  _J. Chem. Phys._ 135, 194105, DOI: 10.1063/1.3660353 (2011).
  * [60] Hsu, H., Umemoto, K., Cococcioni, M. & Wentzcovitch, R.  First-principles study for low-spin LaCoO3 with a structurally consistent Hubbard $U$.  _Phys. Rev. B_ 79, 125124, DOI: 10.1103/PhysRevB.79.125124 (2009).
  * [61] Cai, G. _et al._ Predicting structure-dependent Hubbard U parameters via machine learning.  _Mater. Futures_ 3, 025601, DOI: 10.1088/2752-5724/ad19e2 (2024).
  * [62] Yu, W. _et al._ Active Learning the High-Dimensional Transferable Hubbard $U$ and $V$ Parameters in the DFT+$U$+$V$ Scheme.  _J. Chem. Theory Comput._ 19, 6425, DOI: 10.1021/acs.jctc.2c01116 (2023).
  * [63] Batzner, S. _et al._ E(3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials.  _Nature Communications_ 13, 2453, DOI: 10.1038/s41467-022-29939-5 (2022).
  * [64] Musaelian, A. _et al._ Learning local equivariant representations for large-scale atomistic dynamics.  _Nature Communications_ 14, 579, DOI: 10.1038/s41467-023-36329-y (2023).
  * [65] Bartók, A. P., Kondor, R. & Csányi, G.  On representing chemical environments.  _Physical Review B_ 87, 184115, DOI: 10.1103/PhysRevB.87.184115 (2013).
  * [66] Musil, F. _et al._ Physics-inspired structural representations for molecules and materials.  _Chemical Reviews_ 121, 9759–9815, DOI: 10.1021/acs.chemrev.1c00021 (2021).
  * [67] Uhrin, M.  Through the eyes of a descriptor: Constructing complete, invertible descriptions of atomic environments.  _Phys. Rev. B_ 104, 144110, DOI: 10.1103/PhysRevB.104.144110 (2021).
  * [68] Batatia, I. _et al._ A foundation model for atomistic materials chemistry.  _arXiv:2401.00096_ (2023).
  * [69] Sit, P. H.-L., Car, R., Cohen, M. H. & Selloni, A.  Simple, Unambiguous Theoretical Approach to Oxidation State Determination via First-Principles Calculations.  _Inorg. Chem._ 50, 10259, DOI: 10.1021/ic2013107 (2011).
  * [70] Kondor, R. and Lin, Z. and Trivedi, S.  Clebsch-Gordan Nets: a Fully Fourier Space Spherical Convolutional Neural Network.  _arXiv:1806.09231_ (2018).
  * [71] Thomas, N. Smidt, T. and Kearnes, S. and Yang, L. and Li, L. and Kohlhoff, K. and Riley, P.  Tensor field networks: Rotation- and translation-equivariant neural networks for 3D point clouds.  _arXiv:1802.08219_ (2018).
  * [72] Grisafi, A. and Wilkins, D.M. and Willatt, M.J. and Ceriotti, M.  Atomic-Scale Representation and Statistical Learning of Tensorial Properties.  _ACS Symposium Series_ 1326, 1–21, DOI: 10.1021/bk-2019-1326.ch001 (2019).
  * [73] Geiger, M. _et al._ e3nn/e3nn: 2022-04-13, https://zenodo.org/record/6459381#.YvDuJ-xBxH4, DOI: 10.5281/zenodo.3724963.
  * [74] Merchant, A. _et al._ Scaling deep learning for materials discovery.  _Nature_ 624, 80–85, DOI: 10.1038/s41586-023-06735-9 (2023).
  * [75] Geiger, M. & Smidt, T. E.  e3nn: Euclidean neural networks.  _arXiv:2207.09453_ (2022).
  * [76] Paszke, A. _et al._ Pytorch: An imperative style, high-performance deep learning library.  8024–8035 (Curran Associates, Inc., 2019).
  * [77] Loshchilov, I. & Hutter, F.  Decoupled Weight Decay Regularization.  _arXiv:1711.05101_ (2019).
  * [78] Chakraborty, A., Dixit, M., Aurbach, D. & Major, D.  Predicting accurate cathode properties of layered oxide materials using the SCAN meta-GGA density functional.  _npj Comput. Mater._ 4, 60, DOI: 10.1038/s41524-018-0117-4 (2018).
  * [79] 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. Materials_ 5, 104402, DOI: 10.1103/PhysRevMaterials.5.104402 (2021).
  * [80] Cococcioni, M. & Marzari, N.  Energetics and cathode voltages of Li$M$PO4 olivines ($M$=Fe, Mn) from extended Hubbard functionals.  _Phys. Rev. Materials_ 3, 033801, DOI: 10.1103/PhysRevMaterials.3.033801 (2019).
  * [81] Except Ref. 62 which is based on the fiting of $U$ and $V$ to reproduce the results obtained using hybrid functionals, but it does not provide a ML model that can be used to predict the Hubbard parameters for new materials.
  * [82] Kulik, H. & Marzari, N.  A self-consistent Hubbard U density-functional theory approach to the addition-elimination reactions of hydrocarbons on bare FeO+.  _J. Chem. Phys._ 129, 134314, DOI: 10.1063/1.2987444 (2008).
  * [83] Moore, G. _et al._ High-throughput determination of Hubbard $U$ and Hund $J$ values for transition metal oxides via linear response formalism.  _Phys. Rev. Materials_ 8, 014409 (2024).
  * [84] Huber, S. P. _et al._ AiiDA 1.0, a scalable computational infrastructure for automated reproducible workflows and data provenance.  _Scientific Data_ 7, 300, DOI: 10.1038/s41597-020-00638-4 (2020).
  * [85] 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.  _Computational Materials Science_ 187, 110086, DOI: 10.1016/j.commatsci.2020.110086 (2021).
  * [86] Giannozzi, P. _et al._ Quantum ESPRESSO: A modular and open-source software project for quantum simulations of materials.  _J. Phys.: Condens. Matter._ 21, 395502, DOI: 10.1088/0953-8984/21/39/395502 (2009).
  * [87] Giannozzi, P. _et al._ Advanced capabilities for materials modelling with Quantum ESPRESSO.  _J. Phys.: Condens. Matter._ 29, 465901, DOI: 10.1088/1361-648X/aa8f79 (2017).
  * [88] Giannozzi, P. _et al._ Quantum ESPRESSO toward the exascale.  _J. Chem. Phys._ 152, 154105, DOI: 10.1063/5.0005082 (2020).
  * [89] Perdew, J. _et al._ Restoring the Density-Gradient Expansion for Exchange in Solids and Surfaces.  _Phys. Rev. Lett._ 100, 136406, DOI: 10.1103/PhysRevLett.100.136406 (2008).
  * [90] Prandini, G., Marrazzo, A., Castelli, I. E., Mounet, N. & Marzari, N.  Precision and efficiency in solid-state pseudopotential calculations.  _npj Computational Materials_ 4, 72, DOI: 10.1038/s41524-018-0127-2 (2018).
  * [91] The SSSP library of the Materials Cloud: https://www.materialscloud.org/discover/sssp/table/efficiency.
  * [92] Löwdin, P.-O.  On the Non-Orthogonality Problem Connected with the Use of Atomic Wave Functions in the Theory of Molecules and Crystals.  _J. Chem. Phys._ 18, 365, DOI: 10.1063/1.1747632 (1950).
  * [93] Mayer, I.  On Löwdin’s method of symmetric orthogonalization.  _Int. J. Quant. Chem._ 90, 63, DOI: 10.1002/qua.981 (2002).
  * [94] Fletcher, R.  _Practical Methods of Optimization_ (Wiley, Chichester, 1987), 2nd edn.
  * [95] Timrov, I., Marzari, N. & Cococcioni, M.  HP – A code for the calculation of Hubbard parameters using density-functional perturbation theory.  _Comput. Phys. Commun._ 279, 108455, DOI: 10.24435/materialscloud:v6-zd (2022).