# Localized orbital scaling correction for periodic systems

Aaron Mahler  Duke University, Department of Physics, Durham, NC 27708  Jacob
Williams  Duke University, Department of Chemistry, Durham, NC 27708  Neil
Qiang Su  Department of Chemistry, Key Laboratory of Advanced Energy Materials
Chemistry (Ministry of Education) and Renewable Energy Conversion and Storage
Center (RECAST), Nankai University, Tianjin 300071, China  Duke University,
Department of Chemistry, Durham, NC 27708  Weitao Yang<EMAIL_ADDRESS>Duke University, Department of Chemistry, Durham, NC 27708  Duke University,
Department of Physics, Durham, NC 27708

(August 27, 2024)

###### Abstract

Density functional theory offers accurate structure prediction at acceptable
computational cost, but commonly used approximations suffer from
delocalization error; this results in inaccurate predictions of quantities
such as energy band gaps of finite and bulk systems, energy level alignments,
and electron distributions at interfaces. The localized orbital scaling
correction (LOSC) was developed to correct delocalization error by using
orbitals localized in space and energy. These localized orbitals span both the
occupied and unoccupied spaces and can have fractional occupations in order to
correct both the total energy and the one-electron energy eigenvalues. We
extend the LOSC method to periodic systems, in which the localized orbitals
employed are dually localized Wannier functions. In light of the effect of the
bulk environment on the electrostatic interaction between localized orbitals,
we modify the LOSC energy correction to include a screened Coulomb kernel. For
a test set of semiconductors and large-gap insulators, we show that the
screened LOSC (sLOSC) method consistently improves the band gap compared to
the parent density functional approximation.

##  I Introduction

The cost of solving the electronic Schrödinger equation scales exponentially
with the size of the system, exceeding the computational resources available
on the planet for any system larger than a few tens of electrons. Kohn (1999)
Density functional theory (DFT) sidesteps this exponential cost by treating
the electron density as the fundamental variable instead of computing the
wavefunction directly and by constructing an auxiliary noninteracting
reference system sharing the density of the physical system. Hohenberg and
Kohn (1964); Kohn and Sham (1965) Due to the accuracy attainable at a cost
only cubic in the number of electrons $N$, DFT has become a mainstay of
computational chemistry and materials science. Becke (1993); Lee  _et al._
(1988); Perdew  _et al._ (1996); Van Noorden  _et al._ (2014) While DFT is
exact in theory, the form of the universal exchange-correlation functional is
unknown, and density functional approximations (DFAs) must be used in
practice. Commonly used DFAs suffer from systematic delocalization and static
correlation errors. Cohen  _et al._ (2008a, 2012) The delocalization error
underlies the failure of DFAs to describe energy band gaps of finite and bulk
systems, energy level alignments, and electron distributions at interfaces.
Mori-Sánchez  _et al._ (2008) Overcoming delocalization error remains an
active and challenging research effort.

Connecting single-particle orbital energies $\epsilon$ to observable
quantities was another longstanding question in Kohn–Sham DFT. As a contrast,
Koopmans (1934) showed in 1934 that the Hartree–Fock ionization potential (IP)
and electron affinity (EA) are given under the frozen orbital aproximation by
the negative of the highest occupied and lowest unoccupied molecular orbital
eigenvalues respectively. A series of three results established a rigorous
connection for DFT.

First, Janak (1978) derived a link between the Kohn–Sham orbital energies
$\epsilon_{m}$ and the total energy $E$, viewed as a function of the orbital
occupation numbers $n_{m}$:

$\epsilon_{m}=\frac{\partial E}{\partial n_{m}}.$ (1)

However, $\partial E/\partial n_{m}$ was not yet linked to a physical
observable.

A few years later, Perdew, Parr, Levy, and Balduz Perdew  _et al._ (1982)
showed that $E$ is piecewise linear in the number of electrons $N$ when
computed with the exact functional; that is, for all $\lvert\delta\rvert\leq
1$, we have

$E(N+\delta)=\begin{cases}(1+\delta)E(N)-\delta E(N-1)&\delta<0,\\\
(1-\delta)E(N)+\delta E(N+1)&\delta\geq 0.\end{cases}$ (2)

This relationship, called the PPLB condition, connects the chemical potential
$\mu(N)=\partial E/\partial N$ to the IP and EA; observe that

$\mu(N)=\begin{cases}-I(N)=E(N)-E(N-1)&\partial N<0,\\\
-A(N)=E(N+1)-E(N)&\partial N>0.\end{cases}$ (3)

Finally, Cohen  _et al._ (2008b) proved that the chemical potential is given
by the partial derivative of $E$ with respect to the frontier orbital
eigenvalues,

$\mu(N)=\frac{\partial E}{\partial n_{f}}.$ (4)

Crucially, $f$ labels not only the highest unoccupied molecular orbital (HOMO)
if $\partial N<0$, but also the lowest unoccupied molecular orbital (LUMO) if
$\partial N>0$; this was the first time a physical meaning for the energy of
the Kohn–Sham LUMO was established. This result holds for any local functional
continuous in the electron density, as well as any nonlocal functional
continuous in the Kohn–Sham density matrix; in the latter case, the work also
extends Janak’s theorem to the eigenvalues from the generalized Kohn–Sham
equations.

Combining these three results, we see that

$\mu(N)=\begin{cases}-I(N)=\epsilon_{\tiny{\text{HOMO}}}&\partial N<0,\\\
-A(N)=\epsilon_{\tiny{\text{LUMO}}}&\partial N\geq 0.\end{cases}$ (5)

Thus, the frontier eigenvalues obtained from an $N$-electron DFT calculation
correspond rigorously to physically relevant quantities;Cohen  _et al._
(2008b) if the PPLB condition is obeyed and the functional predicts the exact
energies for $N-1$, $N$, and $N+1$ electrons, the correspondence is exact.

A feature derivable from these quantities is the fundamental or integer gap,
defined as the difference between the IP and the EA:

$\begin{split}E_{\text{gap}}^{\text{integer}}&=I-A\\\
&=E(N-1)-2E(N)+E(N+1).\end{split}$ (6)

$E_{\text{gap}}^{\text{integer}}$ quantifies the difference between positively
and negatively ionizing the system and is a crucial part of the accurate
modeling of semiconductor electronic structure. If the PPLB condition is
obeyed, Eqs. (3) and (5) also allow the gap to be computed from a single
$N$-electron calculation as the discontinuity in the chemical potential; in
this form, it is called the derivative gap, defined as

$\begin{split}E^{\text{deriv}}_{\text{gap}}&=\left.\frac{\partial E}{\partial
N}\right\rvert_{+}-\left.\frac{\partial E}{\partial N}\right\rvert_{-}\\\
&=\epsilon_{\text{\tiny{LUMO}}}-\epsilon_{\text{\tiny{HOMO}}}.\end{split}$ (7)

If the PPLB condition is obeyed, the derivative gap and the integer gap are
equal.Cohen  _et al._ (2008b) In bulk systems with periodic boundary
conditions, the PPLB condition is satisfied by any DFA continuous in the
Kohn–Sham density or density matrix, regardless of systematic errors in its
definition; Mori-Sánchez  _et al._ (2008) thus, the fundamental gap of bulk
systems can be predicted by the (generalized) Kohn–Sham orbital gap, for
functionals continuous in the density (density matrix), as in Eq. (7). In
finite systems, however, the PPLB condition is not in general obeyed, and the
gap computed from Eq. (7) may differ from that computed by calculating the
$(N\pm 1)$-electron energies to obtain the integer gap as in Eq. (6), the
$\Delta$SCF method.

###  I.1 Delocalization error

The delocalization error has a dramatic size-dependent manifestation. In
finite systems, standard DFAs fail to obey the PPLB linearity condition, so
the derivative gap is not equal to the integer gap. This is due to the error
in the approximate exchange-correlation functional, which nearly always yields
$E$ convex in $N$, underestimating the piecewise linearity prescribed by the
PPLB condition. This convex deviation has been identified as the cause for an
unphysical smearing of the electron density in space, as well as
underestimation of the total energy in a delocalized electron density; thus,
we may identify it with delocalization error, as exhibited in small systems.
In bulk systems, the delocalized nature of the orbitals produces a total
energy linear with respect to fractional charge, yielding no deviation from
the PPLB condition; however, delocalization error manifests as an incorrect
slope of the $E(N)$ line at integer $N$. Mori-Sánchez  _et al._ (2008)

The effects of delocalization error include the underestimation of band gaps
and reaction barriers, Merkle  _et al._ (1992) undervaluation of dissociation
curves, Zhang and Yang (1998); Mori-Sánchez  _et al._ (2006a); Perdew  _et
al._ (2007) overestimation of conductance and polarizability, Körzdörfer  _et
al._ (2008) and incorrect energy level alignment and charge transfer across
interfaces. Jones and Gunnarsson (1989); Becke (2014) To capture the full
derivative discontinuity and hence the band gap, it has been shown that the
exact functional, whether local or nonlocal, cannot be a differentiable
functional of the electron density or of the Kohn-Sham density matrix. Mori-
Sánchez  _et al._ (2009); Yang  _et al._ (2012a) To reduce the systematic
delocalization error, many approaches have been developed, including range-
separated functionals, Savin and Flad (1995); Savin (1996); Iikura  _et al._
(2001); Yanai  _et al._ (2004); Vydrov and Scuseria (2006); Chai and Head-
Gordon (2008); Baer  _et al._ (2010) the screened range-separated hybrid
functional, Wing  _et al._ (2021) self-interaction error corrected
functionals, Perdew and Zunger (1981); Mori-Sánchez  _et al._ (2006a, b);
Perdew  _et al._ (2008); Schmidt  _et al._ (2014); Pederson  _et al._ (2014);
Schmidt and Kümmel (2016); Yang  _et al._ (2017) Koopmans-compliant
functionals, Borghi  _et al._ (2014); Colonna  _et al._ (2019) and generalized
transition state methods, Anisimov and Kozhevnikov (2005) along with related
developments using localized Wannier functions. Ma and Wang (2016)

The localized orbital scaling correction (LOSC) method was developed to
eliminate delocalization error systematically. Li  _et al._ (2018); Su  _et
al._ (2020) Previous incarnations of LOSC were implemented for molecular
systems with real orbitals and the boundary condition
$\lim_{|\mathbf{r}|\to\infty}\rho(\mathbf{r})=0$. They accurately model IP,
EA, photoemission spectra, dissociation curves, and polarizabilities, as well
as restore size-consistency. Li  _et al._ (2018); Su  _et al._ (2020); Mei
_et al._ (2020, 2021a) In this work, we extend LOSC to periodic boundary
conditions (PBCs) and complex orbitals. Additionally, we introduce a screened
Coulomb interaction to the LOSC energy correction to enable the accurate
computation of bulk system band structures.

###  I.2 Periodic boundary conditions

In PBCs, the eigenfunctions of the single-particle Hamiltonian are known as
Bloch orbitals; they satisfy
$h_{\text{s}}\lvert\psi_{n}^{\mathbf{k}}\rangle=\epsilon_{n}^{\mathbf{k}}\lvert\psi_{n}^{\mathbf{k}}\rangle$.
The Bloch orbitals are also eigenfunctions of the unit cell translation
operator, so they take the form
$\lvert\psi_{n}^{\mathbf{k}}\rangle=e^{i\mathbf{k}\cdot\mathbf{r}}\lvert
u_{n}^{\mathbf{k}}\rangle$, where $\lvert u_{n}^{\mathbf{k}}\rangle$ has the
periodicity of the unit cell and $\mathbf{k}$ is a point in the Brillouin zone
(the reciprocal-space unit cell). Bloch (1929); Ashcroft and Mermin (2003) The
Bloch orbitals obey the normalization convention
$\langle\psi_{m}^{\mathbf{q}}|\psi_{n}^{\mathbf{k}}\rangle=\delta(\mathbf{k}-\mathbf{q})\delta_{mn}$,
where $\langle
f|g\rangle=\int_{\mathcal{D}}d\mathbf{r}\,\overline{f}(\mathbf{r})g(\mathbf{r})$.
Here, $\delta(\mathbf{k})$ is the Dirac delta distribution, $\delta_{mn}$ is
the Kronecker delta, and $\overline{f}$ is the complex conjugate of $f$. The
domain of integration $\mathcal{D}$ is the periodic unit for the functions
being integrated; for Bloch orbitals, $\mathcal{D}=\mathbb{R}^{3}$. The
$\lvert u_{n}^{\mathbf{k}}\rangle$ are orthonormal in the band index $n$ at a
fixed $\mathbf{k}$-point in reciprocal space: that is, $\langle
u_{m}^{\mathbf{k}}|u_{n}^{\mathbf{k}}\rangle=\delta_{mn}$, where the inner
product integrates over one unit cell. Note that we assume closed-shell
systems in this work.

The single-particle density can be represented in real space by the occupied
Bloch orbitals as

$\rho_{s}(\mathbf{r})=\sum_{n}^{\text{occ}}\frac{V}{(2\pi^{3})}\int_{\text{BZ}}d\mathbf{k}\,\lvert\psi_{n}^{\mathbf{k}}(\mathbf{r})\rvert^{2},$
(8)

where $V$ is the volume of the unit cell and the integral is over the first
Brillouin zone. Since the Hamiltonian is diagonal in $\mathbf{k}$, we can
solve for the Bloch orbitals in reciprocal space, requiring diagonalization
only in one unit cell. In practice, the Brillouin zone is sampled with a
finite number of points; in this work we use a Monkhorst–Pack mesh centered at
the origin $\Gamma$ of the Brillouin zone. Monkhorst and Pack (1976) Thus,
integrals over the Brillouin zone become equally weighted sums over the
$\mathbf{k}$-mesh

$\frac{V}{(2\pi)^{3}}\int
d\mathbf{k}\,f(\mathbf{k})\mapsto\frac{1}{N_{k}}\sum_{\mathbf{k}}f(\mathbf{k}),$
(9)

where $N_{k}$ is the number of $\mathbf{k}$-points in the mesh. Using a
Monkhorst–Pack mesh centered at $\Gamma$ yields Bloch orbitals having the
periodicity of an unfolded supercell comprised of $N_{k}$ primitive unit
cells; this supercell is referred to as the Born–von Karman cell. Ashcroft and
Mermin (2003) The Bloch orbitals then obey the normalization convention
$\langle\psi_{m}^{\mathbf{q}}|\psi_{n}^{\mathbf{k}}\rangle=N_{k}\delta_{\mathbf{k}\mathbf{q}}\delta_{mn}$,
where the integral is over the Born–von Karman cell.

##  II Methods

The LOSC method consists of two steps. First, we find orbitals that are
spatially localized while remaining associated with specific energy ranges.
Next, we compute a curvature matrix modeling the magnitude of the deviation
from linearity. This is combined with the fractionally occupied localized
orbitals to correct the convex deviation of $E(N)$ from linearity at non-
integer $N$, as well as incorrect total energies at integer $N$. Both steps
are implemented as post-processing after a converged self-consistent field
calculation.

###  II.1 Localization

The wave-like nature of the Bloch orbitals prohibits them from being spatially
localized. In order to obtain a state that is localized in space, the discrete
Fourier transform of the Bloch orbitals is used to produce Wannier functions
Wannier (1937)

$\begin{split}\lvert
w_{n}^{\mathbf{R}}\rangle&=\frac{1}{N_{k}}\sum_{\mathbf{k}}e^{-i\mathbf{k}\cdot\mathbf{R}}\lvert\psi_{n}^{\mathbf{k}}\rangle.\end{split}$
(10)

Wannier functions inherit the periodicity of the Born–von Karman cell, and are
indexed by electron bands $n$ and unit cells $\mathbf{R}$ in the supercell.
They are symmetric under translation by unit cell vectors, so that
$w_{n}^{\mathbf{R}}(\mathbf{r})=w_{n}^{\mathbf{0}}(\mathbf{r}-\mathbf{R})$;
$\mathbf{R}=\mathbf{0}$ is referred to as the home unit cell. There is a
unitary, or gauge, freedom at each $\mathbf{k}$-point in the choice of Bloch
orbitals that comprise a Wannier function, so we can define _generalized_
Wannier functions Cloizeaux (1963)

$\begin{split}\lvert
w_{i}^{\mathbf{R}}\rangle&=\frac{1}{N_{k}}\sum_{\mathbf{k}}e^{-i\mathbf{k}\cdot\mathbf{R}}\sum_{n}U_{ni}^{\mathbf{k}}\lvert\psi_{n}^{\mathbf{k}}\rangle\\\
&=\frac{1}{N_{k}}\sum_{\mathbf{k}}e^{-i\mathbf{k}\cdot\mathbf{R}}\lvert\phi_{i}^{\mathbf{k}}\rangle,\end{split}$
(11)

where we refer to $\lvert\phi_{i}^{\mathbf{k}}\rangle$ as a transformed Bloch
orbital (TBO). From now on, we will refer to generalized Wannier functions as
Wannier functions.

The gauge freedom $U^{{}^{\mathbf{k}}}$ in the TBOs can be chosen such that
the resulting set of Wannier functions have advantageous properties. In order
to obtain localization in space, Marzari and Vanderbilt suggested choosing
$U^{{}^{\mathbf{k}}}$ that minimize the Wannier functions’ spatial variance
$\langle\Delta r^{2}\rangle_{i}=\langle
r^{2}\rangle_{i}-\langle\mathbf{r}\rangle^{2}_{i}$, where $\langle
x\rangle_{i}=\langle w_{i}^{\mathbf{0}}\lvert x\rvert
w_{i}^{\mathbf{0}}\rangle$; the resulting orbitals are called maximally
localized Wannier functions. Marzari and Vanderbilt (1997) In molecules, the
scheme of minimizing spatial variance is referred to as Foster–Boys
localization. Foster and Boys (1960) However, constructing maximally localized
Wannier functions from both valence and conduction bands is physically ill-
motivated. Because bands far apart in energy can mix freely and the Bloch
bands form a complete basis, adding more virtual bands will result in
increasing spatial localization of the maximally localized Wannier functions,
with a corresponding loss of information about the energy dispersion of the
bands.

In order to preserve locality in energy while maintaining spatial
localization, enabling simultaneous treatment of the occupied and unoccupied
spaces, we choose the Wannier gauge that minimizes a cost function considering
both energy and spatial variance:

$F=(1-\gamma)\sum_{i}\langle\Delta r^{2}\rangle_{i}+\gamma
C\sum_{i}\langle\Delta h_{\text{s}}^{2}\rangle_{i},$ (12)

where $0\leq\gamma\leq 1$ and in the units used here
$C=$1\text{\,}\mathrm{a_{0}}^{2}\mathrm{/}\mathrm{eV}^{2}$$, where
$\mathrm{a_{0}}$ is the Bohr radius. This cost function was first proposed by
Gygi  _et al._ (2003) for computations sampling the Brillouin zone only at
$\Gamma$ and implemented for such systems by Giustino and Pasquarello (2006).
It was used for LOSC in molecules Su  _et al._ (2020) to treat system
symmetries and degeneracies more robustly than the original localization,
which used soft energy windows; Li  _et al._ (2018) in molecular LOSC, the
localized orbitals are called _orbitalets_. We recently extended $F$ to
systems with $N_{k}\geq 1$; Mahler  _et al._ (2022) we refer to such orbitals
as dually localized Wannier functions (DLWFs). These formulations show how the
combination of occupied and unoccupied spaces can be localized simultaneously
to produce Wannier functions that are localized in both space and energy. This
construction is critical for addressing delocalization error in finite systems
because it allows for dynamic localization in the resulting orbitals; the
orbitals can qualitatively and quantitatively differ depending on the geometry
of the system. Li  _et al._ (2018); Su  _et al._ (2020) In keeping with the
principle of universality in functional development, we use the same mixing
parameter in Eq. (12), setting $\gamma=0.47714$. The value of $\gamma$ has
important implications for the LOSC method; see Section VI of the
Supplementary Information of Su  _et al._ (2020). Setting $\gamma=0$, for
instance, yields maximally localized Wannier functions, Marzari and Vanderbilt
(1997) while $\gamma=1$ yields DLWFs that are pure Fourier transforms (up to
$\mathbf{k}$-dependent phases) of the Kohn–Sham bands.

Note that we have not proven that a unique global minimum of $F$ exists; in
practice, $F$ has a fairly rugged landscape of solutions, and we have observed
multiple local minima. We choose the DLWFs yielding the smallest total cost.
Different DLWFs can produce somewhat different sLOSC corrections, with
eigenvalues varying by up to a few tenths of an $\mathrm{eV}$. The problem of
multiple minima of $F$ has also been observed in molecular LOSC, Mei  _et al._
(2022) but was not found to be the dominant source of error. An additional
question worth exploring is the effect of symmetry breaking, such as that due
to perturbations of the crystal lattice, on the localization procedure.

The compromise between spatial and energy localization and the inclusion of
unoccupied orbitals are key to producing localized orbitals that can address
delocalization error while retaining size-consistency. For example, the DFA
HOMO and LUMO of H${}_{2}^{+}$ at (or near) the dissociation limit are
delocalized over the whole molecule; since they are (nearly) degenerate, there
exists a unitary freedom in the subspace spanned by both. Due to the symmetry
of the system, we expect to obtain two separate H0.5+ fragments; the (small
or) vanishing gap means that any choice of $\gamma<1$ in Eq. (12) will result
in half-occupied orbitals localized on each H atom. This is the physical
motivation for a localization scheme that minimizes the spatial variance of
occupied and unoccupied orbitals while allowing only orbitals that are close
in energy to mix. Li  _et al._ (2018)

###  II.2 Energy corrections

The deviation from energy linearity with respect to fractional charges is
characteristically quadratic in most exchange-correlation functionals. Zheng
_et al._ (2011); Hait and Head-Gordon (2018); Li  _et al._ (2018) To restore
compliance with the PPLB condition for small finite systems, the global
scaling correction (GSC) was developed. GSC corrects the total energy by an
amount quadratic in the occupation numbers of the canonical molecular
orbitals. Zheng  _et al._ (2011); Mei  _et al._ (2021b) This method is
effective at correcting the systematic deviation from the PPLB condition for
systems with fractional charges and leads to accurate prediction of
quasiparticle energies as the eigenvalues from the resulting one-electron
Hamiltonian. However, GSC is applicable only for systems of small and moderate
size; the convex deviation of conventional DFAs from the piecewise linearity
prescribed by the PPLB condition decreases with increasing system size, and
the delocalization error manifests instead as underestimated ground-state
energies for integer systems and incorrect linear $E_{\text{gs}}(N)$ curves
with wrong slopes at the bulk limit. Mori-Sánchez  _et al._ (2008) The
localized orbital scaling correction (LOSC) applies its energy correction
adaptively by the construction of localized orbitals, allowing systematic and
size-consistent correction of delocalization error. Li  _et al._ (2018); Su
_et al._ (2020) In this section, we discuss the extension of LOSC to periodic
systems.

A basic quantity in LOSC is the density matrix in the basis of DLWFs; its
elements are occupations

$\lambda_{ij}^{\mathbf{T}\mathbf{R}}=\langle
w_{i}^{\mathbf{T}}\lvert\rho_{s}\rvert w_{j}^{\mathbf{R}}\rangle.$ (13)

The occupations between all pairs of DLWFs are used to remove quadratic
deviations, while the diagonal terms are used to restore linearity. The energy
correction defined by LOSC for each unit cell is given by

$\Delta
E^{\text{LOSC}}=\frac{1}{2N_{k}}\sum_{\mathbf{T}\mathbf{R}}\sum_{ij}\widetilde{\kappa}_{ij}^{\mathbf{T}\mathbf{R}}{\lambda}_{ij}^{\mathbf{T}\mathbf{R}}(\delta_{ij}^{\mathbf{T}\mathbf{R}}-\overline{\lambda}_{ij}^{\mathbf{T}\mathbf{R}}),$
(14)

where
$\delta_{ij}^{\mathbf{0}\mathbf{R}}=\delta_{ij}\delta_{\mathbf{0}\mathbf{R}}$;
$\widetilde{\kappa}$ models the curvature of the deviation from linearity.

The diagonal terms in the energy correction are proportional to
$\lambda_{ii}^{\mathbf{T}\mathbf{R}}-\lvert\lambda_{ii}^{\mathbf{T}\mathbf{R}}\rvert^{2}$;
thus, if a DLWF has integer occupancy (implying
$\lambda_{ij}^{\mathbf{T}\mathbf{R}}=0$ whenever $i\neq j$ or
$\mathbf{T}\neq\mathbf{R}$), then the energy correction due to that DLWF will
also be zero.

The matrix $[\lambda_{ij}^{\mathbf{T}\mathbf{R}}]$ of occupations between the
DLWFs is the discrete Fourier transform of the occupation matrix between the
TBOs. As such, it is positive semidefinite and Hermitian, and

$\operatorname{tr}_{c}[\lambda_{ij}^{\mathbf{T}\mathbf{R}}]=\frac{1}{N_{k}}\sum_{\mathbf{k}}N_{f}^{\mathbf{k}},$
(15)

where $\operatorname{tr}_{c}$ denotes the trace per unit cell and
$N_{f}^{\mathbf{k}}$ is the number of electrons below the Fermi energy at
$\mathbf{k}$.

Following Su  _et al._ (2020), the elements of the curvature matrix are given
by

$\begin{split}\widetilde{\kappa}_{ij}^{\mathbf{T}\mathbf{R}}=&\operatorname{erf}(8S_{ij}^{\mathbf{T}\mathbf{R}})\sqrt{\kappa_{ii}^{\mathbf{T}\mathbf{T}}\kappa_{jj}^{\mathbf{R}\mathbf{R}}}\\\
+&\operatorname{erfc}(8S_{ij}^{\mathbf{T}\mathbf{R}}){\kappa}_{ij}^{\mathbf{T}\mathbf{R}}.\end{split}$
(16)

Here, $\operatorname{erfc}(r)=1-\operatorname{erf}(r)$ is the complementary
error function. $S_{ij}^{\mathbf{T}\mathbf{R}}$ is the absolute overlap
between DLWFs,

$S_{ij}^{\mathbf{T}\mathbf{R}}=\int
d\mathbf{r}\,\sqrt{\rho_{i}^{\mathbf{T}}(\mathbf{r})\rho_{j}^{\mathbf{R}}(\mathbf{r})},$
(17)

where $\rho_{i}^{\mathbf{T}}(\mathbf{r})=\lvert
w_{i}^{\mathbf{T}}(\mathbf{r})\rvert^{2}$ is a DLWF’s charge density. The
matrix elements ${\kappa}_{ij}^{\mathbf{T}\mathbf{R}}$ in Eq. (16) are given
by

$\kappa_{ij}^{\mathbf{T}\mathbf{R}}=J[\rho_{i}^{\mathbf{T}},\rho_{j}^{\mathbf{R}}]-X[\rho_{i}^{\mathbf{T}},\rho_{j}^{\mathbf{R}}],$
(18)

with

$\displaystyle J[\rho_{i}^{\mathbf{T}},\rho_{j}^{\mathbf{R}}]$
$\displaystyle=\iint
d\mathbf{r}\,d\mathbf{r}^{\prime}\,\rho_{i}^{\mathbf{T}}(\mathbf{r})\rho_{j}^{\mathbf{R}}(\mathbf{r}^{\prime})K(\lvert\mathbf{r}-\mathbf{r}^{\prime}\rvert),$
(19a) $\displaystyle X[\rho_{i}^{\mathbf{T}},\rho_{j}^{\mathbf{R}}]$
$\displaystyle=\tau\frac{2C_{\text{X}}}{3}\int
d\mathbf{r}\,\Big{[}\rho_{i}^{\mathbf{T}}(\mathbf{r})\rho_{j}^{\mathbf{R}}(\mathbf{r})\Big{]}^{2/3}.$
(19b)

In the above, $K(r)=1/r$ is the Coulomb kernel,
$C_{X}=\frac{3}{4}(\frac{6}{\pi})^{1/3}$ is the Dirac exchange constant, Dirac
(1930) and $\tau=6(1-2^{\text{-}1/3})\approx 1.2378$ is a nonempirical
parameter.Li  _et al._ (2018) The derivation of how this correction restores
the PPLB condition can be found in the supplementary data of Li  _et al._
(2018). The use of $\widetilde{\kappa}$ instead of $\kappa$ was introduced
because the cost function in Eq. (12) can induce discontinuous jumps between
localization characters during molecular dissociation.Su  _et al._ (2020) The
diagonal elements of $\widetilde{\kappa}$ and $\kappa$ are equal, so when
$\lambda_{ii}\in\\{0,1\\}$ the corrections from $\widetilde{\kappa}$ and
$\kappa$ are the same. In practice,
$X[\rho_{i}^{\mathbf{T}},\rho_{j}^{\mathbf{R}}]$ term is evaluated using
numerical integration on a grid of real-space points. The Coulomb term
$J[\rho_{i}^{\mathbf{T}},\rho_{j}^{\mathbf{R}}]$ is evaluated in a plane wave
basis, detailed in Sec. II.3.

Applying the extension of Janak’s theoremJanak (1978) to the generalized Kohn-
Sham theory, Cohen  _et al._ (2008b) the LOSC energy correction of Eq. (14)
yields corrections to the Bloch orbital energy eigenvalues
$\epsilon_{n}^{\mathbf{k}}$ given by

$\Delta\epsilon_{n}^{\mathbf{k}}=\sum_{i}\widetilde{\kappa}_{ii}^{\mathbf{0}\mathbf{0}}\left(\frac{1}{2}-\lambda_{ii}^{\mathbf{0}\mathbf{0}}\right)\lvert
U_{ni}^{\mathbf{k}}\rvert^{2}\\\
-\sum_{\mathbf{R}i\neq\mathbf{0}j}\widetilde{\kappa}_{ij}^{\mathbf{0}\mathbf{R}}\mathrm{Re}\left\\{\lambda_{ij}^{\mathbf{0}\mathbf{R}}e^{i\mathbf{k}\cdot\mathbf{R}}\overline{U}_{ni}^{\mathbf{k}}U_{nj}^{\mathbf{k}}\right\\}.$
(20)

Consider the diagonal corrections given by the first summand in Eq. (20).
There is no correction to the eigenvalue when a DLWF is half-occupied
($\lambda_{ii}=\frac{1}{2}$); on the other hand, the correction is maximal
when it is completely occupied or unoccupied. Li  _et al._ (2018) observed
that the slopes of the quadratic DFA and the correct linear $E(N)$ curves
agree at half-integer $N$. Since the frontier orbital energy corresponds to
this slope, we see that accurate frontier orbital energies are given by half-
occupied frontier orbitals. The LOSC correction to the orbital energies
arrives naturally at this conclusion, additionally agreeing with Slater
transition state theory. Slater and Phillips (1974); Slater (1974)

We may also view LOSC as a correction to the Kohn–Sham Hamiltonian. It is
given by the functional derivative of the energy correction with respect to
the density operator under the frozen orbital approximation,

$\Delta v=\left.\frac{\delta\Delta
E^{\text{LOSC}}}{\delta\rho_{s}}\right\rvert_{\\{w_{i}^{\mathbf{R}}\\}},\\\ $
(21)

and can be written in operator form as

$\Delta
v=\frac{1}{2}\sum_{ij,\mathbf{T}\mathbf{R}}\widetilde{\kappa}_{ij}^{\mathbf{T}\mathbf{R}}\,\left(\frac{\delta_{ij}^{\mathbf{T}\mathbf{R}}}{2}-\lambda_{ij}^{\mathbf{T}\mathbf{R}}\right)\lvert
w_{i}^{\mathbf{T}}\rangle\\!\langle w_{j}^{\mathbf{R}}\rvert+\text{h.c.}$ (22)

(See the Supplemental Material for details on this derivation.) The correction
to the $n^{\text{th}}$ Bloch orbital eigenvalue is then given by
$\Delta\epsilon_{n}^{\mathbf{k}}=\langle\psi_{n}^{\mathbf{k}}\lvert\Delta
v\rvert\psi_{n}^{\mathbf{k}}\rangle$.

In practice, the energy corrections are applied to disentangled Bloch
orbitals. The conduction bands of most systems cannot be formed into sets of
bands that do not cross anywhere in the Brillouin zone, a condition referred
to as band entanglement. In order to obtain a finite set of bands for
localization and energy correction, we use the disentanglement procedure
outlined by Souza  _et al._ (2001) This procedure obtains $N_{w}$ bands from a
set of $N_{b}\geq N_{w}$ Bloch orbitals at each $\mathbf{k}$-point, chosen
such that the subspace spanned by the disentangled bands is as smooth as
possible in $\mathbf{k}$. To correct the band gap of semiconductors and
insulators, we include sufficiently many virtual bands in the construction of
the Wannier functions to converge the localization of the frontier bands (that
is, the valence band maximum and conduction band minimum). Mahler  _et al._
(2022) We find that $N_{b}=N_{\text{occ}}+3N_{\text{coord}}$ and
$N_{w}=N_{\text{occ}}+2N_{\text{coord}}$, where $N_{\text{occ}}$ is the number
of occupied bands per unit cell and $N_{\text{coord}}$ is the coordination
number of the lattice, are sufficient; see the Supplemental Material for
details. The $N_{w}$ disentangled Bloch bands yields $N_{w}$ DLWFs per unit
cell, so there are $N_{w}N_{k}$ DLWFs in the Born–von Karman supercell on
which they are periodic. The energy corrections for the $N_{w}$ disentangled
Bloch orbitals at each $\mathbf{k}$-point are implemented using Eq. (20).

In this work, we restrict our attentions to closed-shell systems. Extending
the LOSC method to spin-polarized materials is accomplished by finding the
corrections from the spin-up and spin-down DLWFs independently and summing
them to obtain $\Delta E^{\text{LOSC}}$; this functionality is planned for the
next version of LOSC. However, treating the strong correlation common to open-
shell materials brings its own set of challenges beyond the scope of this
work. We discuss them briefly in Sec. IV below.

###  II.3 Coulomb integrals

Accurate calculation of the Coulomb interaction
$J[\rho_{i}^{\mathbf{T}},\rho_{j}^{\mathbf{R}}]$ is needed for LOSC to restore
the PPLB condition. In PBCs, a plane-wave basis is typically employed. The
Coulomb energy is diagonal in this basis, and the double integral required in
real space collapses to a single sum over basis vectors $\mathbf{G}$:

$J[\rho_{i}^{\mathbf{T}},\rho_{j}^{\mathbf{R}}]=\sum_{\mathbf{G}}\frac{4\pi}{G^{2}}\,\overline{\rho}_{i}^{\mathbf{T}}(\mathbf{G})\rho_{j}^{\mathbf{R}}(\mathbf{G}),$
(23)

where $G=\lvert\mathbf{G}\rvert$. However, this sum converges only for neutral
charge distributions, for which the $\mathbf{G}=\mathbf{0}$ term vanishes. The
DLWF densities are individually charged, so ignoring the divergent term coming
from the net charge will significantly underestimate the Coulomb energy. There
are many methods to evaluate the Coulomb energy for charged densities in the
plane-wave basis accurately, including those of Makov and Payne (1995);
Kantorovich (1999); Dabo  _et al._ (2008), and Li and Dabo (2011). We choose
the spherical cutoff method, Onida  _et al._ (1995); Jarvis  _et al._ (1997);
Rozzi  _et al._ (2006) truncating the Coulomb kernel in Eq. (19a) at a cutoff
radius $R_{c}$; this is taken to be half the length of the shortest Born–von
Karman supercell lattice vector, ensuring that the Coulomb interactions
between the a DLWF density and its images in neighboring supercells are zero.
Thus, the spherical cutoff Coulomb kernel is

$K_{c}(r;R_{c})=\begin{cases}1/r&r<R_{c}\\\ 0&r\geq R_{c},\end{cases}$ (24)

which has Fourier coefficients

$K_{c}(G;R_{c})=\begin{cases}\dfrac{4\pi}{G^{2}}[1-\text{cos}(GR_{c})]&G\neq
0\\\\[6.00006pt] 2\pi R_{c}^{2}&G=0.\end{cases}$ (25)

Observe that $K_{c}(G;R_{c})$ does not diverge for any $\mathbf{G}$. As long
as the pair of DLWF densities in Eq. (19a) lie in a sphere of radius $R_{c}$,
the spherical cutoff method is also accurate in highly anisotropic unit cells,
unlike schemes such as that of Makov and Payne. Kantorovich (1999) We enforce
this containment condition in practice by checking that each DLWF density is
well contained in a volume spanned by half of each Born–von Karman supercell
lattice vector, and only compute curvature elements between pairs of DLWF
densities that have centers closer together than $R_{c}$. We evaluate the
Coulomb integrals on the unfolded supercell in the plane-wave basis, which
requires a fast Fourier transform (FFT) of the DLWF densities on the
supercell.

###  II.4 Screening

Applying LOSC with a bare Coulomb interaction leads to severe overcorrection
of semiconductors’ band gaps in PBCs. However, this is not surprising; we
anticipate an effect of the other electrons in the lattice on
$J[\rho_{i}^{\mathbf{T}},\rho_{j}^{\mathbf{R}}]$. As shown by highly accurate
methods such as $GW$, a screened Coulomb interaction is required to model the
interaction between electrons in a periodic system accurately. Hedin (1965);
Martin  _et al._ (2016) Recently, Mei and coworkers also found that the
deviation from linearity of the total energy as a function of canonical
orbital occupations is given to second order by a screened interaction. Mei
_et al._ (2021b) We model the screening phenomenologically, attenuating the
long-range $1/r$ behavior of the spherical cutoff Coulomb interaction by a
complementary error function. This modifies the Coulomb kernel to read

$K_{s}(r;R_{c},\alpha)=\begin{cases}\operatorname{erfc}(\alpha r)/r&r<R_{c}\\\
0&r\geq R_{c},\end{cases}$ (26)

where $\alpha$ is a screening parameter. We choose the $\alpha$ that best
reproduces the experimental band gaps of a test set of semiconductors and
insulators. For $r$ larger than the screening radius $\alpha^{-1}$,
$K_{s}(r;R_{c},\alpha)$ decays exponentially instead of as $1/r$. The Fourier
coefficients of $K_{s}$ are

$K_{s}(G;R_{c},\alpha)=\begin{cases}\dfrac{4\pi}{G^{2}}\left[1-\cos(GR_{c})\operatorname{erfc}(\alpha
R_{c})-e^{-(G/2\alpha)^{2}}\mathrm{Re}\left\\{\operatorname{erf}\left(\alpha
R_{c}+\dfrac{iG}{2\alpha}\right)\right\\}\right]&G\neq 0\\\\[11.00008pt] 2\pi
R_{c}^{2}+\pi\operatorname{erf}(\alpha
R_{c})\left(\alpha^{-2}-2R_{c}^{2}\right)-2\sqrt{\pi}e^{-(\alpha
R_{c})^{2}}R_{c}/\alpha&G=0.\end{cases}$ (27)

The error function is unbounded for complex arguments, overflowing double-
precision floating-point numbers even for relatively small $G$. Thus, we
evaluate $K_{s}(G;R_{c},\alpha)$ with a scaled form of $\operatorname{erf}z$
called the Faddeeva function, implemented in the numerically stable ACM
Algorithm 916. Zaghloul and Ali (2012); Zaghloul (2016) For details, see the
Supplemental Material.

In principle, the screening is system-dependent. Improved accuracy would be
attainable by setting its value to best reproduce each material’s band gap.
However, the phenomenological screening model of sLOSC does not enable doing
so while retaining predictive ability. This would require a rigorously
screened Coulomb (or Hartree-exchange-correlation) interaction based on the
linear response function
$\chi(\mathbf{r},\mathbf{r}^{\prime})=\delta\rho_{s}(\mathbf{r})/\delta
v(\mathbf{r}^{\prime})$. Ab initio screening of this kind appears in the
extensions of the GSC method to hybrid functionals Zheng  _et al._ (2013) and
in the following exploration of orbital relaxation on GSC, Zhang  _et al._
(2015), the GSC2 method, Mei  _et al._ (2021b) as well as in recent work on
Koopmans-compliant functionals. Colonna  _et al._ (2018); Nguyen  _et al._
(2018); Colonna  _et al._ (2022) For small, finite systems, the delocalization
error is quantified by $\partial^{2}E/\partial n_{i}^{2}$, where $n_{i}$ is
the occupation number of the Kohn–Sham orbital $\lvert\psi_{i}\rangle$;
analytical expressions for $\partial^{2}E/\partial n_{i}^{2}$ were derived in
Yang  _et al._ (2012b) In sLOSC, linear-response screening would very likely
increase the accuracy, but at substantial computational cost to compute the
nonlocal $\chi(\mathbf{r},\mathbf{r}^{\prime})$.

##  III Results

We use the PBE functional Perdew  _et al._ (1996) for the parent DFA
calculations, with optimized norm-conserving Vanderbilt pseudopotentials
Hamann (2013) generated by PseudoDojo. van Setten  _et al._ (2018) Both self-
consistent field (SCF) and non-SCF calculations are carried out in the Quantum
ESPRESSO code suite. Giannozzi  _et al._ (2009, 2017)

The energy cutoff for the fast Fourier transform is set to
$100\text{\,}\mathrm{Ry}$ for wavefunctions and $400\text{\,}\mathrm{Ry}$ for
densities. The Brillouin zone is sampled with Monkhorst–Pack meshes centered
at $\Gamma$, which is necessary for the Wannier functions to be periodic on
the Born–von Karman supercell. For SCF calculations, we use a
16$\times$16$\times$16 $\mathbf{k}$-mesh, while the other calculations are
performed on 6$\times$6$\times$6 grids. The localization step of LOSC is
implemented in a modified fork of the wannier90 code, Mostofi  _et al._
(2014); Pizzi  _et al._ (2020) and the energy correction as module to a fork
of Quantum ESPRESSO.

To determine an optimal screening parameter $\alpha$, we minimize the mean
absolute percent error (MAPE) on the SC/40 set of semiconductors with
experimentally available band gaps Heyd  _et al._ (2005) together with six
additional large-gap insulators. The experimental band gaps studied range from
$0.23\text{\,}\mathrm{eV}21.7\text{\,}\mathrm{eV}$. We find that $\alpha=$
$0.15\text{\,}\mathrm{a_{0}}^{-1}$ achieves the lowest MAPE; coincidentally,
this value is numerically equal to the screening parameter used in the HSE
density functional. Heyd  _et al._ (2005) As shown in Fig. 1, LOSC with
Coulomb screening (sLOSC) yields marked improvement of the band gap for the
test set in comparison with the parent functional. It is also apparent that
unscreened LOSC overcorrects the band gaps; indeed, it is less accurate than
the parent functional. The performance of sLOSC in molecules is better than
the parent functional, but unscreened LOSC achieves the best performance in
molecular systems (see Table 1). The Supplemental Material details the
variation in performance of screened LOSC for both bulk systems and molecules
with the creening parameter $\alpha$.

Figure 1: Comparison of experimental band gaps with those calculated by PBE (o), sLOSC ($\boldsymbol{\times}$), and unscreened LOSC (+). The inset shows systems with an experimental band gap less than $5\text{\,}\mathrm{eV}$.  Table 1: Mean absolute percent error of the band gap for PBC and molecular test sets. For details on the systems tested, see the Supplemental Material.  Method | PBE | LOSC | sLOSC
---|---|---|---
PBC | $47.5\%$ | $158.6\%$ | $19.7\%$
Molecule | $79.8\%$ | $10.1\%$ | $43.6\%$

The band structures of sLOSC and of the parent functional are shown for the
small-gapped semiconductor silicon in Fig. 2 and the larger-gapped insulator
lithium fluoride in Fig. 3. They use the disentangled band structures, which
are numerically indistinguishable from the true band structure at and below
the conduction band minimum for the parent functional. Wannier interpolation
in the same DLWF basis as that used in sLOSC is used to find the energy at the
points in the Brillouin zone not explicitly treated by the localization and
energy correction. The sLOSC correction to the band structure comes largely
from the more localized DLWFs, for which there is a larger Coulomb self-energy
$J[\rho_{i}^{\mathbf{R}},\rho_{i}^{\mathbf{R}}]$. Because of this, sLOSC
mostly corrects the energy of the occupied bands (which we observe to
correspond closely to the occupied DLWFs in semiconductors); it affects the
virtual bands much less. More work on molecule-surface and surface-surface
interactions is required to determine whether LOSC yields correct energy level
alignment and whether the larger correction to the valence bands is physically
meaningful.

Figure 2: Band structure of silicon under the PBE functional (dashes) and
sLOSC (solid). The Fermi energy of the PBE calculation was
$6.23\text{\,}\mathrm{eV}$, while the PBE with LOSC Fermi energy was
$5.49\text{\,}\mathrm{eV}$. Figure 3: Band structure of lithium fluoride under
the PBE functional (dashes) and sLOSC (solid). The Fermi energy of the PBE
calculation was $0.97\text{\,}\mathrm{eV}$, while the PBE with LOSC Fermi
energy was $-3.52\text{\,}\mathrm{eV}$. The core states are not included in
the figure.

##  IV Discussion

We have shown that despite the simple form of its phenomenological Coulomb
screening, sLOSC systematically corrects the band gap error associated with
the parent functional for materials spanning a large range of band gaps.
Screening improves the correction of delocalization error in bulk systems, but
degrades the accuracy of molecular systems’ band gaps relative to unscreened
LOSC; however, sLOSC still offers better band gaps than those computed by the
parent functional. One key remaining challenge is to model the curvature more
accurately for all systems; we expect that linear response of the electron
density, used by Mei  _et al._ (2021b) for accurate screening of the Kohn–Sham
orbitals, could be used to find the exact expression for
$\partial^{2}E/\partial\lambda_{ij}$. This would alleviate the error imposed
by modeling $\kappa$ as a difference between Coulomb repulsion and Dirac
exchange.

We implement the energy correction as a post-processing step to a self-
consistent calculation; such corrections are accurate when the change in
electron density is small and hence the total energy correction is small. For
every system considered in this work, $\Delta E^{\text{LOSC}}$ does not exceed
three parts in $10^{5}$. The corresponding change to the density for such
systems is also expected to be minimal. LOSC can also be implemented self-
consistently; Mei  _et al._ (2020) this can correct the delocalization error
of the total density, improving the accuracy of LOSC for systems with large
total energy corrections. A self-consistent implementation of sLOSC could be
necessary for the accurate computation of heterogeneous and interfacial
systems. Since delocalization error leads to incorrect charge distributions,
the sLOSC energy correction is likely to be larger, and self-consistently
correcting the delocalized density is expected to yield better orbital
energies.

Work is ongoing to implement sLOSC for spin-polarized materials and to
investigate its treatment of metals. The DLWFs of gapless systems constructed
from Bloch bands near the Fermi energy are expected to have occupations
$\lambda_{ii}$ close to $\frac{1}{2}$, which means that the sLOSC correction
to those eigenvalues will be small. While there may be changes to the overall
band structure, it is likely that such systems will remain gapless.

This may not hold in semimetals, whose valence and conduction bands cross only
in a small volume (or a single point) of the Brillouin zone. Metals, on the
other hand, have one band that crosses the Fermi energy. sLOSC can open a gap
in systems the DFA predicts to be semimetals; this occurs with the smallest-
gapped system in our test set, InSb. Thus, it is not certain that true
semimetals would remain so after the sLOSC correction. In addition, the
treatment of strong correlation due to (near) degeneracy of spin states and
the inclusion of topological or spin-orbit effects are beyond the scope of the
current work. A modification of molecular LOSC to include fractional spins was
developed in Su  _et al._ (2018); it could possibly be extended to bulk
materials as well.

###  IV.1 Comparison with other methods

####  IV.1.1 DFT+U(+V)

(s)LOSC is related to the DFT+U Anisimov  _et al._ (1991); Cococcioni and de
Gironcoli (2005) method for correcting delocalization error. The kinship can
be seen in the similarity of the sLOSC energy correction, Eq. (14), to the
rotationally invariant DFT+U correction O’Regan (2012)

$\Delta
E_{\text{DFT+U}}=\frac{1}{2}\sum_{\ell,\sigma}\operatorname{tr}\left[\boldsymbol{n}_{\ell\sigma}(1-\boldsymbol{n}_{\ell\sigma})\right]U_{\ell},$
(28)

where $U_{\ell}$ is the effective Hubbard parameter for orthonormal local
orbital (LO) $\ell$, combining atom and orbital indices, and
$\boldsymbol{n}_{\ell\sigma}$ is the LO occupation matrix. Both offer an
adjustment to the total energy quadratic in the occupation of the LOs. The
energy correction of DFT+U comes only from interactions between LOs on the
same atom, although DFT+U+V Campo and Cococcioni (2010) extends this to
interactions between atoms, analogous to the off-diagonal curvature elements
$\widetilde{\kappa}_{ij}$ of (s)LOSC. However, the LOs of DFT+U(+V) are static
(usually being $d$ and $f$ orbitals on transition metal centers), while the
DLWFs of sLOSC dynamically localize based on the gauge set by the cost
function $F$. Thus, where DFT+U(+V) must recompute the effective Hubbard
parameter for every perturbation of the crystal structure or molecular
geometry, the size of the correction in (s)LOSC follows from the DLWFs.

It is worth noting that, while DFT+U(+V) does not explicitly include energy
localization in its construction, the Hubbard correction applies primarily to
the Kohn–Sham orbitals that have the most overlap with the LOs; viewed another
way, the (spatially) localized orbitals that have the most energy-local
character (via their large overlap with energy eigenstates). Himmetoglu  _et
al._ (2014) In particular, the $d$ and $f$ atomic orbitals of transition
metals correspond closely to flat bands in reciprocal space, which carry some
energy information implicitly. However, they are independent of the of the
system’s geometry. In contrast, the LOs of LOSC are dynamic: the orbitals can
change with the geometric structure of the system. This is key to their
utility in finite systems. In compact structures near equilibrium, the LOSC
LOs can replicate the Kohn–Sham canonical orbitals, while becoming localized
as chemical bonds are stretched. This allows the LOSC total energy correction
to change with the geometry, as seen in Li  _et al._ (2018) and Su  _et al._
(2020)

Neither sLOSC nor DFT+U are suitable for solving the analogue of
delocalization error for systems with fractional spin. Cohen  _et al._
(2008b); Mori-Sánchez  _et al._ (2009) For molecules, fractional-spin LOSC
(FSLOSC) Su  _et al._ (2018) extends the original LOSC method to this case;
the judiciously modified DFT (jmDFT) method Bajaj  _et al._ (2017, 2019) does
the same for DFT+U.

####  IV.1.2 Koopmans-compliant functionals

Koopmans-compliant functionals Dabo  _et al._ (2009, 2010) mitigate
delocalization error by enforcing the PPLB linearity condition directly: in
the Koopmans integral (KI) formulation, Borghi  _et al._ (2014)

$\Delta
E_{\text{KI}}=\sum_{i}\alpha_{i}\left[f_{i}\eta_{i}-\int_{0}^{f_{i}}ds_{i}\,\langle\phi_{i}\lvert
h_{s}(s_{i})\rvert\phi_{i}\rangle\right].$ (29)

Here $\alpha_{i}$ is an orbital-dependent screening function based on the
relaxation of the LOs $\phi_{i}$; $f_{i}$ is the (fractional) occupation of
$\phi_{i}$; $\eta_{i}=\int_{0}^{1}ds_{i}\,\langle\phi_{i}\lvert
h_{\text{PZ}}(s_{i})\rvert\phi_{i}\rangle$, integrating the Perdew–Zunger
self-interaction corrected Kohn–Sham Hamiltonian, Perdew and Zunger (1981)
gives the linearized slope of the energy with respect to $f_{i}$; and the last
term computes the nonlinearity in $E$ that is replaced by $f_{i}\eta_{i}$.
Like (s)LOSC, Koopmans-compliant functionals are dependent on the choice of
localized orbitals. Borghi  _et al._ (2015) In extended systems, localized
orbitals are necessary for a Koopmans-compliant correction to have any effect,
Nguyen  _et al._ (2018) and screening has also been found to effect
improvements in band gap calculation. Colonna  _et al._ (2018, 2022)

An advantage of sLOSC over the Koopmans-compliant functionals for extended
systems is that the DLWFs treat the valence and conduction bands together; on
the one hand, DLWFs are empirically robust to increasing the number of
conduction bands from which they are constructed, and on the other, sLOSC can
in principle be applied to metals without additional modification. For gapped
systems, the energy localization inherent in the DLWF cost function enforces
separation between the occupied and virtual electronic manifolds without
manual input. The system with the smallest gap in our analysis, indium
antimonide (InSb, experimental gap $0.23\text{\,}\mathrm{eV}$), which is
predicted to be gapless by the DFA (and whose sLOSC gap is
$0.260\text{\,}\mathrm{eV}$) has DLWF occupations
$\lambda_{ii}^{\mathbf{R}\mathbf{R}}\geq 0.98$ in the valence manifold and
$\leq 0.0074$ in the conduction manifold. It is conceivable that DLWFs could
also serve as effective LOs for Koopmans-compliant methods, even if their
subspaces corresponding to the valence and conduction bands are not
variational for total Koopmans-compliant energy.

####  IV.1.3 Additional methods

The Fermi–Löwdin orbital (FLO) self-interaction correction (SIC), Pederson
_et al._ (2014); Kao  _et al._ (2017); Jackson  _et al._ (2019) which has its
roots in the Perdew–Zunger (PZ) self-interaction correction method, Perdew and
Zunger (1981) also uses localized orbitals for an energy correction. However,
the self-interaction error treated by both FLOSIC and PZ-SIC is well-defined
only for one-electron systems; LOSC and its derivatives account for the many-
electron nature of delocalization error explicitly. Mori-Sánchez  _et al._
(2006a)

As mentioned in Su  _et al._ (2018), the generalized transition state method
Anisimov and Kozhevnikov (2005) and the Wannier-function method of Ma and Wang
Ma and Wang (2016) are effective at improving band gap predictions. However,
since they do not mix valence and conduction bands to create fractionally
occupied orbitals, they cannot change the total energy of the DFA calculation;
thus, they cannot restore size-consistency to DFAs and will not be able to
capture (for instance) molecular dissociation at the same time as improving
band gap predictions. This problem is shared by early Koopmans-compliant
methods, which used the Kohn–Sham orbitals as the $\phi_{i}$; it underlies the
observation of Nguyen  _et al._ (2018) that localized orbitals such as Wannier
functions are required for Koopmans compliance in extended systems.

###  IV.2 Computational efficiency

The sLOSC method as implemented in this work scales as $O(N_{w}^{3}N_{k})$ for
the localization step, $O(N_{w}^{2}N_{G})$ for the computation of curvature
elements, and $O(N_{w}N_{G}\log N_{G})$ for the FFT of the DLWF densities.
Here, $N_{G}$ is the number of plane waves in the unfolded supercell, which is
$N_{k}$ times the number of plane waves in the unit cell. Calculating the
curvature and energy corrections is the computational bottleneck for the
systems evaluated in this work, with wall times for each system reaching a few
hours using 16 threads on an Intel Xeon E5-2630v3 processor. The systems that
took the longest time were those with the largest number of core states, which
have no effect on frontier state corrections; these could be neglected if only
a correction to the band gap is desired. The running time was divided fairly
evenly between the computation of the matrix elements defined in Eqs. (17,
19a, 19b). We note that the size of the integration domain for these
quantities could be reduced from the full Born–von Karman supercell if the
relevant DLWF densities are contained in a smaller region. This is supported
by the fact that systems that had similar localizations for a
4$\times$4$\times$4 and 6$\times$6$\times$6 $\mathbf{k}$-mesh yielded very
similar energy corrections.

Along the same lines, we find that for the systems in the test set the
correction results are converged with a 6$\times$6$\times$6 $\mathbf{k}$-mesh,
but this is only necessary to achieve a converged localization. Certain
systems exhibit a qualitatively different localization with a smaller
4$\times$4$\times$4 $\mathbf{k}$-mesh; however, by decreasing the value of
$\gamma$ in the localization cost function $F$, a set of DLWFs qualitatively
similar to the 6$\times$6$\times$6 case can be obtained.

Some other methods that attempt to address delocalization error in bulk
calculations, such as the approach of Ma and Wang (2016) and the screened
range-separated hybrid functional, Wing  _et al._ (2021) rely on supercell
self-consistent calculations. These have cubic scaling in the number of
electrons, so an unfolded Born–von Karman supercell arising from $N_{k}$
$\mathbf{k}$-points sampling a unit cell with $N_{w}$ Wannier functions scales
as $O(N_{k}^{3}N_{w}^{3})$. Both of the aforementioned methods use Wannier
functions as a localized charge representation and rely on manually choosing
the Bloch orbitals to comprise the Wannier functions representing the frontier
of the occupied space. The sLOSC method uses DLWFs, which naturally supply
Wannier functions representing the frontier of the occupied and unoccupied
spaces without the need for manual energy windowing.

###### Acknowledgements.

A.M., J.Z.W., N.Q.S., and W.Y. acknowledge support from the National Science
Foundation (Grant No. CHE-1900338); A.M. was additionally supported by the
Molecular Sciences Software Institute Phase-II Software Fellowship, and J.Z.W.
by the National Institutes of Health (Grant No. 5R01GM061870).

## References

  * Kohn (1999) W. Kohn, Rev. Mod. Phys. 71, 1253 (1999), publisher: American Physical Society.
  * Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964), publisher: American Physical Society.
  * Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965), publisher: American Physical Society.
  * Becke (1993) A. D. Becke, The Journal of Chemical Physics 98, 5648 (1993).
  * Lee  _et al._ (1988) C. Lee, W. Yang, and R. G. Parr, Phys. Rev. B 37, 785 (1988), publisher: American Physical Society.
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Van Noorden  _et al._ (2014) R. Van Noorden, B. Maher, and R. Nuzzo, Nature News 514, 550 (2014), cg_type: Nature News Section: News Feature.
  * Cohen  _et al._ (2008a) A. J. Cohen, P. Mori-Sánchez, and W. Yang, Science 321, 792 (2008a), publisher: American Association for the Advancement of Science Section: Special Perspectives.
  * Cohen  _et al._ (2012) A. J. Cohen, P. Mori-Sánchez, and W. Yang, Chem. Rev. 112, 289 (2012), publisher: American Chemical Society.
  * Mori-Sánchez  _et al._ (2008) P. Mori-Sánchez, A. J. Cohen, and W. Yang, Phys. Rev. Lett. 100, 146401 (2008), publisher: American Physical Society.
  * Koopmans (1934) T. Koopmans, Physica 1, 104 (1934).
  * Janak (1978) J. F. Janak, Phys. Rev. B 18, 7165 (1978), publisher: American Physical Society.
  * Perdew  _et al._ (1982) J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz, PHYSICAL REVIEW LETTERS 49, 4 (1982).
  * Cohen  _et al._ (2008b) A. J. Cohen, P. Mori-Sánchez, and W. Yang, Phys. Rev. B 77, 115123 (2008b).
  * Merkle  _et al._ (1992) R. Merkle, A. Savin, and H. Preuss, J. Chem. Phys. 97, 9216 (1992), publisher: American Institute of Physics.
  * Zhang and Yang (1998) Y. Zhang and W. Yang, J. Chem. Phys. 109, 2604 (1998), publisher: American Institute of Physics.
  * Mori-Sánchez  _et al._ (2006a) P. Mori-Sánchez, A. J. Cohen, and W. Yang, J. Chem. Phys. 125, 201102 (2006a), publisher: American Institute of Physics.
  * Perdew  _et al._ (2007) J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, V. N. Staroverov, and J. Tao, Phys. Rev. A 76, 040501 (2007), publisher: American Physical Society.
  * Körzdörfer  _et al._ (2008) T. Körzdörfer, S. Kümmel, and M. Mundt, J. Chem. Phys. 129, 014110 (2008), publisher: American Institute of Physics.
  * Jones and Gunnarsson (1989) R. O. Jones and O. Gunnarsson, Rev. Mod. Phys. 61, 689 (1989), publisher: American Physical Society.
  * Becke (2014) A. D. Becke, J. Chem. Phys. 140, 18A301 (2014), publisher: American Institute of Physics.
  * Mori-Sánchez  _et al._ (2009) P. Mori-Sánchez, A. J. Cohen, and W. Yang, Phys. Rev. Lett. 102, 066403 (2009), publisher: American Physical Society.
  * Yang  _et al._ (2012a) W. Yang, A. J. Cohen, and P. Mori-Sánchez, J. Chem. Phys. 136, 204111 (2012a), publisher: American Institute of Physics.
  * Savin and Flad (1995) A. Savin and H.-J. Flad, International Journal of Quantum Chemistry 56, 327 (1995), _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/qua.560560417.
  * Savin (1996) A. Savin, in  _Theoretical and Computational Chemistry_, Recent Developments and Applications of Modern Density Functional Theory, Vol. 4, edited by J. M. Seminario (Elsevier, 1996) pp. 327–357.
  * Iikura  _et al._ (2001) H. Iikura, T. Tsuneda, T. Yanai, and K. Hirao, J. Chem. Phys. 115, 3540 (2001), publisher: American Institute of Physics.
  * Yanai  _et al._ (2004) T. Yanai, D. P. Tew, and N. C. Handy, Chemical Physics Letters 393, 51 (2004).
  * Vydrov and Scuseria (2006) O. A. Vydrov and G. E. Scuseria, J. Chem. Phys. 125, 234109 (2006), publisher: American Institute of Physics.
  * Chai and Head-Gordon (2008) J.-D. Chai and M. Head-Gordon, Phys. Chem. Chem. Phys. 10, 6615 (2008), publisher: The Royal Society of Chemistry.
  * Baer  _et al._ (2010) R. Baer, E. Livshits, and U. Salzner, Annual Review of Physical Chemistry 61, 85 (2010), _eprint: https://doi.org/10.1146/annurev.physchem.012809.103321.
  * Wing  _et al._ (2021) D. Wing, G. Ohad, J. B. Haber, M. R. Filip, S. E. Gant, J. B. Neaton, and L. Kronik, PNAS 118, 10.1073/pnas.2104556118 (2021), publisher: National Academy of Sciences Section: Physical Sciences.
  * Perdew and Zunger (1981) J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981), publisher: American Physical Society.
  * Mori-Sánchez  _et al._ (2006b) P. Mori-Sánchez, A. J. Cohen, and W. Yang, J. Chem. Phys. 124, 091102 (2006b), publisher: American Institute of Physics.
  * Perdew  _et al._ (2008) J. P. Perdew, V. N. Staroverov, J. Tao, and G. E. Scuseria, Phys. Rev. A 78, 052513 (2008), publisher: American Physical Society.
  * Schmidt  _et al._ (2014) T. Schmidt, E. Kraisler, L. Kronik, and S. Kümmel, Physical Chemistry Chemical Physics 16, 14357 (2014), publisher: Royal Society of Chemistry.
  * Pederson  _et al._ (2014) M. R. Pederson, A. Ruzsinszky, and J. P. Perdew, J. Chem. Phys. 140, 121103 (2014), publisher: American Institute of Physics.
  * Schmidt and Kümmel (2016) T. Schmidt and S. Kümmel, Phys. Rev. B 93, 165120 (2016), publisher: American Physical Society.
  * Yang  _et al._ (2017) Z.-h. Yang, M. R. Pederson, and J. P. Perdew, Phys. Rev. A 95, 052505 (2017), publisher: American Physical Society.
  * Borghi  _et al._ (2014) G. Borghi, A. Ferretti, N. L. Nguyen, I. Dabo, and N. Marzari, Phys. Rev. B 90, 075135 (2014).
  * Colonna  _et al._ (2019) N. Colonna, N. L. Nguyen, A. Ferretti, and N. Marzari, J. Chem. Theory Comput. 15, 1905 (2019), publisher: American Chemical Society.
  * Anisimov and Kozhevnikov (2005) V. I. Anisimov and A. V. Kozhevnikov, Phys. Rev. B 72, 075125 (2005), publisher: American Physical Society.
  * Ma and Wang (2016) J. Ma and L.-W. Wang, Sci Rep 6, 24924 (2016).
  * Li  _et al._ (2018) C. Li, X. Zheng, N. Q. Su, and W. Yang, National Science Review 5, 203 (2018).
  * Su  _et al._ (2020) N. Q. Su, A. Mahler, and W. Yang, J. Phys. Chem. Lett. 11, 1528 (2020).
  * Mei  _et al._ (2020) Y. Mei, Z. Chen, and W. Yang, The Journal of Physical Chemistry Letters 10.1021/acs.jpclett.0c03133 (2020), publisher: American Chemical Society.
  * Mei  _et al._ (2021a) Y. Mei, N. Yang, and W. Yang, J. Chem. Phys. 154, 054302 (2021a), publisher: American Institute of Physics.
  * Bloch (1929) F. Bloch, Z. Physik 52, 555 (1929).
  * Ashcroft and Mermin (2003) N. Ashcroft and D. Mermin,  _Solid State Physics_ (Thomson Press, 2003).
  * Monkhorst and Pack (1976) H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
  * Wannier (1937) G. Wannier, Physical Review 52, 191 (1937).
  * Cloizeaux (1963) J. D. Cloizeaux, Phys. Rev. 129, 554 (1963), publisher: American Physical Society.
  * Marzari and Vanderbilt (1997) N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 (1997), publisher: American Physical Society.
  * Foster and Boys (1960) J. M. Foster and S. F. Boys, Rev. Mod. Phys. 32, 300 (1960), publisher: American Physical Society.
  * Gygi  _et al._ (2003) F. Gygi, J.-L. Fattebert, and E. Schwegler, Computer Physics Communications 155, 1 (2003).
  * Giustino and Pasquarello (2006) F. Giustino and A. Pasquarello, Phys. Rev. Lett. 96, 216403 (2006).
  * Mahler  _et al._ (2022) A. Mahler, J. Z. Williams, N. Q. Su, and W. Yang, arXiv:2201.07751 [cond-mat] (2022), arXiv: 2201.07751.
  * Mei  _et al._ (2022) Y. Mei, J. Yu, Z. Chen, N. Q. Su, and W. Yang, J. Chem. Theory Comput. 10.1021/acs.jctc.1c01058 (2022).
  * Zheng  _et al._ (2011) X. Zheng, A. J. Cohen, P. Mori-Sánchez, X. Hu, and W. Yang, Phys. Rev. Lett. 107, 026403 (2011), publisher: American Physical Society.
  * Hait and Head-Gordon (2018) D. Hait and M. Head-Gordon, J. Phys. Chem. Lett. 9, 6280 (2018), publisher: American Chemical Society.
  * Mei  _et al._ (2021b) Y. Mei, Z. Chen, and W. Yang, arXiv:2106.10358 [physics, physics:quant-ph] (2021b), arXiv: 2106.10358.
  * Dirac (1930) P. a. M. Dirac, Mathematical Proceedings of the Cambridge Philosophical Society 26, 376 (1930), publisher: Cambridge University Press.
  * Slater and Phillips (1974) J. C. Slater and J. C. Phillips, Physics Today 27, 49 (1974), publisher: American Institute of Physics.
  * Slater (1974) J. C. Slater,  _Quantum Theory of Molecules and Solids_ , Vol. IV (McGrow-Hill, 1974).
  * Souza  _et al._ (2001) I. Souza, N. Marzari, and D. Vanderbilt, Phys. Rev. B 65, 035109 (2001), publisher: American Physical Society.
  * Makov and Payne (1995) G. Makov and M. C. Payne, Phys. Rev. B 51, 4014 (1995).
  * Kantorovich (1999) L. N. Kantorovich, Phys. Rev. B 60, 15476 (1999), publisher: American Physical Society.
  * Dabo  _et al._ (2008) I. Dabo, B. Kozinsky, N. E. Singh-Miller, and N. Marzari, Phys. Rev. B 77, 115139 (2008), arXiv: 0709.4647.
  * Li and Dabo (2011) Y. Li and I. Dabo, Phys. Rev. B 84, 155127 (2011), publisher: American Physical Society.
  * Onida  _et al._ (1995) G. Onida, L. Reining, R. W. Godby, R. Del Sole, and W. Andreoni, Phys. Rev. Lett. 75, 818 (1995), publisher: American Physical Society.
  * Jarvis  _et al._ (1997) M. R. Jarvis, I. D. White, R. W. Godby, and M. C. Payne, Phys. Rev. B 56, 14972 (1997), publisher: American Physical Society.
  * Rozzi  _et al._ (2006) C. A. Rozzi, D. Varsano, A. Marini, E. K. U. Gross, and A. Rubio, Phys. Rev. B 73, 205119 (2006).
  * Hedin (1965) L. Hedin, Physical Review series I 139, 796 (1965), number: 3A Publisher: American Institute of Physics (AIP).
  * Martin  _et al._ (2016) R. Martin, L. Reining, and D. Ceperley,  _Interacting Electrons: Theory and Computational Approaches_ (Cambridge University Press, 2016).
  * Zaghloul and Ali (2012) M. R. Zaghloul and A. N. Ali, ACM Trans. Math. Softw. 38, 15:1 (2012).
  * Zaghloul (2016) M. R. Zaghloul, ACM Trans. Math. Softw. 42, 26:1 (2016).
  * Zheng  _et al._ (2013) X. Zheng, T. Zhou, and W. Yang, J. Chem. Phys. 138, 174105 (2013).
  * Zhang  _et al._ (2015) D. Zhang, X. Zheng, C. Li, and W. Yang, J. Chem. Phys. 142, 154113 (2015).
  * Colonna  _et al._ (2018) N. Colonna, N. L. Nguyen, A. Ferretti, and N. Marzari, J. Chem. Theory Comput. 14, 2549 (2018).
  * Nguyen  _et al._ (2018) N. L. Nguyen, N. Colonna, A. Ferretti, and N. Marzari, Phys. Rev. X 8, 021051 (2018).
  * Colonna  _et al._ (2022) N. Colonna, R. De Gennaro, E. Linscott, and N. Marzari, arXiv (2022), arXiv:2202.08155 [cond-mat, physics:physics] .
  * Yang  _et al._ (2012b) W. Yang, A. J. Cohen, F. De Proft, and P. Geerlings, J. Chem. Phys. 136, 144110 (2012b).
  * Hamann (2013) D. R. Hamann, Phys. Rev. B 88, 085117 (2013).
  * van Setten  _et al._ (2018) M. J. van Setten, M. Giantomassi, E. Bousquet, M. J. Verstraete, D. R. Hamann, X. Gonze, and G. M. Rignanese, Computer Physics Communications 226, 39 (2018).
  * Giannozzi  _et al._ (2009) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. d. Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, J. Phys.: Condens. Matter 21, 395502 (2009), publisher: IOP Publishing.
  * Giannozzi  _et al._ (2017) P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. d. Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. Otero-de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, J. Phys.: Condens. Matter 29, 465901 (2017), publisher: IOP Publishing.
  * Mostofi  _et al._ (2014) A. A. Mostofi, J. R. Yates, G. Pizzi, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, Computer Physics Communications 185, 2309 (2014).
  * Pizzi  _et al._ (2020) G. Pizzi, V. Vitale, R. Arita, S. Blügel, F. Freimuth, G. Géranton, M. Gibertini, D. Gresch, C. Johnson, T. Koretsune, J. Ibañez-Azpiroz, H. Lee, J.-M. Lihm, D. Marchand, A. Marrazzo, Y. Mokrousov, J. I. Mustafa, Y. Nohara, Y. Nomura, L. Paulatto, S. Poncé, T. Ponweiser, J. Qiao, F. Thöle, S. S. Tsirkin, M. Wierzbowska, N. Marzari, D. Vanderbilt, I. Souza, A. A. Mostofi, and J. R. Yates, J. Phys.: Condens. Matter 32, 165902 (2020).
  * Heyd  _et al._ (2005) J. Heyd, J. E. Peralta, G. E. Scuseria, and R. L. Martin, J. Chem. Phys. 123, 174101 (2005), publisher: American Institute of Physics.
  * Su  _et al._ (2018) N. Q. Su, C. Li, and W. Yang, PNAS 115, 9678 (2018).
  * Anisimov  _et al._ (1991) V. I. Anisimov, J. Zaanen, and O. K. Andersen, Phys. Rev. B 44, 943 (1991), publisher: American Physical Society.
  * Cococcioni and de Gironcoli (2005) M. Cococcioni and S. de Gironcoli, Phys. Rev. B 71, 035105 (2005).
  * O’Regan (2012) D. D. O’Regan, in  _Optimised Projections for the Ab Initio Simulation of Large and Strongly Correlated Systems_, edited by D. D. O’Regan (Springer, Berlin, Heidelberg, 2012) pp. 89–123.
  * Campo and Cococcioni (2010) V. L. J. Campo and M. Cococcioni, J. Phys.: Condens. Matter 22, 055602 (2010).
  * Himmetoglu  _et al._ (2014) B. Himmetoglu, A. Floris, S. de Gironcoli, and M. Cococcioni, Int. J. Quantum Chem. 114, 14 (2014).
  * Bajaj  _et al._ (2017) A. Bajaj, J. P. Janet, and H. J. Kulik, J. Chem. Phys. 147, 191101 (2017).
  * Bajaj  _et al._ (2019) A. Bajaj, F. Liu, and H. J. Kulik, J. Chem. Phys. 150, 154115 (2019).
  * Dabo  _et al._ (2009) I. Dabo, M. Cococcioni, and N. Marzari, arXiv:0901.2637 [cond-mat] (2009), arXiv:0901.2637 [cond-mat] .
  * Dabo  _et al._ (2010) I. Dabo, A. Ferretti, N. Poilvert, Y. Li, N. Marzari, and M. Cococcioni, Phys. Rev. B 82, 115121 (2010).
  * Borghi  _et al._ (2015) G. Borghi, C.-H. Park, N. L. Nguyen, A. Ferretti, and N. Marzari, Phys. Rev. B 91, 155112 (2015).
  * Kao  _et al._ (2017) D.-y. Kao, K. Withanage, T. Hahn, J. Batool, J. Kortus, and K. Jackson, J. Chem. Phys. 147, 164107 (2017), publisher: American Institute of Physics.
  * Jackson  _et al._ (2019) K. A. Jackson, J. E. Peralta, R. P. Joshi, K. P. Withanage, K. Trepte, K. Sharkas, and A. I. Johnson, J. Phys.: Conf. Ser. 1290, 012002 (2019).