# DFT+$U$ within the framework of linear combination of numerical atomic
orbitals

Xin Qu  Rocket Force University of Engineering, Xi’an, Shaanxi 710025, China
CAS Key Laboratory of Quantum Information, University of Science and
Technology of China, Hefei, Anhui 230026, China  Peng Xu  Rocket Force
University of Engineering, Xi’an, Shaanxi 710025, China  Hong Jiang
<EMAIL_ADDRESS>Beijing National Laboratory for Molecular Sciences,
College of Chemistry and Molecular Engineering, Peking University, 100871
Beijing, China  Lixin He<EMAIL_ADDRESS>CAS Key Laboratory of Quantum
Information, University of Science and Technology of China, Hefei, Anhui
230026, China  Xinguo Ren<EMAIL_ADDRESS>Institute of Physics, Chinese
Academy of Sciences, Beijing 100190, China  Songshan Lake Materials
Laboratory, Dongguan 523808, Guangdong, China

###### Abstract

We present a formulation and implementation of the DFT+U method within the
framework of linear combination of numerical atomic orbitals (NAO). Our
implementation not only enables single-point total energy and electronic-
structure calculations but also provides access to atomic forces and stresses,
hence allowing for full structure relaxations of periodic systems.
Furthermore, our implementation allows one to deal with non-collinear spin
texture, with the spin-orbit coupling (SOC) effect treated self-consistently.
The key aspect behind our implementation is a suitable definition of the
correlated subspace when multiple atomic orbitals with the same angular
momentum are used, and this is addressed via the “Mullken charge projector”
constructed in terms of the first (most localized) atomic orbital within the
$d/f$ angular momentum channel. The important Hubbard $U$ and Hund $J$
parameters can be estimated from a screened Coulomb potential of the Yukawa
type, with the screening parameter either chosen semi-empirically or
determined from the Thomas-Fermi screening model. Benchmark calculations are
performed for four late transition metal monoxide bulk systems, i.e., MnO,
FeO, CoO, and NiO, and for the 5$d$-electron compounds IrO2. For the former
type of systems, we check the performance of our DFT+$U$ implementation for
calculating band gaps, magnetic moments, electronic band structures, as well
as forces and stresses; for the latter, the efficacy of our DFT+$U$+SOC
implementation is assessed. Systematic comparisons with available experimental
results, and especially with the results from other implementation schemes are
carried out, which demonstrate the validity of our NAO-based DFT+$U$ formalism
and implementation.

##  I Introduction

The density functional theory (DFT) developed by Hohenberg, Kohn and Sham [1,
2] uses the charge density as the basic variable to determine the ground state
of interacting many-particle systems, which allows one to simulate systems of
sizes that are prohibitively expensive for wave function based methods [3].
Within DFT, suitable approximations can be designed to achieve excellent
balance between accuracy and efficiency, and as such DFT has been widely used
for electronic structure calculations for a large variety of physical and
chemical problems. However, within the popular local (spin-) density
approximation (L(S)DA) or generalized gradient approximations (GGAs), DFT
encounters failures in several known situations – in particular for strongly
correlated materials, usually characterized by partially filled d/f electron
states. These include transition metals (TM) and their oxides, rare-earth
compounds, and lanthanides, to name a few, where L(S)DA/GGAs typically yield
quantitatively or even qualitatively wrong results.

The failure of L(S)DA or GGAs in strongly correlated materials can be traced
back to the large delocalization error [4, 5, 6, 7], more often termed as
many-body self-interaction errors [8], and the static correlation error [4, 9,
10] associated with these functionals. Many approaches have been proposed to
address these deficiencies, such as self-interaction corrected DFT [11],
hybrid functionals [12, 13], the localized orbital scaling correction [5, 7],
fractional spin correction [6], and so on. Among these, the most popular
approaches in solid-state physics are the combination of L(S)DA and GGAs with
the non-perturbative many-body technique – dynamic mean-field theory (DMFT)
[14, 15, 16, 17, 18, 19] and with a simpler mean-field-type correction based
on the Hubbard model [20, 21, 22]. The latter approach, commonly known as
DFT+U [23, 24, 25], inherits the efficiency of L(S)DA/GGA, but gains the
strength of the Hubbard model in describing the physics of strongly correlated
systems. Owing to its success in describing certain non-trivial properties
arising from strong correlation, e.g. band-gap opening in transition metal
oxides (TMO) and/or rare-earth compounds at a similar cost as L(S)DA/GGAs,
DFT+U has become one of the most widely used first-principles approaches for
strongly correlated (especially insulating) systems.

Being enormously successful in dealing with open-shell systems with partially
filled $d$/$f$ states, DFT+U has become a standard module available in many
DFT code packages, based on various basis-set frameworks. For instance, it has
been implemented within the full-potential all-electron linearized augmented
plane-wave (LAPW) framework, as exemplified by the WIEN2k [26] and ELK [27]
codes, and within the projector-augmented-wave (PAW) or norm-conserving
pseudopotential based plane-wave framework, with VASP [28], Abinit [29], and
QUANTUM ESPRESSO [30, 31] as prominent examples. In recent years, the linear
combination of numerical atomic orbitals (NAOs) has emerged as a versatile
basis set framework for implementing both conventional local and semi-local
density functional approximations [32, 33, 34, 35, 36, 37], hybrid functionals
[38, 39, 40, 41], and many-body perturbation theories [38, 42, 43]. Compared
to other basis set choices, NAOs are considered to be advantageous for
simulating large-scale systems, due to their compact size and strict locality
in real space. Thus the implementation of the DFT+U method within the NAO
framework is of great interest. Similar to other numerical schemes, the key
aspect for the DFT+U implementation is to define a suitable projector, which
maps the full-orbital space where the first-principles DFT calculations are
performed to a local, correlated subspace where the local orbitals behind the
(generalized) Hubbard model are defined. The NAO-based DFT+$U$ implementations
have been reported for SIESTA [34, 44], OpenMX [45], and FHI-aims codes [46],
where the emphasis has been placed on the choice of the projector function, as
well as their influence on the suitable $U$ values and the obtained results.
Experiences strongly suggest that the DFT+$U$ cannot be taken as a black-box
method and a proper use of this approach requires a good understanding of the
underlying technical aspects of a given implementation.

In this work, we report yet another implementation of the DFT+$U$ method
within ABACUS [37], which is a first-principles code package based on the
norm-conserving pseudopotentials and NAO basis functions. ABACUS allows one to
use plane-wave basis functions as well, but our current DFT+$U$ implementation
is based on the NAO basis set framework, and hence shares similarities with
the OpenMX [45] and FHI-aims implementations [46]. However, due to the
different strategies adopted for basis set generation, the resultant NAOs
differ in shapes and spatial extent, which further affects the projector and
other implementation details. Our implementation not only supports usual self-
consistent DFT+$U$ electronic structure calculations but also allows for force
and stress computations, thus enabling full structural relaxations.
Furthermore, non-collinear spin configurations and SOC effects can be treated
within our implementation, which is instrumental for systems containing heavy
elements. Last but not least, we have made attempts to compute the Hubbard $U$
value on the fly from a Yukawa-type screened Coulomb potential. The
performance of such a scheme for determining the $U$ value will be examined.
Considering all these aspects, we believe that a concise description of the
formulation and numerical details of our implementation should be not only
useful for the users and developers of the ABACUS code, but also of general
interest to the electronic-structure community using the NAO basis sets.

The paper is organized as follows. In Sec. II we present the detailed
formulation behind our implementation, including the rotationally invariant
DFT+$U$ energy functional, our choice of the local projector, the force and
stress evaluations, as well as the incorporation of SOC in the NAO-based
DFT+$U$ formalism and the determination of the $U$, $J$ values from the Yukawa
screened Coulomb potential. Sec. III presents the computational details in
this paper. In Sec IV, the efficacy of our formalism and the validity of our
implementation will be examined. This is done by comparing the results of our
implementation to those of experiments and particularly of other well-tested
codes. Finally, we conclude this work in Sec. V.

##  II Formalism

###  II.1 General DFT+U functional

The basic idea of DFT+U is to treat strongly correlated d/f electrons in terms
of multi-orbital Hubbard model at the level of static Hartree-Fock mean-field
theory, whereas all the rest electrons are described at the level of
conventional density functional approximations (DFAs) like L(S)DA and GGA.
Since the $U$ correction term in DFT+U corresponds to a Hartree-Fock
approximation of the multi-orbital Hubbard model, the unphysical self-
interactions among strongly correlated d/f electrons present in the Hartree
energy are canceled out. This is considered as the main reason responsible for
the success of DFT+U in producing more reliable insulating gaps, magnetic
moments, and other properties for TMOs [25]. The many-body complexity is
encoded into the screened Coulomb interaction among strongly correlated d/f
electrons, usually parameterized in terms of the Hubbard U for describing the
direct Coulomb interaction and the Hund J for describing the exchange
interaction.

The formulation of the DFT+U approach begins with the following energy
functional

$E_{\text{DFT}+U}=E_{\text{DFA}}+E_{U}-E_{\text{dc}}~{},$ (1)

whereby $E_{\text{DFA}}$ is the energy of density functional approximations at
the level of L(S)DA or GGA, and $E_{U}$ is the Coulomb interaction energy due
to strongly correlated electrons given by the Hartree-Fock approximation to
multi-orbital Hubbard model. The double counting term $E_{\text{dc}}$ is
subtracted here to discount the Coulomb interaction energy that is already
included in DFAs at an average level.

Now we discuss the last two terms on the right-hand side of Eq. (1). Using the
second-quantization language, the full electron-electron interaction term can
be written as

$\begin{aligned}
\hat{H}_{\text{Hub}}=&\frac{1}{2}\sum_{\\{m\\}}\sum_{\sigma\sigma^{\prime}}\langle
mm^{\prime}\left|v_{sc}\right|m^{\prime\prime}m^{\prime\prime\prime}\rangle\\\
&\times\hat{c}^{\sigma\dagger}_{m}\hat{c}^{\sigma^{\prime}\dagger}_{m^{\prime}}\hat{c}_{m^{\prime\prime\prime}}^{\sigma^{\prime}}\hat{c}_{m^{\prime\prime}}^{\sigma}\end{aligned}~{},$
(2)

where the Coulomb interaction matrix elements are

$\displaystyle\langle
mm^{\prime}\left|v_{sc}\right|m^{\prime\prime}m^{\prime\prime\prime}\rangle$
(3) $\displaystyle=$ $\displaystyle\int d\textbf{r}\int
d\textbf{r}^{\prime}\varphi_{m}^{*}\left(\textbf{r}\right)\varphi_{m^{\prime}}^{*}\left(\textbf{r}^{\prime}\right)v_{sc}\left(\textbf{r},\textbf{r}^{\prime}\right)\varphi_{m^{\prime\prime}}\left(\textbf{r}\right)\varphi_{m^{\prime\prime\prime}}\left(\textbf{r}^{\prime}\right)~{}.$

In Eqs. (2) and (3), $v_{sc}({\bf r},{\bf r^{\prime}})$ is the (statically)
screened Coulomb potential, $\\{m\\}$ the local orbital indices for d or f
subshell, $\sigma$ the spin index and $\hat{c}_{m}^{\sigma\dagger}$ and
$\hat{c}_{m}^{\sigma}$ the creation and annihilation operators associated with
the local correlated orbitals.

For simplicity, we neglect the spin-orbit coupling (SOC) effect at this stage.
The generalized DFT+$U$ formalism that incorporates the SOC effect will be
presented in Sec. II.4. In this case, each Kohn-Sham (KS) spin orbital is a
product of a spatial function $\psi_{n\textbf{k}}\left(\textbf{r}\right)$ and
a spin function. The ground state of the KS system is a Slater determinant
formed by occupied KS spin-orbitals [47]. The Hartree-Fock approximation to
the Hubbard Hamiltonian Eq. (2) can be obtained by evaluating its expectation
value within the KS ground state $|0\rangle$, yielding the energy contribution
as

$\displaystyle E_{U}=$ $\displaystyle\langle 0|\hat{H}_{\text{Hub}}|0\rangle$
(4) $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{\\{m\\},\sigma}\left\\{\left(\left\langle
mm^{\prime}\left|v_{sc}\right|m^{\prime\prime}m^{\prime\prime\prime}\right\rangle\right.\right.$
$\displaystyle\left.-\left\langle
mm^{\prime}\left|v_{\mathrm{sc}}\right|m^{\prime\prime\prime}m^{\prime\prime}\rangle\right)\right.n_{m^{\prime\prime}m}^{\sigma}n_{m^{\prime\prime\prime}m^{\prime}}^{\sigma}$
$\displaystyle\left.+\left\langle
mm^{\prime}\left|v_{\mathrm{sc}}\right|m^{\prime\prime}m^{\prime\prime\prime}\right\rangle
n_{m^{\prime\prime}m}^{\sigma}n_{m^{\prime\prime\prime}m^{\prime}}^{-\sigma}\right\\}\,.$

In Eq. (4), $n_{mm^{\prime}}^{\sigma}$ is the local occupation matrix given by

$\displaystyle n_{mm^{\prime}}^{\sigma}$ $\displaystyle=\langle
0|\hat{c}^{\sigma\dagger}_{m^{\prime}}\hat{c}^{\sigma}_{m}|0\rangle$ (5)
$\displaystyle=\frac{1}{N_{\textbf{k}}}\sum_{n\textbf{k}}f_{n\textbf{k}}^{\sigma}\langle\psi_{n\textbf{k}}^{\sigma}\sigma|m\sigma\rangle\langle
m^{\prime}\sigma|\psi_{n\textbf{k}}^{\sigma}\sigma\rangle\,$

with $f_{n\textbf{k}}^{\sigma}$ being the occupation number of KS orbitals
$|\psi_{n\textbf{k}}^{\sigma}\sigma\rangle$. Here
$|\psi_{n\textbf{k}}^{\sigma}\sigma\rangle=|\psi_{n\textbf{k}}^{\sigma}\rangle|\sigma\rangle$
is a product of the spatial wavevector $|\psi_{n\textbf{k}}^{\sigma}\rangle$,
which is the Kohn-Sham wavefunction of the $\sigma$-spin component in the
position space, and the spin function $|\sigma\rangle$. Similarly the local
spin orbital $|m\sigma\rangle=|\phi_{m}\rangle|\sigma\rangle$. Furthermore,
$N_{\textbf{k}}$ is the number of ${\bf k}$ points in the Brillouin zone (BZ),
which equals the number of unit cells in the Born-Von-Kármén (BvK) supercell
under the periodic boundary condition. Equation (4) is the well known
rotationally invariant form of DFT+U firstly proposed by Lichtenstein et al.
[48] in 1995, whereby the local occupation number matrix
$n_{mm^{\prime}}^{\sigma}$ is the key quantity.

Since the local occupation matrix is symmetric, one can always introduce an
unitary transformation to diagonalize it and arrives at

$E_{U}=\frac{1}{2}U\sum_{mm^{\prime},\sigma}\textbf{n}_{m}^{\sigma}\textbf{n}_{m^{\prime}}^{-\sigma}+\frac{1}{2}\left(U-J\right)\sum_{m\neq
m^{\prime},\sigma}\textbf{n}_{m}^{\sigma}\textbf{n}_{m^{\prime}}^{\sigma}.$
(6)

Here we use bold n to denote the vector comprising the eigenvalues of the
local occupation matrix. Furthermore, $U\equiv\langle
mm^{\prime}\left|v_{sc}\right|mm^{\prime}\rangle$ and $J\equiv\langle
mm^{\prime}\left|v_{sc}\right|m^{\prime}m\rangle$ are the direct Coulomb and
exchange integrals of the electrons in the correlated subspace, respectively,
which are assumed to be isotropic, i.e., independent of the magnetic quantum
number $m$. Note that the self-interaction is absent both in Eq. (4) and (6).
Theoretically, $\langle mm^{\prime}\left|v_{sc}\right|mm^{\prime}\rangle$ can
be evaluated through Slater integrals, but the detailed from of the screened
Coulomb interaction $v_{\text{sc}}$ remains unknown, and therefore in
practical calculations $U$ and $J$ are most commonly treated as adjustable
parameters or obtained via pragmatic schemes like constrained DFT [49, 50, 51,
23], constrained random-phase approximation (RPA) [52, 53, 54, 55], or linear-
response approach [31].

The double-counting term $E_{\text{dc}}$ in Eq. (1) is an important portion of
the DFT+$U$ theory and needs to be properly treated. Unfortunately, there are
uncertainties for a rigorous definition of this term. This difficulty arises
from by the fact that local/semi-local DFAs are not orbital-resolved theories,
and contributions from individual orbitals cannot be separated from one
another. By now there are two main double counting schemes used in practical
DFT+U calculations. One is the so-called “around mean field (AMF)” scheme and
another is the “fully localized limit (FLL)” scheme [23, 24, 56]. Both schemes
are physically motivated. It’s generally accepted that the former gives a
better description of metallic systems while the latter one is more suitable
for insulating systems [10]. The FLL double counting term is given by

$E_{\text{dc}}=\frac{1}{2}UN\left(N-1\right)-\frac{1}{2}J\sum_{\sigma}N^{\sigma}\left(N^{\sigma}-1\right),$
(7)

where $N=\sum_{\sigma}N^{\sigma}$, and $N^{\sigma}$ is the total number of
correlated d or f electrons of spin $\sigma$. The FLL double counting term can
be derived by assuming integer occupations of correlated d or f electrons in
the atomic limit. In our implementation, the FLL scheme is used.

Subtracting the double counting term $E_{\text{dc}}$ from $E_{U}$ and making
some simple derivation, the DFT+U energy correction can be explicitly
expressed as

$\Delta
E_{\text{DFT+}U}=E_{U}-E_{\text{dc}}=\frac{1}{2}\left(U-J\right)\sum_{m\sigma}\left(\textbf{n}_{m}^{\sigma}-\textbf{n}_{m}^{\sigma}\textbf{n}_{m}^{\sigma}\right)~{}.$
(8)

Since the trace of an arbitrary matrix remains unchanged after the unitary
transform, $\Delta E_{\text{DFT+}U}$ can also be rewritten as

$\displaystyle\Delta
E_{\text{DFT+}U}=\frac{1}{2}\left(U-J\right)\sum_{\sigma}\left[\sum_{m}n_{mm}^{\sigma}-\sum_{mm^{\prime}}n_{mm^{\prime}}^{\sigma}n_{m^{\prime}m}^{\sigma}\right].$
(9)

The above energy correction functionals, Eq. (8) and (9), are the simplified
form of the rotationally invariant scheme proposed by Dudarev et al. [57],
while still retaining the rotational invariance as the energy correction stays
unchanged under unitary transformations of the given set of correlated
orbitals. Within this functional, the DFT+U total energy reduces to standard
L(S)DA/GGA in the case of empty or full (0 or 1) occupation of local orbitals.

###  II.2 Mulliken charge projector

The DFT+$U$ formalism presented in Sec. II.1 only applies to the single-site
case, i.e., only one correlated atom in the cell. To deal with the multi-site
cases, it is necessary to introduce an extra correlated atomic index $I$ to
label the local occupation matrix and the parameters U and J. Furthermore, in
Sec. II.1 it is implicitly assumed that the local correlated orbitals
$\\{|m\rangle\\}$ are orthonormal to each other, and these are not necessarily
satisfied for practically chosen local orbitals. Considering these
complexities, it is convenient to introduce a local projection operator
$P_{mm^{\prime}}^{\sigma}$, called “projector”, such that the spin-dependent
local occupation matrix is given by

$n_{I,mm^{\prime}}^{\sigma}=\frac{1}{N_{\textbf{k}}}\sum_{n\textbf{k}}f_{n\textbf{k}}^{\sigma}\langle\psi_{n\textbf{k}}^{\sigma}\sigma|\hat{P}_{I,mm^{\prime}}^{\sigma}|\psi_{n\textbf{k}}^{\sigma}\sigma\rangle\,.$
(10)

where $I$ denotes a correlated atom to which the $U$ correction needs to be
applied. The key issue in the implementation of DFT+U is to construct such a
projector that maps the full Kohn-Sham orbital space into the correlated
subspace. The choice of the projector depends on the underlying computational
frameworks, ranging from the linear muffin-tin orbital method [25], the LAPW
method [26], to the PAW method [29, 28] and pseudopotential-based plane-wave
method [31]. Within the NAO-based framework, the most straightforward way is
to utilize the local $d/f$-type atomic-orbital basis functions to construct
the projector. However, in practical calculations the NAOs centering on
neighboring atoms have finite overlaps, i.e., they are non-orthogonal to each
other. This non-orthogonality has to be taken into account when defining a
suitable projector. In this regard, we follow the previous work of Han et al.
[45] where the so-called “Mulliken charge projector” is used. This projector
has the nice property that the sum rule is satisfied, in the sense that the
total electronic charges are conserved when summing up partial charges over
all projected channels.

Specifically, one needs to define dual orbitals associated with the original
atomic orbitals as

$\tilde{\phi}_{\textbf{k}\mu}({\bf r})=\sum_{\nu}\phi_{\textbf{k}\nu}({\bf
r})S_{\nu\mu}^{-1}\left(\textbf{k}\right),$ (11)

where

$\phi_{\textbf{k}\nu}({\bf
r})=\sum_{\textbf{R}}e^{i\textbf{k}\textbf{R}}\phi_{\nu}^{{\bf R}}\left({\bf
r}\right)$ (12)

is the Bloch summation of NAOs. And $\phi_{\nu}^{{\bf R}}\left({\bf
r}\right)=\phi_{\nu}\left(\textbf{r}-{\bm{\tau}}_{a}-\textbf{R}\right)$
denotes a NAO centering on the a-th atom within the unit cell R. The orbital
indices $\mu$ and $\nu$ are a combination of $\\{a,l,\zeta,m\\}$, with $a$
labeling the atomic site, $l$, $m$ the angular and magnetic momentum, and
$\zeta$ the multiplicity (different radial functions) for a given $l$,
respectively. Furthermore, $S_{\mu\nu}\left(\textbf{k}\right)$ is the overlap
matrix in reciprocal space [37]

$S_{\mu\nu}\left(\textbf{k}\right)=\langle\phi_{\textbf{k}\mu}|\phi_{\textbf{k}\nu}\rangle=\sum_{\textbf{R}}e^{-i\textbf{k}\textbf{R}}\langle\phi_{\mu}^{{\bf
R}}|\phi_{\nu}^{\textbf{0}}\rangle\,.$ (13)

Here we follow the convention that the extended Bloch orbitals are normalized
within the BvK supercell cell and the real-space integration indicated by the
braket goes over the supercell cell. It can be readily shown that the dual and
original Bloch orbitals satisfy the following biorthogonality relation

$\langle\phi_{\textbf{q}\mu}|\tilde{\phi}_{\textbf{k}\nu}\rangle=\delta_{\textbf{k}\textbf{q}}\delta_{\mu\nu}\,.$
(14)

Making use of the dual orbitals, we define the projector used in the present
work as

$\begin{aligned}
\hat{P}_{I,mm^{\prime}}^{\sigma}=&\frac{1}{4N_{\textbf{k}}}\sum_{\textbf{k}}\left(|\tilde{\phi}_{\textbf{k},\beta
m^{\prime}}\sigma\rangle\langle\phi_{\textbf{k},\beta m}\sigma|\right.\\\
&+|\phi_{\textbf{k},\beta
m^{\prime}}\sigma\rangle\langle\tilde{\phi}_{\textbf{k},\beta
m}\sigma|+|\tilde{\phi}_{\textbf{k},\beta
m}\sigma\rangle\langle\phi_{\textbf{k},\beta m^{\prime}}\sigma|\\\
&+\left.|\phi_{\textbf{k},\beta
m}\sigma\rangle\langle\tilde{\phi}_{\textbf{k},\beta
m^{\prime}}\sigma|\right)\end{aligned}~{},$ (15)

whereby the index $\beta$ groups together the indices $\\{I,l,\zeta\\}$ with
$l$ and $\zeta$ belonging to the correlated channel of the correlated atom
$I$. Note that this projector is slightly different from that introduced in
the work of Han et al [45], which is essentially an average of the first two
terms in Eq. (15). The local occupation matrices yielded by the projector used
in Ref. [45] are Hermitian but not necessarily real symmetric. For convenience
and numerical simplicity, we symmetrize the projector as is done in Eq. (15),
and then the resultant local occupation matrix is guaranteed to be real and
symmetric. This is consistent with the feature that the “on-site” global KS
density matrix is also real and symmetric. Inserting Eq. (15) into Eq. (10),
we arrive at

$\displaystyle n_{I,mm^{\prime}}^{\sigma}=$
$\displaystyle\frac{1}{4N_{\textbf{k}}}\sum_{\textbf{k}}\sum_{\mu}\left(S_{\beta
m,\mu}\left(\textbf{k}\right)\rho_{\mu,\beta
m^{\prime}}^{\sigma}\left(\textbf{k}\right)\right.$ (16) $\displaystyle+$
$\displaystyle\left.\rho_{\beta
m,\mu}^{\sigma}\left(\textbf{k}\right)S_{\mu,\beta
m^{\prime}}\left(\textbf{k}\right)+S_{\beta
m^{\prime},\mu}\left(\textbf{k}\right)\rho_{\mu,\beta
m}^{\sigma}\left(\textbf{k}\right)\right.$ $\displaystyle+$
$\displaystyle\left.\rho_{\beta
m^{\prime},\mu}^{\sigma}\left(\textbf{k}\right)S_{\mu,\beta
m}\left(\textbf{k}\right)\right)\,,$

where $\rho_{\mu\nu}^{\sigma}\left(\textbf{k}\right)$ is the spin-dependent KS
density matrix

$\rho_{\mu\nu}^{\sigma}\left(\textbf{k}\right)=\sum_{n}f_{n\textbf{k}}^{\sigma}c_{n\textbf{k},\mu}^{\sigma}c_{n\textbf{k},\nu}^{\sigma*}$
(17)

with $c_{n\textbf{k},\mu}^{\sigma}$ being the KS eigenvectors, satisfying the
generalized orthogonality relationship,

$\sum_{\mu,\nu}c_{n\textbf{k},\mu}^{\sigma}S_{\mu\nu}({\bf
k})c_{n^{\prime}\textbf{k},\nu}^{\sigma\ast}=\delta_{nn^{\prime}}\,.$ (18)

If we project the KS density matrix to all local atomic orbital channels and
sum the traces of the resultant local occupation matrices up, i.e., requiring
that $\beta,m$ go over all the basis indices $\\{a,l,\zeta,m\\}$, one then
obtains

$\displaystyle\sum_{\sigma,\nu}n_{\nu\nu}^{\sigma}$ (19)
$\displaystyle=\frac{1}{4N_{\textbf{k}}}\sum_{\sigma}\sum_{\textbf{k}}\sum_{\mu\nu}\left(S_{\nu,\mu}\left(\textbf{k}\right)\rho_{\mu,\nu}^{\sigma}\left(\textbf{k}\right)+\rho_{\nu,\mu}^{\sigma}\left(\textbf{k}\right)S_{\mu,\nu}\left(\textbf{k}\right)\right.$
$\displaystyle~{}~{}~{}~{}\quad\left.+S_{\nu,\mu}\left(\textbf{k}\right)\rho_{\mu,\nu}^{\sigma}\left(\textbf{k}\right)+\rho_{\nu,\mu}^{\sigma}\left(\textbf{k}\right)S_{\mu,\nu}\left(\textbf{k}\right)\right)$
$\displaystyle=\frac{1}{N_{\bf k}}\sum_{n,\sigma}\sum_{{\bf k}}f_{n{\bf
k}}^{\sigma}$ $\displaystyle=N_{e}\,,$

where $N_{e}$ is the total number of electrons in one unit cell. In deriving
Eq. (19), the orthogonality relationship Eq. (18) is used. Equation (19) is
the above-mentioned sum rule satisfied by the “Mulliken charge operator”.

###  II.3 Effective potential, force and stress

To perform self-consistent DFT+$U$ calculations, and to enable structure
relaxations, one needs to derive the expressions of the effective single-
particle potential and the forces and stresses corresponding to the DFT+$U$
energy functional. To this end, we first generalize the single-site DFT+$U$
energy correction as given by Eq. (9) to the multi-site case,

$\Delta
E_{\text{DFT+}U}=\frac{1}{2}\sum_{I}\bar{U}_{I}\sum_{\sigma}\\{\text{Tr}\left(n_{I}^{\sigma}\right)-\text{Tr}\left(n_{I}^{\sigma}n_{I}^{\sigma}\right)\\}\,,$
(20)

where $\bar{U}_{I}=U_{I}-J_{I}$ is the effective interaction parameter on the
correlated atom $I$. Again the isotropy of the interaction parameters is
assumed. The contribution of the energy correction $\Delta E_{\text{DFT+}U}$
to the KS effective potential operator is given by its derivative with respect
to the ${\bf k}$-dependent density matrix operator,

$\hat{\rho}_{\textbf{k}}^{\sigma}=\frac{1}{N_{\bf
k}}\sum_{n}f^{\sigma}_{n\textbf{k}}|\psi_{n\textbf{k}}^{\sigma}\sigma\rangle\langle\psi_{n\textbf{k}}^{\sigma,\dagger}\sigma|~{}.$
(21)

That is

$\begin{aligned}
\Delta\hat{V}_{\text{eff}}^{\sigma}\left(\textbf{k}\right)&=\frac{\delta\Delta
E_{\text{DFT+}U}}{\delta\hat{\rho}_{\textbf{k}}^{\sigma}}\\\
&=\sum_{I}\sum_{mm^{\prime}}\Delta
V_{I,mm^{\prime}}^{\sigma}\hat{P}_{I,mm^{\prime}}^{\sigma}\left(\textbf{k}\right)\end{aligned}~{},$
(22)

where

$\Delta
V_{I,mm^{\prime}}^{\sigma}=\bar{U}_{I}\left(1/2\delta_{mm^{\prime}}-n_{I,mm^{\prime}}^{\sigma}\right)$
(23)

is the correction to the effective single-particle potential in the local
subspace, and

$\displaystyle\hat{P}_{I,mm^{\prime}}^{\sigma}\left(\textbf{k}\right)=$
$\displaystyle\frac{1}{4}\left(|\tilde{\phi}_{\textbf{k},\beta
m^{\prime}}\sigma\rangle\langle\phi_{\textbf{k},\beta
m}\sigma|+|\phi_{\textbf{k},\beta
m^{\prime}}^{\sigma}\rangle\langle\tilde{\phi}_{\textbf{k},\beta
m}\sigma|\right.$ (24) $\displaystyle\left.+|\tilde{\phi}_{\textbf{k},\beta
m}\sigma\rangle\langle\phi_{\textbf{k},\beta
m^{\prime}}\sigma|+|\phi_{\textbf{k},\beta
m}\sigma\rangle\langle\tilde{\phi}_{\textbf{k},\beta
m^{\prime}}\sigma|\right)\,.$

is the ${\bf k}$-dependent Mulliken projector. The matrix form of the
effective potential within the full NAO basis set is given by

$\displaystyle\Delta V_{\text{eff},\mu\nu}^{\sigma}\left(\textbf{k}\right)=$
$\displaystyle\langle\phi_{{\bf
k}\mu}|\Delta\hat{V}_{\text{eff}}^{\sigma}\left(\textbf{k}\right)|\phi_{{\bf
k}\nu}\rangle$ (25) $\displaystyle=$
$\displaystyle\sum_{\textbf{R}}e^{-i\textbf{k}\textbf{R}}\langle\phi_{\mu}^{\sigma}\left(\textbf{R}\right)|\Delta\hat{V}_{\text{eff}}^{\sigma}\left(\textbf{k}\right)|\phi_{\nu}^{\sigma}\left(\textbf{0}\right)\rangle$
$\displaystyle=$ $\displaystyle\frac{1}{4}\sum_{I}\sum_{mm^{\prime}}\Delta
V_{I,mm^{\prime}}^{\sigma}\left\\{S_{\beta
m,\nu}\left(\textbf{k}\right)\delta_{\mu,\beta m^{\prime}}\right.$
$\displaystyle+S_{\mu,\beta m^{\prime}}\left(\textbf{k}\right)\delta_{\beta
m,\nu}+S_{\beta m^{\prime},\nu}\left(\textbf{k}\right)\delta_{\mu,\beta m}$
$\displaystyle+\left.S_{\mu,\beta m}\left(\textbf{k}\right)\delta_{\beta
m^{\prime},\nu}\right\\}\,$

which is to be added to DFA Hamiltonian matrix to obtain the DFT+$U$ one. As
shown in Eq. (23), in the case of diagonal half-integer occupations the DFT+U
Hamiltonian reduces to standard DFAs.

The contribution of the energy correction $\Delta E_{\text{DFT+}U}$ to the
force on the $a$-th atom can be evaluated by its derivative with respect to
the atomic coordinate ${\bm{\tau}}_{a}$

$\displaystyle F_{a}=$ $\displaystyle\frac{d\Delta
E_{\text{DFT+}U}}{d{\bm{\bm{\tau}}}_{a}}=\sum_{\sigma}\sum_{I}\sum_{mm^{\prime}}\Delta
V_{I,mm^{\prime}}^{\sigma}\frac{dn_{I,mm^{\prime}}^{\sigma}}{d{\bm{\tau}}_{a}},$
(26)

which implies that the force due to the DFT+$U$ energy correction stems
entirely from the change of the local occupation matrix in response to the
atomic displacement. According to Eq. (16), the change of the local occupation
matrix can arise either from the change of the overlap matrix $S$, or from
that of the KS density matrix $\rho$, namely,

$\displaystyle\frac{dn_{I,mm^{\prime}}^{\sigma}[S,\rho]}{d{\bm{\tau}}_{a}}=$
$\displaystyle\frac{\partial n_{I,mm^{\prime}}^{\sigma}[S,\rho]}{\partial
S}\bigg{|}_{\rho}\frac{\partial S}{\partial{\bm{\tau}}_{a}}+$ (27)
$\displaystyle\frac{\partial
n_{I,mm^{\prime}}^{\sigma}[S,\rho]}{\partial\rho}\bigg{|}_{S}\frac{\partial\rho}{\partial{\bm{\tau}}_{a}}\,.$

Thus, the force correction brought by DFT+$U$ is also split into two
contributions. The first contribution, arising from the change of the overlap
matrix, is given by

$\begin{aligned} &\sum_{\sigma}\sum_{I}\sum_{mm^{\prime}}\Delta
V_{I,mm^{\prime}}^{\sigma}\frac{\partial
n_{I,mm^{\prime}}^{\sigma}[S,\rho]}{\partial S}\bigg{|}_{\rho}\frac{\partial
S}{\partial{\bm{\tau}}_{a}}\\\
&=\sum_{\sigma}\frac{1}{4N_{\textbf{k}}}\sum_{\textbf{k}}\sum_{I}\sum_{mm^{\prime}}\Delta
V_{I,mm^{\prime}}^{\sigma}\sum_{\mu\in a}\\\
&\quad\left(\frac{\text{d}S_{\beta m,\mu}\left(\textbf{k}\right)}{\text{d}{\bf
r}_{\beta m,\mu}}\rho_{\mu,\beta
m^{\prime}}^{\sigma}\left(\textbf{k}\right)+\rho_{\mu,\beta
m}^{\sigma*}\left(\textbf{k}\right)\frac{\text{d}S_{\beta
m^{\prime},\mu}^{*}\left(\textbf{k}\right)}{\text{d}{\bf r}_{\beta
m^{\prime},\mu}}\right.\\\ &\quad+\left.\frac{\text{d}S_{\beta
m^{\prime},\mu}\left(\textbf{k}\right)}{\text{d}{\bf r}_{\beta
m^{\prime},\mu}}\rho_{\mu,\beta
m}^{\sigma}\left(\textbf{k}\right)+\rho_{\mu,\beta
m^{\prime}}^{\sigma*}\left(\textbf{k}\right)\frac{\text{d}S_{\beta
m,\mu}^{*}\left(\textbf{k}\right)}{\text{d}{\bf r}_{\beta m,\mu}}\right)\\\
&\quad+\sum_{\sigma}\frac{1}{4N_{\textbf{k}}}\sum_{\textbf{k}}\sum_{I=a}\sum_{mm^{\prime}}\Delta
V_{I,mm^{\prime}}^{\sigma}\sum_{\mu}\\\ &\quad\left(\frac{\text{d}S_{\mu,\beta
m}^{*}\left(\textbf{k}\right)}{\text{d}{\bf r}_{\mu,\beta m}}\rho_{\beta
m^{\prime},\mu}^{\sigma*}\left(\textbf{k}\right)+\rho_{\beta
m,\mu}^{\sigma}\left(\textbf{k}\right)\frac{\text{d}S_{\mu,\beta
m^{\prime}}\left(\textbf{k}\right)}{\text{d}{\bf r}_{\mu,\beta
m^{\prime}}}\right.\\\ &\quad+\left.\frac{\text{d}S_{\mu,\beta
m^{\prime}}^{*}\left(\textbf{k}\right)}{\text{d}{\bf r}_{\mu,\beta
m^{\prime}}}\rho_{\beta m,\mu}^{\sigma*}\left(\textbf{k}\right)+\rho_{\beta
m^{\prime},\mu}^{\sigma}\left(\textbf{k}\right)\frac{\text{d}S_{\mu,\beta
m}\left(\textbf{k}\right)}{\text{d}{\bf r}_{\mu,\beta
m}}\right)\end{aligned}~{},$ (28)

where ${\bf r}_{\nu,\mu}={\bm{\tau}}_{a_{\mu}}-{\bm{\tau}}_{a_{\nu}}$ with
$a_{\nu}$ and $a_{\mu}$ refer to the atoms that the NAO basis functions $\nu$
and $\mu$ are centering on, respectively. In the above derivation, we have
used the Hermiticity of the overlap and density matrix, and the following
relationship for the two-center integrals

$\frac{\partial\langle\phi_{\nu}\mid\phi_{\mu}\rangle}{\partial{\bm{\tau}}_{a}}=\left\\{\begin{aligned}
&\langle\phi_{\nu}\mid\frac{\phi_{\mu}}{d{\bm{\tau}_{a}}}\rangle=\frac{d\langle\phi_{\nu}\mid\phi_{\mu}\rangle}{d{\bf
r}_{\nu,\mu}},\mu\in a,\nu\notin a&\\\
&\langle\phi_{\mu}\mid\frac{d\phi_{\nu}}{d{\bm{\tau}_{a}}}\rangle=\frac{d\langle\phi_{\mu}\mid\phi_{\nu}\rangle}{d{\bf
r}_{\mu,\nu}},\mu\notin a,\nu\in a&\\\
&0,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\mu\in a,\nu\in
a~{}\text{or}~{}\mu\notin a,\nu\notin a&\end{aligned}\right.~{}.$ (29)

The second part of the force consists in the contribution from the change of
KS density matrix. Similarly it is denoted as $\frac{\partial
n_{I,mm^{\prime}}^{\sigma}}{\partial{\bm{\tau}}_{a}}|_{S}$ which means that
overlap matrix is fixed. Its contribution to the total force is given by

$\displaystyle\sum_{\sigma}\sum_{I}\sum_{mm^{\prime}}\Delta
V_{I,mm^{\prime}}^{\sigma}\frac{\partial
n_{I,mm^{\prime}}^{\sigma}[S,\rho]}{\partial\rho}\bigg{|}_{S}\frac{\partial\rho}{\partial{\bm{\tau}}_{a}}$
(30)
$\displaystyle=\sum_{\sigma}\frac{1}{4N_{\textbf{k}}}\sum_{\textbf{k}}\sum_{I}\sum_{mm^{\prime}}\Delta
V_{I,mm^{\prime}}^{\sigma}\sum_{\mu}$ $\displaystyle\quad\left(S_{\beta
m,\mu}\left(\textbf{k}\right)\frac{\text{d}\rho_{\mu,\beta
m^{\prime}}^{\sigma}\left(\textbf{k}\right)}{\text{d}{\bm{\tau}}_{a}}+\frac{\text{d}\rho_{\beta
m,\mu}^{\sigma}\left(\textbf{k}\right)}{\text{d}{\bm{\tau}}_{a}}S_{\mu,\beta
m^{\prime}}\left(\textbf{k}\right)\right.$ $\displaystyle\quad+\left.S_{\beta
m^{\prime},\mu}\left(\textbf{k}\right)\frac{\text{d}\rho_{\mu,\beta
m}^{\sigma}\left(\textbf{k}\right)}{\text{d}{\bm{\tau}}_{a}}+\frac{\text{d}\rho_{\beta
m^{\prime},\mu}^{\sigma}\left(\textbf{k}\right)}{\text{d}{\bm{\tau}}_{a}}S_{\mu,\beta
m}\left(\textbf{k}\right)\right)$
$\displaystyle=\frac{1}{N_{\textbf{k}}}\sum_{\textbf{k}\sigma}\sum_{\mu\nu}\Delta
V_{\text{eff},\nu\mu}^{\sigma}\left({\bf
k}\right)\frac{\text{d}\rho_{\mu\nu}}{\text{d}{\bm{\tau}}_{a}}\left({\bf
k}\right)\,.$

This term is the so-called orthogonality force:
$\sum_{\mu\nu}H_{\nu\mu}\partial\rho_{\mu\nu}/\partial{\bm{\tau}}_{a}$.
Readers who are interested in the detailed derivation of this term are
referred to Refs. [34, 37]. However, this term in fact requires no additional
treatment here since it has already been properly included in orthogonality
force term via the usual DFA calculation procedure.

The stress is defined as the derivative of the total energy with respect to
the strain tensor. An efficient evaluation of the stress tensor is essential
for relaxing the shape and size of the unit cell of periodic systems. The
stress contribution from the DFT+U energy correction can be expressed as

$\begin{aligned} \sigma_{\gamma\eta}=\frac{\partial\Delta
E_{\text{DFT+}U}}{\partial\epsilon_{\gamma\eta}}=\sum_{\sigma}\sum_{I}\sum_{mm^{\prime}}V_{I,mm^{\prime}}^{\sigma}\frac{\partial
n_{I,mm^{\prime}}^{\sigma}}{\partial\epsilon_{\gamma\eta}}\end{aligned},$ (31)

where $\gamma$ and $\eta$ denote the Cartesian coordinate indices, and
$\epsilon$ the strain tensor. As pointed out in Ref. [34], stress calculations
require very little extra effort beyond multiplying the counterpart of the
force by ${\bf r}_{\mu,\nu}^{\eta}$, where ${\bf r}_{\mu,\nu}$ is the vector
connecting the atoms where the NAO basis functions $\phi_{\mu}$ and
$\phi_{\nu}$ are centering on. Therefore, similar to the case of force
calculations (cf. Eq. 27), the contributions tothe stress can also be
decomposed into two parts. The first part $\frac{\partial
n_{I,mm^{\prime}}^{\sigma}}{\partial\epsilon_{\gamma\eta}}|_{\rho}$, arising
from the derivative of the overlap matrix with respect to the strain tensor,
is given by

$\displaystyle\sum_{\sigma}\sum_{I}\sum_{mm^{\prime}}V_{I,mm^{\prime}}^{\sigma}\frac{\partial
n_{I,mm^{\prime}}^{\sigma}}{\partial\epsilon_{\gamma\eta}}|_{\rho}$ (32)
$\displaystyle=\sum_{\sigma}\frac{1}{4N_{\textbf{k}}}\sum_{\textbf{k}}\sum_{I}\sum_{mm^{\prime}}V_{I,mm^{\prime}}^{\sigma}\sum_{\mu}$
$\displaystyle\left(\rho_{\mu,\beta
m^{\prime}}^{\sigma}\left(\textbf{k}\right)\sum_{\textbf{R}}\frac{\text{d}S_{\beta
m,\mu}\left(\textbf{R}\right)}{\text{d}{\bf r}_{\beta m,\mu}^{\gamma}}{\bf
r}_{\beta m,\mu}^{\eta}e^{-i\textbf{k}\textbf{R}}\right.$
$\displaystyle\quad+\rho_{\beta
m,\mu}^{\sigma}\left(\textbf{k}\right)\sum_{\textbf{R}}\frac{\text{d}S_{\mu,\beta
m^{\prime}}\left(\textbf{R}\right)}{\text{d}{\bf r}_{\mu,\beta
m^{\prime}}^{\gamma}}{\bf r}_{\mu,\beta
m^{\prime}}^{\eta}e^{-i\textbf{k}\textbf{R}}$
$\displaystyle\quad+\rho_{\mu,\beta
m}^{\sigma}\left(\textbf{k}\right)\sum_{\textbf{R}}\frac{\text{d}S_{\beta
m^{\prime},\mu}\left(\textbf{R}\right)}{\text{d}{\bf r}_{\beta
m^{\prime},\mu}^{\gamma}}{\bf r}_{\beta
m^{\prime},\mu}^{\eta}e^{-i\textbf{k}\textbf{R}}$
$\displaystyle\quad+\left.\rho_{\beta
m^{\prime},\mu}^{\sigma}\left(\textbf{k}\right)\sum_{\textbf{R}}\frac{\text{d}S_{\mu,\beta
m}\left(\textbf{R}\right)}{\text{d}{\bf r}_{\mu,\beta m}^{\gamma}}{\bf
r}_{\mu,\beta m}^{\eta}e^{-i\textbf{k}\textbf{R}}\right)\,.$

The second part of the stress is the counterpart of the orthogonality force
term (cf. Eq. (30)). For the same reason as in force calculations, this term
has been included automatically in the total orthogonality stress in usual DFA
calculations and requires no additional treatment. Hence, for brevity, its
explicit expression is not given here.

###  II.4 Spin-orbit coupling

In the discussion of the NAO-based DFT+$U$ formalism presented above, the SOC
effect is neglected. Physically the SOC stems from the interaction between the
intrinsic magnetic moment of the electrons and the magnetic field induced by
their orbital angular momenta and is a consequence of the relativistic effect.
The magnitude of the SOC increases with the atomic number; for systems
containing heavy elements, such effect must be taken into account to obtain
physically meaningful results. Below we discuss how SOC is incorporated in our
DFT$+U$ implementation.

When the SOC is present, the KS eigenstates $|\psi_{n\textbf{k}}\rangle$
become two-component spinors,

$|\psi_{n\textbf{k}}\rangle=\sum_{\mu}\left(\begin{array}[]{c}c_{n\mu,\textbf{k}}^{\uparrow}\\\
c_{n\mu,\textbf{k}}^{\downarrow}\end{array}\right)|\phi_{{\bf
k}\mu}\rangle=|\psi_{n\textbf{k}}^{\uparrow}\rangle|\uparrow\rangle+|\psi_{n\textbf{k}}^{\downarrow}\rangle|\downarrow\rangle$
(33)

where —$\uparrow\rangle$ and —$\downarrow\rangle$ are the up- and down-channel
spin states, and $\psi_{n\textbf{k}}^{\uparrow}({\bf r})$ and
$\psi_{n\textbf{k}}^{\downarrow}({\bf r})$ are the associated spatial wave
functions. The two-component spinor can be regarded as a superposition state
of its two components: state $|\Psi_{n\textbf{k}}^{\uparrow}\rangle$
multiplied with spin-up function and state
$|\Psi_{n\textbf{k}}^{\downarrow}\rangle$ with spin-down function. The form of
such two-component eigenstates is different from the one-component eigenstates
$\Psi_{n\textbf{k}}^{\sigma}\left(\textbf{r}\right)$ of the non-SOC
Hamiltonian, whereby the variables of spin $\sigma$ and position r are
independent. In the spirit of the method of separation of variables, the
eigenstates $\Psi_{n\textbf{k}}^{\sigma}\left(\textbf{r}\right)$ is product of
a function of the position $\Psi_{n\textbf{k}}\left(\textbf{r}\right)$ and a
eigenstate of spin variables $\alpha_{n}\left(\sigma\right)$ with eigenvalue
$\sigma$. Obviously the $|\Psi_{n\textbf{k}}\rangle$ in Eq. (33) cannot be
reduced to such a form.

The DFT+$U$ energy functional is given by the expectation value of the second-
quantized Hubbard Hamiltonian, i.e., Eq. (2), in the local subspace within the
Hartree-Fock ground-state, given by the Slater determinant that is formed by
the lowest $N$ spinors. After some simple derivations, we obtain

$\begin{aligned}
E_{U}=&\frac{1}{2}\sum_{I}\sum_{\\{m\\}}\sum_{\sigma\sigma^{\prime}}\left\langle
mm^{\prime}\left|v_{sc}\right|m^{\prime\prime}m^{\prime\prime\prime}\right\rangle
n_{I,mm^{\prime\prime}}^{\sigma\sigma}n_{I,m^{\prime}m^{\prime\prime\prime}}^{\sigma^{\prime}\sigma^{\prime}}\\\
&-\frac{1}{2}\sum_{I}\sum_{\\{m\\}}\sum_{\sigma\sigma^{\prime}}\langle
mm^{\prime}\left|v_{\mathrm{sc}}\right|m^{\prime\prime\prime}m^{\prime\prime}\rangle
n_{I,mm^{\prime\prime}}^{\sigma\sigma^{\prime}}n_{I,m^{\prime}m^{\prime\prime\prime}}^{\sigma^{\prime}\sigma}\end{aligned}~{}.$
(34)

The $n_{I,mm^{\prime}}^{\sigma\sigma^{\prime}}$ is the local occupation matrix
within the SOC scheme that will be addressed later. Following a similar
procedure as Eq. (4)$\sim$ Eq. (6), one can also introduce a unitary
transformation to Eq. (34) and then arrive at

$\begin{aligned}
E_{U}=&\frac{1}{2}\sum_{I}U_{I}\sum_{mm^{\prime}}\sum_{\sigma}\textbf{n}_{I,\sigma
m}\textbf{n}_{I,-\sigma m^{\prime}}\\\
&+\frac{1}{2}\sum_{I}\left(U_{I}-J_{I}\right)\sum_{m\neq
m^{\prime}}\sum_{\sigma}\textbf{n}_{I,\sigma m}\textbf{n}_{I,\sigma
m^{\prime}}\end{aligned}~{}.$ (35)

In this derivation, the $n_{I,mm^{\prime}}^{\sigma\sigma^{\prime}}$ is
regarded as a local occupation matrix in the spin-orbit representation rather
a four-order tensor, i.e., $\sigma m$ is a compact index of spin-orbit. The
$\textbf{n}_{I,\sigma m}$ is the diagonalized local occupation matrix of
$n_{I,mm^{\prime}}^{\sigma\sigma^{\prime}}$. For simplicity, here we have
assumed that the “off-diagonal blocks” of
$n_{I,mm^{\prime}}^{\sigma\sigma^{\prime}}$ with $\sigma\neq\sigma^{\prime}$
are much smaller in magnitude than the “diagonal blocks” with
$\sigma=\sigma^{\prime}$, and that the eigenvalues obtained by diagonalizing
$n_{I,mm^{\prime}}^{\sigma\sigma}$ separately for $\sigma=\uparrow,\downarrow$
don’t differ appreciably from those obtained by diagonalizing the full local
occupation matrix. Test calculations indicate that this is a rather good
approximation. Subtracting the same double counting term as Eq. (7) we get the
energy correction within SOC scheme as

$\displaystyle\Delta E_{\text{DFT+}U}=$
$\displaystyle\frac{1}{2}\sum_{I}\left(U_{I}-J_{I}\right)\sum_{m\sigma}\left(\textbf{n}_{I,\sigma
m}-\textbf{n}_{I,\sigma m}\textbf{n}_{I,\sigma m}\right)$ (36)
$\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{I}\bar{U}_{I}\sum_{\sigma}\left[\sum_{m}n_{I,mm}^{\sigma\sigma}\right.$
$\displaystyle\left.-\sum_{mm^{\prime},\sigma^{\prime}}n_{I,mm^{\prime}}^{\sigma\sigma^{\prime}}n_{I,m^{\prime}m}^{\sigma^{\prime}\sigma}\right].$

The effective single-particle potential is given by

$V_{I,mm^{\prime}}^{\sigma\sigma^{\prime}}=\bar{U}_{I}\left(1/2\delta_{mm^{\prime}}\delta_{\sigma\sigma^{\prime}}-n_{I,mm^{\prime}}^{\sigma\sigma^{\prime}}\right).$
(37)

Particularly, in this case, the Mulliken charge projector in Eq. (15) is
generalized to a $2\times 2$ tensor in the spin space,

$\begin{aligned}
\hat{P}_{I,mm^{\prime}}^{\sigma,\sigma^{\prime}}=&\frac{1}{4N_{\textbf{k}}}\sum_{\textbf{k}}\left(|\tilde{\phi}_{\textbf{k},\beta
m^{\prime}}\sigma^{\prime}\rangle\langle\phi_{\textbf{k},\beta
m}\sigma|\right.\\\ &+|\phi_{\textbf{k},\beta
m^{\prime}}\sigma^{\prime}\rangle\langle\tilde{\phi}_{\textbf{k},\beta
m}\sigma|+|\tilde{\phi}_{\textbf{k},\beta
m}\sigma\rangle\langle\phi_{\textbf{k},\beta m^{\prime}}\sigma^{\prime}|\\\
&\left.+|\phi_{\textbf{k},\beta
m}\sigma\rangle\langle\tilde{\phi}_{\textbf{k},\beta
m^{\prime}}\sigma^{\prime}|\right)\end{aligned}~{}.$ (38)

Thus the local occupation matrix in SOC case becomes

$\begin{aligned}
n_{I,mm^{\prime}}^{\sigma\sigma^{\prime}}=&\frac{1}{N_{\textbf{k}}}\sum_{n\textbf{k}}f_{n\textbf{k}}\langle\psi_{n\textbf{k}}|\hat{P}_{I,mm^{\prime}}^{\sigma,\sigma^{\prime}}|\psi_{n\textbf{k}}\rangle\\\
=&\frac{1}{4N_{\textbf{k}}}\sum_{\textbf{k}}\sum_{\mu}\left(S_{\beta
m,\mu}\left(\textbf{k}\right)\rho_{\mu,\beta
m^{\prime}}^{\sigma\sigma^{\prime}}\left(\textbf{k}\right)\right.\\\
&+\left.\rho_{\beta
m,\mu}^{\sigma\sigma^{\prime}}\left(\textbf{k}\right)S_{\mu,\beta
m^{\prime}}\left(\textbf{k}\right)+S_{\beta
m^{\prime},\mu}\left(\textbf{k}\right)\rho_{\mu,\beta
m}^{\sigma^{\prime}\sigma}\left(\textbf{k}\right)\right.\\\
&+\left.\rho_{\beta
m^{\prime},\mu}^{\sigma^{\prime}\sigma}\left(\textbf{k}\right)S_{\mu,\beta
m}\left(\textbf{k}\right)\right)\end{aligned}~{}.$ (39)

Within the scheme of two-component spinors, the density matrix becomes a
$2\times 2$ tensor in the spin space

$\rho\left(\textbf{k}\right)=\left[\begin{array}[]{c
c}\rho_{\mu\nu}^{\uparrow\uparrow}\left(\textbf{k}\right)&\rho_{\mu\nu}^{\uparrow\downarrow}\left(\textbf{k}\right)\\\
\rho_{\mu\nu}^{\downarrow\uparrow}\left(\textbf{k}\right)&\rho_{\mu\nu}^{\downarrow\downarrow}\left(\textbf{k}\right)\end{array}\right],$
(40)

where

$\rho_{\mu\nu}^{\sigma\sigma^{\prime}}\left(\textbf{k}\right)=\sum_{n}f_{n\textbf{k}}c_{n\textbf{k},\mu}^{\sigma}c_{n\textbf{k},\nu}^{\sigma^{\prime}*}~{}.$
(41)

Then the effective single-particle potential operator becomes

$\hat{V}_{\text{eff}}^{\sigma\sigma^{\prime}}\left(\textbf{k}\right)=\sum_{I}\sum_{mm^{\prime}}V_{I,mm^{\prime}}^{\sigma\sigma^{\prime}}\hat{P}_{I,mm^{\prime}}^{\sigma\sigma^{\prime}}\left(\textbf{k}\right),$
(42)

where

$\begin{aligned}
\hat{P}_{I,mm^{\prime}}^{\sigma,\sigma^{\prime}}\left(\textbf{k}\right)=&\frac{1}{4}\left(|\tilde{\phi}_{\textbf{k},\beta
m^{\prime}}\sigma^{\prime}\rangle\langle\phi_{\textbf{k},\beta
m}\sigma|+|\phi_{\textbf{k},\beta
m^{\prime}}\sigma^{\prime}\rangle\langle\tilde{\phi}_{\textbf{k},\beta
m}\sigma|\right.\\\ &\left.+|\tilde{\phi}_{\textbf{k},\beta
m}\sigma\rangle\langle\phi_{\textbf{k},\beta
m^{\prime}}\sigma^{\prime}|+|\phi_{\textbf{k},\beta
m}\sigma\rangle\langle\tilde{\phi}_{\textbf{k},\beta
m^{\prime}}\sigma^{\prime}|\right)\end{aligned}~{}.$ (43)

###  II.5 $U$ and $J$ parameters determined from the Yukawa potential

As shown in the previous sections, the Hubbard U and Hund J are two key
parameters in DFT+U calculations. They can be taken as empirical parameters,
or determined from pragmatic schemes like constrained DFT [49, 50, 51, 23],
constrained RPA [52, 53, 54, 55], or linear-response approach [31]. In
interacting many-electron systems, the Coulomb interaction between electrons
is screened, resulting weaker and often shorter-ranged effective interactions.
In principle, one could directly model such a screened Coulomb potential, and
directly use Eq. (3) to compute the $U$ and $J$ parameters. Previously, a
simple form of such a screened potential – the Yukawa potential has been
employed in the literature, and demonstrated to work reasonably well [58, 59,
60]. In the present work, we follow such an approach and check how it works
within our NAO-based DFT+$U$ scheme.

The Yukawa potential reads

$v_{\text{sc}}\left(\mathbf{r},\mathbf{r}^{\prime}\right)=\frac{e^{-\lambda|\mathbf{r}-\textbf{r}^{\prime}|}}{|\mathbf{r}-\mathbf{r}^{\prime}|}$
(44)

where $\lambda$ is a screening parameter. Compared to the bare Coulomb
potential, the Yukawa potential decays exponentially fast to zero for large
separations of two spatial points. Mathematically screened Coulomb interaction
matrix elements based on the Yukawa potential given by Eq. (3) can be
decomposed into two parts, i.e., an angular integral part involving spherical
harmonics and a radial integral part called the Slater integrals [61], namely,

$\langle
mm^{\prime}\left|v_{sc}\right|m^{\prime\prime}m^{\prime\prime\prime}\rangle=\sum_{k=0}^{2l}a_{k}\left(m,m^{\prime},m^{\prime\prime}m^{\prime\prime\prime}\right)F^{\left(k\right)}\,.$
(45)

In Eq. (45),

$\displaystyle
a_{k}\left(m,m^{\prime},m^{\prime\prime}m^{\prime\prime\prime}\right)=$ (46)
$\displaystyle\frac{4\pi}{2k+1}\sum_{q=-k}^{k}\int
d\hat{\mathbf{r}}Y_{lm}^{*}\left(\hat{\mathbf{r}}\right)Y_{kq}\left(\hat{\mathbf{r}}\right)Y_{lm^{\prime\prime}}\left(\hat{\mathbf{r}}\right)$
$\displaystyle\times\int
d\hat{\mathbf{r}}^{\prime}Y_{lm^{\prime}}^{*}\left(\hat{\mathbf{r}}^{\prime}\right)Y_{kq}^{*}\left(\hat{\mathbf{r}}^{\prime}\right)Y_{lm^{\prime\prime\prime}}\left(\hat{\mathbf{r}}^{\prime}\right)$

is the angular part which can be easily evaluated by Gaunt coefficients, and

$F^{\left(k\right)}=\int\int
R^{2}_{l}\left(r\right)v_{sc}^{\left(k\right)}\left(r,r^{\prime}\right)R_{l}^{2}\left(r^{\prime}\right)r^{2}r^{\prime
2}drdr^{\prime}$ (47)

is the radial part. The $v_{sc}^{\left(k\right)}\left(r,r^{\prime}\right)$ is
the $k$-order coefficient of the expansion of
$v_{sc}\left(\mathbf{r},\mathbf{r}^{\prime}\right)$ by spherical harmonics
[60], i.e.,

$v_{sc}\left(\mathbf{r},\mathbf{r}^{\prime}\right)=\sum_{k=0}^{\infty}\frac{4\pi}{2k+1}v_{sc}^{\left(k\right)}\left(r,r^{\prime}\right)\sum_{q=-k}^{k}Y_{kq}^{*}\left(\hat{\mathbf{r}}\right)Y_{kq}\left(\hat{\mathbf{r}}^{\prime}\right)~{}.$
(48)

In case of the Yukawa potential form, $F^{\left(k\right)}$ in Eq. (47) is
further reduced to

$\displaystyle F^{\left(k\right)}=-\left(2k+1\right)\lambda$ (49)
$\displaystyle\times\int\int R^{2}_{l}\left(r\right)j_{k}\left(i\lambda
r_{<}\right)h_{k}^{(1)}\left(i\lambda
r_{>}\right)R_{l}^{2}\left(r^{\prime}\right)r^{2}r^{\prime 2}drdr^{\prime},$

where $j_{k}$ and $h_{k}^{(1)}$ are the spherical Bessel function and the
spherical Hankel function of the first kind at order $k$, respectively.
$r_{>}$ and $r_{<}$ are the smaller and the larger radius entering in the
integral, i.e. $r_{>}=\text{max}\left(r,r^{\prime}\right)$ and
$r_{<}=\text{min}\left(r,r^{\prime}\right)$. In this formulation, for a given
set of local orbitals, the matrix elements of the screened Coulomb potential
as given by Eq. (45), from which the $U$ and $J$ parameters can be extracted,
depends only on the screening parameter $\lambda$. Thus, the problem of
determining the $U$, $J$ values becomes one for determining $\lambda$.

In the theory of Thomas-Fermi model, the screening parameter $\lambda$ is a
function of the charge density of the system,

$\lambda=2\left[\frac{3\rho}{\pi}\right]^{1/6}~{},$ (50)

where $\rho$ is the electron density. In this work, we use the effective
screening parameter $\bar{\lambda}$ which is obtained by averaging the space-
dependent screening parameter $\lambda({\bf r})$

$\bar{\lambda}=\frac{\int\lambda\left({\bf r}\right)\rho\left({\bf
r}\right)d{\bf r}}{\int\rho\left({\bf r}\right)d{\bf r}}~{},$ (51)

where the $\rho\left({\bf r}\right)$ is selected as the pseudo charge density
in the KS self-consistent iteration. This model provides us with a possible
scheme of parameter-free DFT+U with U and J determined in a self-consistent
way.

For practical DFT+U calculations, determining the $U$ and $J$ values based on
the Yukawa potential brings simplification. Firstly, in standard DFT+U
calculations, two parameters, i.e. Hubbard U and Hund J, need to be
determined, while in the Yukawa-potential approach only one parameter
$\lambda$ is needed. The screening parameter $\lambda$ can be evaluated self-
consistently with the help of the Thomas-Fermi screening model. Thus a
parameter-free DFT+U scheme is in principle achievable. Secondly, for systems
that have not been well investigated and no reference results are available,
it is highly nontrivial to obtain the appropriate U and J values. Under such
circumstances, a universal screening parameter $\lambda$ in this approach can
provide initial information for further investigation. Thirdly, for accurate
NAO-based calculations, there are often more than one radial function used for
each $d/f$ angular momentum channel. This raises the question if the Hubbard
$U$ correction need to be applied to all these $d/f$ function channels and how
to determine the $U$, $J$ parameters for each individual orbital. Based on Eq.
(45), the Yukawa potential approach allows one to conveniently determine the
orbital-dependent $U$, $J$ parameters, from which one can readily decide the
relevant orbital channels where the Hubbard $U$ correction is necessary.

##  III Computational details

Our DFT+$U$ implementation is carried out within the ABACUS code package [37].
In the present work, we use the SG15 optimized norm-conserving Vanderbilt
(ONCV) multi-projector pseudo-potentials [62, 63, 64] to describe the ion
cores, and optimized double-$\zeta$ plus polarization (DZP) atomic basis sets
[65, 41] to expand the Kohn-Sham eigenfunctions. Within DZP, $4s2p2d1f$ atomic
functions with a cutoff radius of 9 Bohr are used for the TM atoms and
$2s2p1d$ atomic functions with $r_{\text{cutoff}}$ of
$7\text{\,}\mathrm{B}\mathrm{o}\mathrm{h}\mathrm{r}$ for oxygen (O) atoms. In
certain occasions, the triple-$\zeta$ plus polarization (TZDP) basis set with
$5s3p3d2f$ for TMs and $3s3p2d$ for O are also used. For Brillouin zone
sampling, a $\Gamma$-inclusive $6\times 6\times 6$ k-mesh is used. In band
structure calculations, we set up k-point paths explicitly along specified
high-symmetry directions of the Brillouin zone. We consider ten high-symmetry
points in k-space for the rhombohedral Bravais lattice of type-II
antiferromagnetic (AFM) structures (e.g. late TM monoxides MnO, FeO, CoO and
NiO). The corresponding k-path is $\Gamma-L-B_{1}|B-Z-\Gamma-X|Q-F-
P_{1}-Z|L-P$ [66]. As for the case of simple tetragonal IrO2, the ${\bf
k}$-path is $\Gamma-X-M-\Gamma-Z-R-A-M$ [67]. For all calculations, the
Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation (GGA) is used
as the exchange-correlation functional [68] at the level of DFA, i.e., the
PBE+$U$ scheme is employed in our test calculations below.

To validate our scheme and implementation we make comparisons to several
widely used softwares, including the LAPW method based code WIEN2k, Vienna ab
initio simulation package (VASP) which is based on the PAW method, the
Quantum-ESPRESSO (QE) package which uses pseudopotential and plane-wave basis
set [30, 31], and LCAO based OpenMX code [69, 70, 71, 45]. For all
calculations, we use the isotropic version of the simplified DFT+U method, in
which only the difference between Hubbard U and Hund J matters [57]. For all
WIEN2k based calculations, we set the muffin-tin radii of TM atoms as
$2.1\text{\,}\mathrm{B}\mathrm{o}\mathrm{h}\mathrm{r}$ and that for O atoms as
$1.4\text{\,}\mathrm{B}\mathrm{o}\mathrm{h}\mathrm{r}$. The convergence of the
calculations is controlled by $R_{\text{MT}}\times K_{\text{max}}=7.0$. In
VASP calculations, we set the energy cutoff as $700\text{\,}\mathrm{eV}$, and
non-spherical contribution in gradient corrections inside the PAW sphere is
included (LASPH=.True.). In order to make comparisons, we use the same SG15
ONCV multi-projector pseudopotentials in QE calculations. The kinetic energy
cutoffs for wavefunctions and charge density are 60 and
$480\text{\,}\mathrm{R}\mathrm{y}$ respectively. For OpenMX results, we simply
cite those reported in Ref. [45].

The introduction of multiple radial functions (multiple-$\zeta$) with the same
angular momentum to the basis set in the LCAO framework makes it intricate to
define the projector, and hence the local subspace. The definition of the
projector in section II.2 indicates that the correlated sub-shell is not
uniquely specified with only the correlated site $I$ and angular momentum $l$,
because there is still a $\zeta$ index. In the present work, the DZP basis
sets we are using contain two $d$-type functions and there is naturally the
question regarding how to define the $\beta$ index when constructing the
projector (cf. Eq. (15) and (34)). In our current procedure [65, 41] to
optimize the NAO basis sets, the first generated $d$ function of TMs is most
localized and represents best the local correlated subspace around the Fermi
level. Also, the estimated $U$ and $J$ parameters for these orbitals based on
the Yukawa potential as described in Sec. II.5 also agree fairly well with
those reported in the literature. In contrast, the second $d$ function is
rather delocalized and contributes little to the correlated subspace. As such,
in the calculations reported below, the first $d$ orbital is used to construct
the projector. Further discussions about our choice and its efficacy will be
discussed in section IV.5 and appendix A.

##  IV Results and discussion

###  IV.1 Band gaps and magnetic moment

One of the most noteworthy successes of the DFT+U approach is the reproducing
of the observed insulating behavior of TMOs. In the standard DFT scheme, LDA
and GGA either predict conducting behavior or yield gaps that are much smaller
than the experimental values. In TMOs, correlated valence d electrons tend to
be localized due to strong on-site Coulomb interaction. Local and semi-local
DFAs fail to describe this localization behavior. This deficiency is largely
rectified by DFT+$U$. Here, we apply our NAO-based DFT+$U$ implementation to
prototypical Mott insulators MnO, FeO, CoO and NiO, and check how the obtained
band gaps vary with the $\bar{U}$ value. We also compare our results with
those produced by established computer codes, whereby the validity of our
implementation can be demonstrated.

In our calculations, the NaCl-type structure is assumed for the four late TMOs
– MnO, FeO, CoO and NiO – where the chosen lattice constants are
$4.445\text{\,}\mathrm{\SIUnitSymbolAngstrom}$,
$4.334\text{\,}\mathrm{\SIUnitSymbolAngstrom}$,
$4.254\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ and
$4.171\text{\,}\mathrm{\SIUnitSymbolAngstrom}$, respectively [72]. For all
these four systems we consider the type-II AFM structure, where the AFM order
is along the $\left(111\right)$ direction [31, 73, 60]. The SOC effect is not
considered for these systems.

Table 1 presents the calculated gaps as a function of the effective on-site
Coulomb interaction parameter $\bar{U}$. Without the $U$ correction, PBE
yields zero gaps for FeO and CoO, and significantly underestimates the band
gaps for MnO and NiO. When the $U$ correction is applied, the obtained band
gap opens up and steadily increases with an increasing $\bar{U}$ value, as
expected. When $\bar{U}$ reaches a physically meaningful value of about 6 eV,
our calculated PBE+$U$ band gaps show a fairly good agreement with the
experimental values for all four TMOs.

Table 1:  Band gaps (in eV) of MnO, FeO, CoO and NiO as a function of
effective on-site Coulomb energy $\bar{U}$ (in eV). The experimental values
are presented in the last row.

$\bar{U}$ (eV) |  MnO |  FeO |  CoO |  NiO
---|---|---|---|---
0.0 |  1.13 |  0.00 |  0.00 |  0.82
2.0 |  1.65 |  0.76 |  1.84 |  1.88
4.0 |  1.98 |  2.44 |  2.57 |  2.60
6.0 |  2.18 |  2.62 |  3.07 |  3.21
Expt. |  3.6-3.8${}^{\text{a}}$ |  2.4${}^{\text{b}}$ |  2.4${}^{\text{c}}$ |  4.0${}^{\text{d}}$, 4.3${}^{\text{e}}$

  * 1

${}^{\text{a}}$Reference [74]

  * 2

${}^{\text{b}}$Reference [75]

  * 3

${}^{\text{c}}$Reference [76]

  * 4

${}^{\text{d}}$Reference [77]

  * 5

${}^{\text{e}}$Reference [78]

Since different projectors and numerical frameworks are used in different
computer codes, the calculated DFT+$U$ band gaps could differ noticeably even
with the same $U$ value. However, the results should be qualitatively similar.
For instance, the calculated DFT+$U$ band gaps should follow a similar trend
as the $U$ value increases. To check this, in Fig. 1 we present the band gaps
for the four TMOs as a function of $\bar{U}$, as obtained by five computer
codes: ABACUS, OpenMX, VASP, WIEN2k, and QE. The OpenMX results for all the
four systems are taken from Ref. [45], where LDA functional was used. For all
other codes, the PBE functional is used in DFT+$U$ calculations. We don’t
expect using PBE instead of LDA will produce significant difference for the
curves presented in Fig. 1. For MnO, the results of VASP, WIEN2k and QE are
taken from Ref. [73], whereas the results for other TMOs are calculated in the
present work. Figure 1 indicates the band gaps obtained from different codes
follow the same trend as $\bar{U}$ increases, but the obtained band gaps can
differ by as large as 1 eV for the same $\bar{U}$ value. This means that the
“best” $\bar{U}$ value to reproduce the experimental results varies for
different numerical frameworks. Compared to other codes, ABACUS displays a
rather “normal” behavior, which is a strong evidence of the validity of our
projector scheme and numerical implementations.

Figure 1: Band gaps of MnO, FeO, CoO and NiO as a function of $\bar{U}$ as
calculated by different computer codes. The rhombohedral unit cell with type-
II AFM structure is used. For FeO, we couldn’t obtain stable self-consistent
results for all $\bar{U}$ values with QE, and hence the QE results are not
included in FeO panel. The LDA+$U$ results taken from Ref. ([45]) are
presented for OpenMX, whereas for all other codes the PBE+$U$ results are
presented.

The local magnetic moment of TM atoms within MnO, FeO, CoO and NiO are also
calculated by ABACUS and the results are presented in Table 2. Table 2 shows
the expected trend that the local magnetic moments get enhanced with
increasing $\bar{U}$. With physical $\bar{U}$ values, the obtained PBE+$U$
magnetic moments show an overall good agreement with experimental results and
previously reported theoretical results [23, 45].

Table 2:  Magnetic moments (in $\mu_{B}$) of TM atoms within the four TMOs as
a function of effective on-site Coulomb energy $\bar{U}$.

$\bar{U}$ (eV) |  MnO |  FeO |  CoO |  NiO
---|---|---|---|---
0.0 |  4.39 |  3.44 |  2.47 |  1.21
2.0 |  4.64 |  3.60 |  2.56 |  1.53
4.0 |  4.74 |  3.70 |  2.69 |  1.64
6.0 |  4.80 |  3.77 |  2.74 |  1.71
Expt. |  4.58${}^{\text{a}}$, 4.79${}^{\text{b}}$ |  3.32${}^{\text{c}}$ |  3.8${}^{\text{c}}$, 3.35${}^{\text{d}}$ |  1.90${}^{\text{a}}$, 1.77${}^{\text{b}}$
Theo. |  4.61${}^{\text{e}}$ |  3.62${}^{\text{e}}$ |  2.63${}^{\text{e}}$ |  1.69${}^{\text{e}}$, 1.74${}^{\text{f}}$

  * 1

${}^{\text{a}}$Reference [79]

  * 2

${}^{\text{b}}$Reference [80]

  * 3

${}^{\text{c}}$Reference [81]

  * 4

${}^{\text{d}}$Reference [82]

  * 5

${}^{\text{e}}$Reference [23]

  * 6

${}^{\text{f}}$Reference [45] ($U$=$6.0\text{\,}\mathrm{eV}$)

###  IV.2 Projected density of states and band structures

To gain more insights about the performance of our DFT+U scheme, we calculate
the projected density of states (PDOS) and band structures of MnO, FeO, CoO
and NiO at a fixed effective on-site Coulomb interaction parameter $\bar{U}$
of $5.0\text{\,}\mathrm{eV}$, and the obtained results are plotted in Fig. 2
and Fig. 3, respectively . In the type–II AFM unit cells, the spin
polarizations of the two neighboring TM atoms are of the same magnitude but
anti-parallel in direction. In Fig. 2, we project the total density of states
(TDOS) to the spin-up polarized TM atoms and their nearest oxygen atoms.
Figure 2(a) indicates that 3d electrons of the projecting site Mn are nearly
fully spin-up polarized, which suggests all 3d electrons occupy the spin up
states and this result agrees with the previous work using WIEN2k, VASP and QE
codes [73]. As the number of $d$ electrons increases from Mn to Ni, the spin-
down $d$ states also get populated, resulting in a decrease of the total
magnetic moments (cf. Table 2). In all four mono-oxides, there are significant
hybridizations between the O 2p and TM 3d orbitals over a wide range of
energy. In particular, the top valence states have a predominant contribution
from the O 2$p$ orbitals, suggesting a charge transfer character of these
TMOs. Such behaviors agree with the chemistry bond theory and previous
theoretical results [23, 31, 45, 73].

Figure 2: The TDOS and PDOS of spin-up polarized TM atom and its nearest O
atom of MnO, FeO, CoO and NiO calculated with $\bar{U}$ of
$5\text{\,}\mathrm{eV}$ are depicted in (a), (b), (c), and (d), respectively.
The zeroes of the energy axes are set to be the Fermi level (vertical dashed
line). The upper and lower panels correspond to the up-spin and down-spin
channels.

We further present the calculated band structures of the four TMOs in Fig. 3,
which display the typical character of strongly correlated systems. The
occupied manifold of the KS states is a mixture of the O $2p$ and TM $3d$
characters. Energetically the TM 3$d$ dominating bands sit below the O 2$p$
dominating ones, due to the fact that the strong Coulomb repulsion pushes the
occupied TM 3$d$ states deeply down in energy. Furthermore, when going from
MnO to NiO, a part of the low conduction bands transfers into the valence
bands, while maintaining an insulating band gap of about 3 eV (cf. Fig.
3(a)-(d)). Both the top valence bands and bottom conduction bands show a very
small dispersion of the order of $1.0\text{\,}\mathrm{eV}$, which is much
smaller when compared to the effective on-site Coulomb interaction energy
$\bar{U}$. Theoretically, the hoping amplitude of electrons between
neighboring lattice sites is proportional to the bandwidth. In the cases where
the bandwidth is much smaller than the on-site Coulomb interactions, the
transport process of the valence electrons gets hampered and it is very
difficult for them to hop between neighboring sites, so than 3d electrons
become localized, leading to the insulating behavior. The narrow widths of the
top valence bands agree well with this physics picture of strongly correlated
systems. All these properties suggest that our DFT+U implementation
successfully captures the fundamental physical mechanism of strongly
correlated TM monoxides.

Figure 3: The ABACUS PBE+$U$ band structures of MnO, FeO, CoO, and NiO
calculated with $\bar{U}$ of $5.0\text{\,}\mathrm{eV}$ are depicted in panels
(a), (b), (c), and (d), respectively. The energy zero is set to be the top of
the occupied bands (horizontal dashed line).

As a comparison to another implementation scheme, Fig. 4 presents the PBE+$U$
band structures of NiO calculated by ABACUS and the plane-wave based code QE
at $\bar{U}$ of $5.0\text{\,}\mathrm{eV}$. The DFT+U bands of NiO given by the
two codes are fairly close, with the occupied $d$ bands of ABACUS slightly
lower in energy compared to the QE results. This difference may be caused by
the different underlying DFT+U schemes of the two codes. It should be noted
that the same $\bar{U}$ value can lead to somewhat different effects within
different DFT+U schemes.

Figure 4: The PBE+$U$ band structures of NiO calculated by ABACUS and QE with
$\bar{U}=$$5\text{\,}\mathrm{eV}$. The zeroes of energy is set to be the top
of highest occupied bands .

###  IV.3 Validation of the force and stress calculations

To assess the validity of our DFT+$U$ force and stress implementations, we
compare the results calculated via the analytical derivative formalism in Eq.
(26) and (31) with those obtained by the finite-difference (FD) method, taking
NiO in the conventional cell as the testing system. We first calculate the
atomic forces for a series of structures with the Ni atom moving along the
$z$-axis while other atoms staying at their original positions, and the
results are reported in Table 3. As clearly shown in Table 3, the forces based
on analytical gradients and those determined by the FD method agree fairly
well, with the remaining discrepancy below $2\%$.

Table 3: Benchmark force calculations for NiO. The 1st column is the $z$ coordinate of the Ni atom, and the 2nd column shows the energy increase for each move of the Ni atom. The 3rd and 4th columns present the FD and analytical forces (in eV/Å) experienced by the Ni atom, respectively. The last column presents the relative deviations (in percentage) of the results computed by the two approaches.  $z$ (Å) |  $\Delta E$ (eV) |  $\Delta E/\Delta z$ |  Force |  Deviation
---|---|---|---|---
2.297820 |  |  |  -2.563364 |
2.301998 |  0.010869 |  -2.601461 |  -2.637573 |  1.37%
2.306175 |  0.011193 |  -2.679244 |  -2.712991 |  1.24%
2.310353 |  0.011472 |  -2.745873 |  -2.789429 |  1.56%
2.314531 |  0.011831 |  -2.831828 |  -2.868811 |  1.29%

For benchmark stress calculations, we fix the lattice constants $a$ and $b$
and vary the length of $c$ of the NiO (conventional cubic) cell, and the
results are presented in Table 4. Comparison of the FD and analytical stress
results shows that the relative deviations are in the order of 2-3$\%$. Such
level of accuracy is adequate for relaxing the lattice structures.

Table 4: Benchmark stress calculations for NiO. The 1st column is the length of the lattice constant $c$ of the conventional unit cell and the 2nd column shows the energy change of each increase of $c$. The 3rd and 4th columns present the FD and analytical stress results (in kbar), respectively. In the 3rd column, the FD stress is calculated by $\Delta E/{\epsilon\Omega}$ where $\epsilon$ is the strain component of c direction, i.e., $\epsilon_{cc}$ and $\Omega$ is the cell volume. $c$ (Bohr) |  $\Delta E$ (eV) |  $\Delta E/{\epsilon\Omega}$ |  Stress |  Deviations
---|---|---|---|---
7.895000 |  |  |  -47.484882 |
7.902895 |  0.002266 |  -52.926873 |  -51.879179 |  -2.02%
7.910790 |  0.002464 |  -57.497658 |  -56.164020 |  -2.37%
7.918685 |  0.002666 |  -62.141370 |  -60.327858 |  -3.01%
7.926580 |  0.002842 |  -66.188954 |  -64.360639 |  -2.84%

###  IV.4 DFT+$U$+SOC band structure

The strength of SOC scales as $O(Z^{4})$, where $Z$ is the atomic number. The
atomic number of Iridium (Ir) is 77, which is nearly three times larger than
that of Mn (atomic number 25). Thus theoretically the SOC effect of Ir-based
compounds is much stronger than 3d TM materials. In this subsection, we
examine the band structure of IrO2 to assess the performance of our DFT+U+SOC
implementation. We consider the AFM structure with the magnetic momenta of two
neighboring Ir atoms in the conventional cell antiparallel along the $z$-axis.
The DFT+U correction is applied to the 5d electrons of Ir atoms with
$\bar{U}=$2.0\text{\,}\mathrm{eV}$$. The calculated PBE+$U$ and PBE+$U$+SOC
band structures are presented in Fig. 5(a) and Fig. 5(b), respectively.

Figure 5: Band structures of IrO2 calculated by (a) ABACUS DFT+U, (b)ABACUS
DFT+U+SOC, and (c) VASP DFT+U+SOC schemes.

The PDOS analysis (not shown) suggests that the main components of the bands
located in the energy window included in Fig. 5 are 5d electrons of the Ir
atom. Comparing Fig. 5(a) and 5(b), the originally degenerate bands along the
${\bf k}$-path from Z to A are split into subbands after the introduction of
SOC, which is the typical effects of SOC on heavy elements.

For comparison, in Fig. 5(c) the PBE+$U$+SOC band structure calculated by VASP
are also presented. As can be seen, ABACUS and VASP give similar band
structures for IrO2. The small difference in details between ABACUS (5(b)) and
VASP (5(c)) results may stem from different basis sets, different
pseudopotentials and the different DFT+U and SOC schemes. In the latter case,
in VASP the SOC is included in a non-self-consistent way based on a
perturbation approach, whereas in ABACUS, the SOC is treated in a self-
consistent way. In addition, the results in Fig. 5 are in good agreement with
previous theoretical work using all-electron full-potential DFT+U+SOC and
LDA+DMFT method [67].

###  IV.5 Yukawa-potential based $U$ and $J$ parameters

In the above discussions, we focused on the validation of our DFT+$U$
implementation whereby $U$ and $J$ are treated as given parameters. To apply
DFT+$U$ to novel materials where no reference results are available, it is
crucial to be able to determine these parameters from first principles
calculations. As mentioned above, in the literature, different theoretical
schemes have been developed to determine these parameters [51, 23, 52, 53, 54,
55, 31]. In this subsection, we check how the scheme based on the Yukawa
potential, as described in Sec. II.5, works for determining $U$, $J$
parameters.

As discussed in Sec. II.5, for a given set of NAO basis functions, the
calculated interaction parameters based on Eq. (45) depend solely on the
screening parameter $\lambda$. Here, we first examine how the Coulomb
interaction parameter $U$ and the exchange parameter $J$ change with the
$\lambda$ value for the local $d$/$f$-type orbitals We take three TMOs, i.e.,
MnO, FeO, and NiO as test examples, and use NAO DZP basis sets in the
calculations. For the TMs Mn, Fe, and Ni, the DZP basis set contains two $d$
functions and one $f$ function. The Yukawa potential scheme allows one to
readily compute the $U$, $J$ parameters for all these orbital functions.

Table 5 presents the calculated U and J parameters for the two $d$ and one $f$
orbitals with $\lambda$ varying from 0.80 to
$1.08\text{\,}\mathrm{B}\mathrm{o}\mathrm{h}\mathrm{r}^{-1}$. First, one can
see that the $U$, $J$ values for the first $d$ orbital are one order of
magnitude larger than those for the second one. This is because the first $d$
orbital is localized and has no node, whereas the second orbital, designed to
be orthogonalized to the first one, is much more extended and has a node (see
the Appendix for further details). Thus, it is not surprising that the $U$,
$J$ corresponding to the second $d$ orbital are much smaller than those of the
first one. According to Table 5, one can see that the second $d$ orbital can
be safely excluded from the correlated subspace in the DFT+$U$ calculations.
Furthermore, for these TM compounds, the $f$ functions are “polarized
orbitals” and represent orbital space that is well separated from the
correlated subspace around the Fermi level, and consequently require no
special treatment. Consistent with this, the $U$, $J$ parameters for the $f$
orbital are also significantly smaller than those of the first $d$ orbital.

Table 5 shows that the obtained $U$, $J$ values steadily decrease upon
increasing the $\lambda$ value, which is easy to understand since a larger
$\lambda$ means a stronger screening of the Coulomb interaction. Focusing on
the first $d$ orbital ($d1$ in Table 5) and comparing to the results given by
constrained DFT [83], constrained RPA [55] and local screened Coulomb
correction (LSCC) scheme [60], we find that a screening parameter $\lambda$
around $1.00\text{\,}\mathrm{B}\mathrm{o}\mathrm{h}\mathrm{r}^{-1}$ yields
$U$, $J$ values that are fairly close to those reported in the literature.
Remarkably, it seems that
$\lambda\approx$1\text{\,}\mathrm{B}\mathrm{o}\mathrm{h}\mathrm{r}^{-1}$$ is a
sensible choice that works well for all three TMOs. Therefore, at least the TM
compounds, our DFT+$U$ scheme together with the Yukawa-potential scheme for
determining $U$, $J$ parameters can be viewed as a semi-empirical approach
that has predictive power; namely, the only input parameter is fixed and does
not vary for different materials.

Table 5: On-site Coulomb energy $U$ and exchange energy $J$ (in eV) as a
function of screening parameter $\lambda$ (in Bohr-1) for MnO, FeO and NiO in
DZP basis.Here d and f are angular momentum indexes while the number behind
them are multiplicity indexes $\zeta$.

$\lambda$ | MnO |  | FeO |  | NiO
---|---|---|---|---|---
d1 |  |  d2 |  |  f1 |  |  d1 |  |  d2 |  |  f1 |  |  d1 |  |  d2 |  |  f1
$U$ | $J$ |  | $U$ | $J$ |  | $U$ | $J$ |  | $U$ | $J$ |  | $U$ | $J$ |  | $U$ | $J$ |  | $U$ | $J$ |  | $U$ | $J$ |  | $U$ | $J$
0.80 | 6.62 | 1.00 |  | 0.62 | 0.20 |  | 2.93 | 0.49 |  | 7.41 | 1.08 |  | 0.57 | 0.20 |  | 2.80 | 0.48 |  | 8.04 | 1.20 |  | 0.49 | 0.20 |  | 2.29 | 0.44
0.84 | 6.35 | 1.00 |  | 0.57 | 0.19 |  | 2.77 | 0.48 |  | 7.12 | 1.07 |  | 0.53 | 0.19 |  | 2.65 | 0.48 |  | 7.74 | 1.19 |  | 0.45 | 0.19 |  | 2.16 | 0.44
0.88 | 6.09 | 0.99 |  | 0.53 | 0.18 |  | 2.62 | 0.47 |  | 6.84 | 1.06 |  | 0.49 | 0.18 |  | 2.51 | 0.47 |  | 7.44 | 1.18 |  | 0.42 | 0.18 |  | 2.04 | 0.43
0.92 | 5.85 | 0.98 |  | 0.50 | 0.18 |  | 2.48 | 0.47 |  | 6.58 | 1.05 |  | 0.46 | 0.18 |  | 2.38 | 0.46 |  | 7.16 | 1.17 |  | 0.39 | 0.18 |  | 1.93 | 0.42
0.96 | 5.62 | 0.97 |  | 0.47 | 0.17 |  | 2.36 | 0.46 |  | 6.33 | 1.04 |  | 0.43 | 0.17 |  | 2.26 | 0.46 |  | 6.90 | 1.16 |  | 0.36 | 0.17 |  | 1.82 | 0.42
1.00 | 5.40 | 0.96 |  | 0.44 | 0.16 |  | 2.24 | 0.45 |  | 6.10 | 1.03 |  | 0.40 | 0.16 |  | 2.14 | 0.45 |  | 6.65 | 1.15 |  | 0.34 | 0.16 |  | 1.73 | 0.41
1.04 | 5.20 | 0.95 |  | 0.41 | 0.16 |  | 2.13 | 0.45 |  | 5.87 | 1.03 |  | 0.38 | 0.16 |  | 2.04 | 0.44 |  | 6.41 | 1.15 |  | 0.32 | 0.16 |  | 1.64 | 0.40
1.08 | 5.00 | 0.94 |  | 0.39 | 0.15 |  | 2.03 | 0.44 |  | 5.66 | 1.02 |  | 0.35 | 0.15 |  | 1.94 | 0.44 |  | 6.18 | 1.14 |  | 0.30 | 0.15 |  | 1.56 | 0.40
cDFTa |  $U$=4.7, $J$=0.8 |  |  $U$=4.8, $J$=0.9 |  |  $U$=5.2, $J$=0.9
cRPAb |  $U$=5.5, $J$=0.6 |  |  $U$=5.7, $J$=0.7 |  |  $U$=6.6, $J$=0.7

  * 1

${}^{\text{a}}$The results of constrained DFT implemented in LAPW framework
from Ref.[83]

  * 3

${}^{\text{b}}$The results of constrained RPA in maximally localized Wannier
functions from Ref.[55]

Although a
$\lambda$=$1.00\text{\,}\mathrm{B}\mathrm{o}\mathrm{h}\mathrm{r}^{-1}$ seems
to work well for all TMOs, we don’t expect it to work in general cases,
because a fixed $\lambda$ value means that the obtained $U$, $J$ values only
depend on the atomic species and the chosen NAO basis sets, but not on the
chemical environment. For a generally applicable scheme, $\lambda$ should
reflect the chemical environment of the system. The averaged $\bar{\lambda}$
introduced in Eq. (51) depends on the electron density of the system via the
Thomas-Fermi screening model, and thus accounts for the chemical environment
in a natural way. In Ref. [60], it has been shown that such a scheme (termed
as LSCC there) yields rather good $U$, $J$ values within the LAPW framework.
Note that, within such a scheme, the averaged screening parameter
$\bar{\lambda}$ varies during the self-consistent iterations and hence so do
the $U$, $J$ values, until the convergence is reached. We also implemented the
LSCC scheme in ABACUS, and the self-consistently determined $\bar{\lambda}$
values for MnO, FeO, CoO, and NiO are 1.581,1.627,1.649 and 1.677 Bohr-1,
respectively. However, the $\bar{\lambda}$ parameters lead to too short
screening lengths for our NAO basis sets and the calculated $U$, $J$ values
are too small. In practice, we find that introducing a scaling factor of 0.625
can reduce the $\bar{\lambda}$ value to a range (around 1.0 Bohr-1) that
yields physically reasonable $U$, $J$ values. The reason that an additional
scaling factor is needed here, compared to the original LAPW-based LSCC
implementation, is that the local NAOs used here is more extended than the
local orbitals in the LAPW framework, which are restricted within the muffin-
tin sphere. Hence the computed interaction parameters within the NAO framework
will be smaller than the LAPW case if the same screening parameter is used.

In Table 6, we present the energy differences between the AFM and
ferromagnetic (FM) states of four materials, as computed by PBE and PBE+$U$.
The experimental results and the results reported in Ref. [60] are included
for comparison. The PBE+$U$ with a scaled $\bar{\lambda}$ parameter, obtained
from the Thomas Fermi model, yields results that show satisfactory agreement
with experiment and the LSCC results. The accuracy of the results is a factor
of two better than that of PBE.

Table 6:  The energy differences (in meV) between the AFM and FM states of four materials calculated by PBE and PBE+$U$ with $U$ and $J$ determined with a scaled screening parameter (see the text), and compared to the experimental values and the calculated results reported in Ref. [60]. In LSCC, the $U$, $J$ parameters are obtained without rescaling the $\bar{\lambda}$ parameter whereas in our work, a scaling factor of 0.625 is used. |  PBE |  PBE [60] |  LSCC [60] |  PBE+$U$ |  Expt.
---|---|---|---|---|---
|  (this work) |  |  |  (this work) |
MnO |  -157.3 |  -152 |  -99.8 |  -84.5 |  -62[84]
NiO |  -258.2 |  -261 |  -107 |  -145.4 |  -112[85]
MnF2 |  -64.6 |  -60.5 |  -28.3 |  -22.2 |  -15.2[86]
NiF2 |  -85.9 |  -69.6 |  -20.1 |  -31.7 |  -13.8[86]

##  V Summary

We present a detailed formulation of the DFT+U method within the framework of
NAO basis set. The key in this formulation is to use a symmetrized Mulliken
charge projector, constructed in terms of the most localized $d$ or $f$
orbital basis functions, to project a correlated local subspace out of the
full KS orbital space. We implemented such a scheme within the ABACUS code
package, and our implementation allows not only self-consistent electronic
structure calculations with or without including the SOC effect, but also
enables force and stress calculations. The efficacy of our formalism and
implementation has been demonstrated for the prototypical TMOs and IrO2.
Furthermore, we tested the scheme for computing the $U$, $J$ parameters based
on a screened Yukawa potential, and found that, while a fixed screening
parameter works for all TMOs, determining such a parameter from the electron
density via the Thomas-Fermi model leads to an underestimation of the $U$, $J$
values. However, this issue can be fixed by introducing a universal rescaling
parameter to increase the screening length. We believe that the experience
gained in the present work will be very helpful for developing Hubbard-type
local correction scheme within the atomic-orbital basis set framework.

###### Acknowledgements.

This work is supported by National Natural Science Foundation of China (Grant
Nos. 12134012, 11874335, 21873005) and the Strategic Priority Research Program
of Chinese Academy of Sciences (Grant No. XDPB25). We thank Dr. Wenshuai Zhang
for generating the optimized NAO basis sets used in the present work,
available at the official ABACUS website [87]. Parts of the calculations are
done on the supercomputing system in the Supercomputing Center of USTC. *

##  Appendix A The influence of the choice of local orbitals on DFT+U

In this appendix, we investigate two issues about employing the local atomic
orbitals to construct the projector within the NAO basis set framework, which
have not been elaborated in the main text. The first issue, which has already
been mentioned in Sec IV.5, is whether all basis orbitals belonging to
characteristic d/f angular moment channel need to be included in the DFT+U
correction. The second issue is the influence of the shape of the local
orbitals on the DFT+U results.

Table 7:  Partial occupation numbers (sum of spin-up and spin-down channels) of Mn d orbitals of MnO as calculated by PBE. Results obtained using DZP (two $d$ orbitals) and TZDP (three $d$ orbitals) are presented. |  |  $d_{xy}$ |  $d_{xz}$ |  $d_{z^{2}}$ |  $d_{yz}$ |  $d_{x^{2}-y^{2}}$
---|---|---|---|---|---|---
DZP |  zeta=1 |  1.372 |  1.372 |  0.773 |  1.372 |  0.773
|  zeta=2 |  0.012 |  0.012 |  0.007 |  0.012 |  0.007
TZDP |  zeta=1 |  1.355 |  1.343 |  0.728 |  1.474 |  0.728
|  zeta=2 |  0.030 |  0.032 |  0.028 |  0.002 |  0.032
|  zeta=3 |  0.002 |  0.002 |  -0.006 |  0.002 |  -0.006

To address the first question, we perform a partial occupation analysis of the
$d$ orbitals, taking MnO as an example. Table 7 presents the local occupation
numbers of all Mn d orbitals of MnO, as given by Eq. (16). For completeness
both DZP and TZDP basis sets have been used, and in the latter case there are
three $d$ orbitals for Mn. From Table 7, it can be clearly seen that local
occupation numbers of the second or third d orbital are significantly smaller
than the first d one. This means that the second or third $d$ orbitals
contribute little to the top valence states, and the majority of the
correlated local subspace is described by the first $d$ orbital. From the
point of view of real-space locality, the first (and innermost) d orbital,
i.e. zeta 1 in Table 7, is most localized and has no node. The other d
orbitals are designed to be orthogonal to the first one, so that they have
nodes and are more delocalized (cf. Fig. 6). Such behavior is also reflected
in two-electron integrals of the screened Coulomb potential. As is shown in
Table 6, the on-site Coulomb interaction and exchange energies of the second d
orbital are smaller by one order of magnitude compared to the first d one.
Both the orbital occupation and interaction parameter analyses suggest that we
can most likely neglect the on-site Coulomb correction to the second or higher
d orbital.

Figure 6: Projected DOS of NiO where different on-site Coulomb interactions
are employed to different d/f orbital:(a) and (b) DZP basis, (c) and (d) TZDP
basis. The upper and lower panels correspond to up-spin and down-spin. The
zero of energy is set to be Fermi level (vertical dashed line).

Despite the observation that most probably one only needs to apply the Hubbard
$U$ correction to the most localized orbital for multi-zeta NAO basis sets, we
nevertheless also checked what if the $U$ correction is added to all $d$
orbitals. In Fig. 6, the PDOS of NiO by using DZP and TZDP basis set on the
spin-up polarized nickel atoms and their nearest oxygen atoms are shown.
Figure. 6(a) and 6(c) are the results of standard DFT+U, where $\bar{U}$ is
set to be 5.5 eV and the on-site Coulomb interaction correction is only
applied to the first d orbitals. For comparison, We then apply on-site Coulomb
interaction corrections to all d and polarization f orbitals where the
corresponding $\bar{U}$ value are determined by the Yukawa potential scheme as
described in Sec. II.5 with screened parameter $\lambda$ fixed at
$1.00\text{\,}\mathrm{B}\mathrm{o}\mathrm{h}\mathrm{r}^{-1}$. The obtained
results are depicted in Fig. 6(b) and Fig. 6(d) for the DZP and TZDP basis
sets, respectively. In Fig. 6 the label d/f-i means the i-th d/f orbital and
the values behind them are the corresponding $\bar{U}$ values. As clearly
shown in Fig. 6, the PDOS results undergo little changes for both DZP and TZDP
basis sets if the $U$ corrections are added to all d and polarization f
orbitals. This result validates our DFT+$U$ projection scheme that only the
first, most localized correlated orbital needs to be included in the
construction of the projector.

Figure 7: Radial functions of 3d basis of Mn atom in ABACUS

Finally, we briefly discuss the possible influence of the shape of the local
atomic orbitals used in the projector construction on the DFT+U calculations.
From the discussion in Sec. II.1, one can see that the KS wavefunctions and
the local orbitals control the local occupation matrix and hence directly
affect the results of DFT+U calculations. Since all properly worked DFT codes
must yield nearly the same wavefunctions in real space, provided that the same
pseudopotential is used, the DFT+U results mainly depend on the choice of
local correlated orbitals. In Fig. 7, we plot the radial functions of Mn $d$
orbitals of the DZP basis sets, employed in ABACUS (upper panel) and OpenMX
(lower panel) calculations. For both codes, the first Mn d orbitals (the
nodeless ones) are localized within 1.5 Å around the nucleus, with a sharp
peak positioned around $r=$0.5\text{\,}\mathrm{\SIUnitSymbolAngstrom}$$. The
local occupation numbers given by the Mulliken charge projector in OpenMX are
1.447 for $t_{2g}$ orbitals and 0.412 for $e_{g}$ orbitals, in comparison to
1.372 and 0.773 as given by ABACUS (cf. Table 7). Such difference is a
manifestation of the difference in the radial shape of the first $d$ functions
between the two codes, as plotted in Fig. 7. As such, the DFT+$U$
implementations in ABACUS and OpenMX yield qualitatively similar but
quantitatively noticeably different results, as demonstrated in Fig. 1.

## References

  * [1] P. Hohenberg and W. Kohn.  Inhomogeneous electron gas.  Phys. Rev., 136:B864–B871, 1964.
  * [2] W. Kohn and L. J. Sham.  Self-Consistent Equations Including Exchange and Correlation Effects.  Phys. Rev., 140(4A):A1133–A1138, 1965.
  * [3] W. Kohn.  Nobel Lecture: Electronic structure of matter—wave functions and density functionals.  Rev. Mod. Phys., 71(5):1253–1266, oct 1999.
  * [4] Aron J. Cohen, P. Mori-Sanchez, and Weitao Yang.  Insights into Current Limitations of Density Functional Theory.  Science, 321(5890):792–794, 2008.
  * [5] Chen Li, Xiao Zheng, Neil Qiang Su, and Weitao Yang.  Localized orbital scaling correction for systematic elimination of delocalization error in density functional approximations.  Natl. Sci. Rev., 5(2):203–215, 2018.
  * [6] Neil Qiang Su, Chen Li, and Weitao Yang.  Describing strong correlation with fractional-spin correction in density functional theory.  Proc. Natl. Acad. Sci., 115(39):9678–9683, 2018.
  * [7] Neil Qiang Su, Aaron Mahler, and Weitao Yang.  Preserving Symmetry and Degeneracy in the Localized Orbital Scaling Correction Approach.  J. Phys. Chem. Lett., 11(4):1528–1535, 2020.
  * [8] J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz.  Density-functional theory for fractional particle number - derivative discontinuities of the energy.  Phys. Rev. Lett., 49(23):1691, 1982.
  * [9] Vladimir Anisimov and Yuri Izyumov.  Electronic Structure of Strongly Correlated Materials.  Springer Series in Solid-State Sciences. Springer Berlin Heidelberg, Berlin, Heidelberg, 2010.
  * [10] Burak Himmetoglu, Andrea Floris, Stefano de Gironcoli, and Matteo Cococcioni.  Hubbard-corrected DFT energy functionals: The LDA+U description of correlated systems.  Int. J. Quantum Chem., 114(1):14–49, jan 2014.
  * [11] J. P. Perdew and Alex Zunger.  Self-interaction correction to density-functional approximations for many-electron systems.  Phys. Rev. B, 23(10):5048, 1981.
  * [12] Axel D. Becke.  Density‐functional thermochemistry. iii. the role of exact exchange.  J. Chem. Phys., 98(7):5648, 1993.
  * [13] Jochen Heyd, Gustavo E. Scuseria, and Matthias Ernzerhof.  Hybrid functionals based on a screened coulomb potential.  J. Chem. Phys., 118(18):8207, 2003.
  * [14] Walter Metzner and Dieter Vollhardt.  Correlated Lattice Fermions in $d=\infty$ Dimensions.  Phys. Rev. Lett., 62(3):324–327, jan 1989.
  * [15] Antoine Georges, Gabriel Kotliar, Werner Krauth, and Marcelo J. Rozenberg.  Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions.  Rev. Mod. Phys., 68(1):13–125, 1996.
  * [16] Vladimir I. Anisimov, F. Aryasetiawan, and A. I. Lichtenstein.  First-principles calculations of the electronic structure and spectra of strongly correlated systems: the LDA + U method.  J. Phys. Condens. Matter, 9(4):767–808, 1997.
  * [17] A. I. Lichtenstein and M. I. Katsnelson.  Ab initio calculations of quasiparticle band structure in correlated systems: LDA++ approach.  Phys. Rev. B, 57(12):6884–6895, 1998.
  * [18] G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti.  Electronic structure calculations with dynamical mean-field theory.  Rev. Mod. Phys., 78(3):865–951, 2006.
  * [19] K. Held.  Electronic structure calculations using dynamical mean field theory.  Adv. Phys., 56(6):829–926, 2007.
  * [20] J. Hubbard.  Electron correlations in narrow energy bands.  J. Proc. R. Soc. A, 276(1365):238, 1963.
  * [21] Martin C. Gutzwiller.  Effect of correlation on the ferromagnetism of transition metals.  Phys. Rev. Lett., 10(5):159, 1963.
  * [22] Junjiro Kanamori.  Electron correlation and ferromagnetism of transition metals.  Prog. Theor. Phys., 30(3):275, 1963.
  * [23] Vladimir I. Anisimov, Jan Zaanen, and Ole K. Andersen.  Band theory and Mott insulators: Hubbard U instead of Stoner I.  Phys. Rev. B, 44(3):943–954, 1991.
  * [24] V. I. Anisimov, I. V. Solovyev, M. A. Korotin, M. T. Czyżyk, and G. A. Sawatzky.  Density-functional theory and NiO photoemission spectra.  Phys. Rev. B, 48(23):16929–16934, 1993.
  * [25] Vladimir I. Anisimov, F. Aryasetiawan, and A. I. Lichtenstein.  First-principles calculations of the electronic structure and spectra of strongly correlated systems: the LDA+U method.  J. Phys. Condens. Matter, 9(4):767–808, 1997.
  * [26] A. B. Shick, A. I. Liechtenstein, and W. E. Pickett.  Implementation of the LDA+U method using the full-potential linearized augmented plane-wave basis.  Phys. Rev. B, 60(15):10763–10769, 1999.
  * [27] https://elk.sourceforge.io/.
  * [28] O. Bengone, M. Alouani, P. Blöchl, and J. Hugel.  Implementation of the projector augmented-wave LDA+U method: Application to the electronic structure of NiO.  Phys. Rev. B, 62(24):16392–16401, dec 2000.
  * [29] B. Amadon, F. Jollet, and M. Torrent.  $\gamma$ and $\beta$ cerium: LDA+U calculations of ground-state parameters.  Phys. Rev. B, 77(15):155104, 2008.
  * [30] Paolo Giannozzi, Stefano Baroni, Nicola Bonini, and Matteo Calandra.  QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials.  J. Phys. Condens. Matter, 21(39):395502, 2009.
  * [31] Matteo Cococcioni and Stefano de Gironcoli.  Linear response approach to the calculation of the effective interaction parameters in the LDA+U method.  Phys. Rev. B, 71(3):035105, 2005.
  * [32] B. Delley.  From molecules to solids with the dmol3 approach.  J. Chem. Phys., 113(18):7756, 2000.
  * [33] K. Koepernik and H. Eschrig.  Full-potential nonorthogonal local-orbital minimum-basis band-structure scheme.  Phys. Rev. B, 59(3):1743, 1999.
  * [34] José M. Soler, Emilio Artacho, Julian D. Gale, Alberto García, Javier Junquera, Pablo Ordejón, and Daniel Sánchez-Portal.  The SIESTA method for ab initio order- N materials simulation.  J. Phys. Condens. Matter, 14(11):2745–2779, mar 2002.
  * [35] T. Ozaki, H. Kino, J. Yu, M. Han, N. Kobayashi, M. Ohfuti, F. Ishii, and T. Ohwaki.  User’s manual of OpenMX, http://www.openmx-square.org, 2008.
  * [36] Volker Blum, Ralf Gehrke, Felix Hanke, Paula Havu, Ville Havu, Xinguo Ren, Karsten Reuter, and Matthias Scheffler.  Ab initio molecular simulations with numeric atom-centered orbitals.  Comput. Phys. Commun., 180(11):2175, 2009.
  * [37] Pengfei Li, Xiaohui Liu, Mohan Chen, Peize Lin, Xinguo Ren, Lin Lin, Chao Yang, and Lixin He.  Large-scale ab initio simulations based on systematically improvable atomic basis.  Comput. Mater. Sci., 112:503–517, 2016.
  * [38] Xinguo Ren, Patrick Rinke, Volker Blum, Jürgen Wieferink, Alexandre Tkatchenko, Andrea Sanfilippo, Karsten Reuter, and Matthias Scheffler.  Resolution-of-identity approach to hartree–fock, hybrid density functionals, rpa, mp2 andgwwith numeric atom-centered orbital basis functions.  New J. Phys., 14(5), 2012.
  * [39] Sergey V. Levchenko, Xinguo Ren, Jürgen Wieferink, Rainer Johanni, Patrick Rinke, Volker Blum, and Matthias Scheffler.  Hybrid functionals for large periodic systems in an all-electron, numeric atom-centered basis framework.  Comput. Phys. Commun., 192:60, 2015.
  * [40] P. Lin, X. Ren, and L. He.  Accuracy of localized resolution of the identity in periodic hybrid functional calculations with numerical atomic orbitals.  J. Phys. Chem. Lett., 11(8):3082, 2020.
  * [41] P. Lin, X. Ren, and L. He.  Efficient hybrid density functional calculations for large periodic systems using numerical atomic orbitals.  J. Chem. Theory Comput., 17(1):222, 2021.
  * [42] Muhammad N. Tahir and Xinguo Ren.  Comparing particle-particle and particle-hole channels of the random phase approximation.  Phys. Rev. B, 99(19):195149, 2019.
  * [43] Xinguo Ren, Florian Merz, Hong Jiang, Yi Yao, Markus Rampp, Hermann Lederer, Volker Blum, and Matthias Scheffler.  All-electron periodic $g_{0}w_{0}$ implementation with numerical atomic orbital basis functions: Algorithm and benchmarks.  Phys. Rev. Mater., 5(1):013807, 2021.
  * [44] https://siesta-project.org/siesta/.
  * [45] Myung Joon Han, Taisuke Ozaki, and Jaejun Yu.  O$\left(N\right)$ LDA+U electronic structure calculation method based on the nonorthogonal pseudoatomic orbital basis.  Phys. Rev. B, 73(4):045110, 2006.
  * [46] M. Kick, K. Reuter, and H. Oberhofer.  Intricacies of dft+u, not only in a numeric atom centered orbital framework.  J. Chem. Theory Comput., 15(3):1705, 2019.
  * [47] Richard M. Martin.  Electronic Structure.  Cambridge University Press, 2004.
  * [48] A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen.  Density-functional theory and strong interactions: Orbital ordering in Mott-Hubbard insulators.  Phys. Rev. B, 52(8):R5467–R5470, 1995.
  * [49] P. H. Dederichs, S. Blügel, R. Zeller, and H. Akai.  Ground States of Constrained Systems: Application to Cerium Impurities.  Phys. Rev. Lett., 53(26):2512–2515, dec 1984.
  * [50] M. R. Norman and A. J. Freeman.  Model supercell local-density calculations of the 3d excitation spectra in NiO.  Phys. Rev. B, 33(12):8896–8898, jun 1986.
  * [51] O. Gunnarsson, O. K. Andersen, O. Jepsen, and J. Zaanen.  Density-functional calculation of the parameters in the Anderson model: Application to Mn in CdTe.  Phys. Rev. B, 39(3):1708–1722, jan 1989.
  * [52] F. Aryasetiawan, M. Imada, A. Georges, G. Kotliar, S. Biermann, and A. I. Lichtenstein.  Frequency-dependent local interactions and low-energy effective models from electronic structure calculations.  Phys. Rev. B, 70(19):195104, nov 2004.
  * [53] Takashi Miyake and F. Aryasetiawan.  Screened Coulomb interaction in the maximally localized Wannier basis.  Phys. Rev. B, 77(8):085122, feb 2008.
  * [54] Takashi Miyake, Ferdi Aryasetiawan, and Masatoshi Imada.  Ab initio procedure for constructing effective models of correlated materials with entangled band structure.  Phys. Rev. B, 80(15):155134, oct 2009.
  * [55] R. Sakuma and F. Aryasetiawan.  First-principles calculations of dynamical screened interactions for the transition metal oxides MO ( M=Mn, Fe, Co, Ni).  Phys. Rev. B, 87(16):165118, 2013.
  * [56] Erik R. Ylvisaker, Warren E. Pickett, and Klaus Koepernik.  Anisotropy and magnetism in the LSDA+U method.  Phys. Rev. B, 79(3):035103, 2009.
  * [57] S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton.  Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+U study.  Phys. Rev. B, 57(3):1505–1509, 1998.
  * [58] Fredrik Bultmark, Francesco Cricchio, Oscar Grånäs, and Lars Nordström.  Multipole decomposition of LDA+U energy and its application to actinide compounds.  Phys. Rev. B, 80(3):035121, 2009.
  * [59] M. R. Norman.  Calculation of effective Coulomb interaction for Pr3+, U4+, and UPt3 .  Phys. Rev. B, 52(3):1421–1424, 1995.
  * [60] Yue Chao Wang and Hong Jiang.  Local screened Coulomb correction approach to strongly correlated d -electron systems.  J. Chem. Phys., 150(15):16–20, 2019.
  * [61] A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen.  Density-functional theory and strong interactions: Orbital ordering in Mott-Hubbard insulators.  Phys. Rev. B, 52(8):R5467–R5470, aug 1995.
  * [62] D. R. Hamann.  Optimized norm-conserving Vanderbilt pseudopotentials.  Phys. Rev. B, 88(8):085117, 2013.
  * [63] Martin Schlipf and François Gygi.  Optimization algorithm for the generation of ONCV pseudopotentials.  Comput. Phys. Commun., 196:36–44, 2015.
  * [64] Peter Scherpelz, Marco Govoni, Ikutaro Hamada, and Giulia Galli.  Implementation and Validation of Fully Relativistic GW Calculations: Spin–Orbit Coupling in Molecules, Nanocrystals, and Solids.  J. Chem. Theory Comput., 12(8):3523–3544, 2016.
  * [65] Mohan Chen, G-C Guo, and Lixin He.  Systematically improvable optimized atomic basis sets for ab initio calculations.  J. Phys. Condens. Matter, 22(44):445501, 2010.
  * [66] Wahyu Setyawan and Stefano Curtarolo.  High-throughput electronic band structure calculations: Challenges and tools.  Comput. Mater. Sci., 49(2):299–312, 2010.
  * [67] S. K. Panda, S. Bhowal, A. Delin, O. Eriksson, and I. Dasgupta.  Effect of spin orbit coupling and Hubbard U on the electronic structure of IrO2.  Phys. Rev. B, 89(15):155102, 2014.
  * [68] John P. Perdew, Kieron Burke, and Matthias Ernzerhof.  Generalized Gradient Approximation Made Simple.  Phys. Rev. Lett., 77(18):3865–3868, 1996.
  * [69] T. Ozaki.  Variationally optimized atomic orbitals for large-scale electronic structures.  Phys. Rev. B, 67(15):155108, 2003.
  * [70] T. Ozaki and H. Kino.  Numerical atomic basis orbitals from H to Kr.  Phys. Rev. B, 69(19):195113, may 2004.
  * [71] T. Ozaki and H. Kino.  Numerical atomic basis orbitals from H to Kr.  Phys. Rev. B, 69(19):1–19, 2004.
  * [72] Fabien Tran, Peter Blaha, Karlheinz Schwarz, and Pavel Novák.  Hybrid exchange-correlation energy functionals for strongly correlated electrons: Applications to transition-metal monoxides.  Phys. Rev. B, 74(15):155108, 2006.
  * [73] Yue-Chao Wang, Ze-Hua Chen, and Hong Jiang.  The local projection in the density functional theory plus U approach: A critical assessment.  J. Chem. Phys., 144(14):144106, apr 2016.
  * [74] L. Messick, W. C. Walker, and R. Glosser.  Direct and Temperature-Modulated Reflectance Spectra of MnO, CoO, and NiO.  Phys. Rev. B, 6(10):3941–3949, 1972.
  * [75] H.K. Bowen, D. Adler, and B.H. Auker.  Electrical and optical properties of FeO.  J. Solid State Chem., 12(3-4):355–359, jan 1975.
  * [76] R. J. Powell and W. E. Spicer.  Optical Properties of NiO and CoO.  Phys. Rev. B, 2(6):2182–2193, 1970.
  * [77] G. A. Sawatzky and J. W. Allen.  Magnitude and Origin of the Band Gap in NiO.  Phys. Rev. Lett., 53(24):2339–2342, 1984.
  * [78] S. Hüfner, J. Osterwalder, T. Riesterer, and F. Hulliger.  Photoemission and inverse photoemission spectroscopy of NiO.  Solid State Commun., 52(9):793–796, 1984.
  * [79] A. K. Cheetham and D. A O Hope.  Magnetic ordering and exchange effects in the antiferromagnetic solid solutions MnxNi1-xO.  Phys. Rev. B, 27(11):6964–6967, 1983.
  * [80] B. E. F. Fender, A. J. Jacobson, and F. A. Wedgwood.  Covalency Parameters in MnO, $\alpha$‐MnS, and NiO.  J. Chem. Phys., 48(3):990–994, 1968.
  * [81] W. L. Roth.  Magnetic Structures of MnO, FeO, CoO, and NiO.  Phys. Rev., 110(6):1333–1341, 1958.
  * [82] D. C. Khan and R. A. Erickson.  Magnetic Form Factor of Co++ Ion in Cobaltous Oxide.  Phys. Rev. B, 1(5):2243–2249, 1970.
  * [83] Hong Jiang, Ricardo I. Gomez-Abal, Patrick Rinke, and Matthias Scheffler.  First-principles modeling of localized d states with the GW@LDA+U approach.  Phys. Rev. B, 82(4):045108, 2010.
  * [84] G. Pepy.  Spin waves in MnO; from 4K to temperatures close to $T_{N}$.  J. Phys. Chem. Solids, 35(3):433–444, jan 1974.
  * [85] R. Shanker and R. A. Singh.  Analysis of the Exchange Parameters and Magnetic Properties of NiO.  Phys. Rev. B, 7(11):5000–5005, jun 1973.
  * [86] Xiaobing Feng and N. M. Harrison.  Magnetic coupling constants from a hybrid density functional with 35% Hartree-Fock exchange.  Phys. Rev. B, 70(9):092402, sep 2004.
  * [87] The ABACUS software webpage: http://abacus.ustc.edu.cn.