# Automated mixing of maximally localized Wannier functions into target
manifolds

Junfeng Qiao<EMAIL_ADDRESS>Theory and Simulations of Materials
(THEOS), and National Centre for Computational Design and Discovery of Novel
Materials (MARVEL), École Polytechnique Fédérale de Lausanne, 1015 Lausanne,
Switzerland  Giovanni Pizzi  Theory and Simulations of Materials (THEOS), and
National Centre for Computational Design and Discovery of Novel Materials
(MARVEL), École Polytechnique Fédérale de Lausanne, 1015 Lausanne, Switzerland
Laboratory for Materials Simulations (LMS), Paul Scherrer Institut (PSI),
CH-5232 Villigen PSI, Switzerland  Nicola Marzari  Theory and Simulations of
Materials (THEOS), and National Centre for Computational Design and Discovery
of Novel Materials (MARVEL), École Polytechnique Fédérale de Lausanne, 1015
Lausanne, Switzerland  Laboratory for Materials Simulations (LMS), Paul
Scherrer Institut (PSI), CH-5232 Villigen PSI, Switzerland

###### Abstract

Maximally localized Wannier functions (MLWFs) are widely used to construct
first-principles tight-binding models that accurately reproduce the electronic
structure of materials. Recently, robust and automated approaches to generate
these MLWFs have emerged, leading to natural sets of atomic-like orbitals that
describe both the occupied states and the lowest-lying unoccupied ones (when
the latter can be meaningfully described by bonding/anti-bonding combinations
of localized orbitals). For many applications, it is important to instead have
MLWFs that describe only certain target manifolds separated in energy between
them—the occupied states, the empty states, or certain groups of bands. Here,
we start from the full set of MLWFs describing simultaneously all the target
manifolds, and then mix them using a combination of parallel transport and
maximal localization to construct orthogonal sets of MLWFs that fully and only
span the desired target submanifolds. The algorithm is simple and robust, and
it is applied to some paradigmatic but non-trivial cases (the valence and
conduction bands of silicon, the top valence band of MoS2, the $3d$ and
$t_{2g}$/$e_{g}$ bands of SrVO3) and to a mid-throughput study of 77
insulators.

##  I Introduction

Maximally localized Wannier functions (MLWFs) [1, 2, 3, 4] are accurate
reduced-order models [5] for the electronic structures of periodic crystals.
The generation of MLWFs from Bloch wavefunctions typically requires a choice
of initial projectors, which are often conjectured from chemical intuition
with trial and error. For metals, or when considering both the valence and
conduction bands (VCB) of insulators and semiconductors, one typically deals
with bands that overlap (i.e., are entangled) [2] with higher-energy bands. In
such cases, since the low-energy electronic structure can often be well
described by a tight-binding model of atomic-like orbitals, the initial
projectors are usually chosen from the hydrogenic $s,p,d,f$ orbitals. However,
when it comes to the case of identifying the optimal projectors for valence
bands (VB) alone, or especially conduction bands (CB) which are mixed with
higher-energy bands, it might become difficult to find good initial guesses.
Indeed, the VB/CB often consists of bonding/anti-bonding orbitals, or
combination of atomic orbitals which are more challenging to guess or to
describe, unless the crystal offers a very simple chemical picture. Achieving
separate Wannierization of target manifolds is also advantageous for many
applications. Some physical properties (such as the electric polarization)
depend only on the Wannier functions (WFs) of the occupied manifold (sum of
Wannier centers of all the valence WFs). Moreover, using dedicated MLWFs means
that one can obtain smaller tight-binding models that are thus more efficient
when computing, e.g., transport properties of large systems. Koopmans spectral
functionals also require separate occupied and unoccupied manifolds [6]. Last,
low-energy models, such as those used in correlated-electrons calculations [7,
8, 9, 10], require a description of the correlated manifold.

Several approaches have been developed in the past few years to simplify the
construction of MLWFs. The selected columns of the density matrix (SCDM)
algorithm [11] uses QR decomposition with column pivoting on the density
matrix to automatically generate initial projection orbitals, and a sensible
choice of the density matrix can be obtained from the projectability of Bloch
states onto pseudo-atomic orbitals from pseudopotentials [12]. The optimal
projection functions method [13] starts with a larger manifold and generates
the MLWFs of the valence manifold by a single rotation matrix, which is
computed by a product of a series of Givens rotations. The dually localized
Wannier functions method [14] adds an additional term to spread functional, to
localize the WFs in both space and energy, achieving a separation of VB and
CB.

Here, we propose a different approach to automatically mix optimal MLWFs
spanning valence and conduction into several submanifolds, provided that these
submanifolds are gapped in their energy spectrum. This naturally applies to
the case of separate Wannierizations of valence and conduction manifolds, but
more generally extends to arbitrary groups of bands separated in energy. We
start from the Wannierization of a larger manifold (e.g., the VCB manifold),
that we do not discuss here since robust methods already exist: in addition to
hydrogenic $s,p,d,f$ initial projections, partly occupied WF method [15, 16],
the fully automated SCDM method [11, 12], or the projectability-disentangled
Wannier function (PDWF) that we recently introduced are available (in
particular, the latter appears as a very general and remarkably robust
approach allowing, e.g., to construct $\sim 1.3$ million PDWFs for $\sim 22$
thousands materials [17]). Once these MLWFs are obtained, we then diagonalize
the Wannier Hamiltonian at every $k$-point and partition the states into
submanifolds (e.g., valence, conduction): they are grouped together if they
fall inside the desired energy interval. Next, we fix the gauge randomness of
the submanifolds using parallel transport [18]. Finally, the MLWFs for each
submanifold are generated by maximally localizing their spread functionals,
independently. Since the submanifolds are already isolated in energy (i.e., no
disentanglement is needed), and parallel transport provides a continuous
gauge, the final maximal localization converges effortlessly. In the case of
separating VB and CB, the final two groups of MLWFs span the fully occupied
valence and the fully unoccupied conduction manifolds, and their shapes
closely resemble bonding and anti-bonding orbitals, respectively. Compared
with SCDM, the present method works fully in reciprocal space, reducing the
computer memory requirements and also being computationally faster. Compared
with the optimal projection functions method [13] or the dually localized
Wannier functions method [14], we do not change the spread functional but use
the original one in Marzari and Vanderbilt [1], thus the resulting WFs are
maximally-localized in their original definition; moreover, the parallel
transport step is non-iterative and always quickly provides a good starting
point for the final maximally-localization step, avoiding potential
convergence issues that might occur in an iterative method.

In the following, we first discuss and validate the present method, which we
name manifold-remixed Wannier function (MRWF), on the VCB of a 3D material
(silicon), the VCB of a 2D material (MoS2), the top VB of MoS2, and the $3d$
manifold of SrVO3. We also discuss the bonding/anti-bonding character of the
resulting MLWFs, as well as band interpolation accuracy. To analyze statistics
of band interpolation quality and demonstrate the robustness of the present
approach, we Wannierize the VB and CB of a diverse set of 77 insulators, with
the number of atoms between 1 and 45.

##  II Results

###  II.1 The manifold separation algorithm

While obtaining the starting WFs is not the focus of this paper, we remind
here that the standard Wannierization algorithm [1, 2] requires initial
projection orbitals $\ket{g_{n}}$ to guide the spread minimization and find
the most meaningful minimum and the related unitary transformation matrices
$U_{\mathbf{k}}$ at each $k$-point $\mathbf{k}$. The projectors $\ket{g_{n}}$
are used to rotate the original Bloch wavefunctions $\ket{\psi_{m\mathbf{k}}}$
into

$\ket{\tilde{\psi}_{n\mathbf{k}}}=\sum_{m=1}^{M}\ket{\psi_{m\mathbf{k}}}\braket{\psi_{m\mathbf{k}}}{g_{n}},$
(1)

where $n$ and $m$ are the indices for WFs and Bloch bands, respectively; $M$
is the total number of Bloch bands; and $\mathbf{k}$ is the Bloch quasi-
momentum. Note that $\ket{\tilde{\psi}}$ are independent of any arbitrary
rotation gauge for the $\ket{\psi_{n\mathbf{k}}}$. For metals or for VCB of
insulators, one typically starts with hydrogenic $s,p,d,f$ orbitals [1] as the
initial guesses for all the corresponding valence electrons. Then, the MLWFs
can be generated using either the standard disentanglement [2] and maximal
localization algorithms [1] or minimizing directly the total spread, such as
the partly occupied WF method [15] or a variational formulation [19]. Instead
of hydrogenic orbitals, one can use SCDM [11] or the recently introduced
projectability disentanglement [17] for a fully automated Wannierization.
Irrespective of the approach taken to obtain MLWFs describing the VCB, these
MLWFs will be the starting point of the present algorithm, with the next step
to separate e.g. the VB and CB manifolds from the disentangled MLWFs that span
both simultaneously. Note that while in the following we use the separation of
VB and CB as an example to illustrate the method for clarity and simplicity,
the present approach is not limited to the case of two submanifolds, but can
be applied to any groups of bands separated in energy.

Since the disentanglement procedure aims at obtaining the lowest-possible
spreads, it typically achieves this goal by mixing states originating from all
the submanifolds (e.g., VB and CB) of interest. To decompose the manifold into
two orthogonal submanifolds, we diagonalize the Wannier-gauge Hamiltonian
${H}^{W}_{\mathbf{k}}$ (the superscript $W$ indicates the Wannier gauge),

$H^{W}_{\mathbf{k}}=V_{\mathbf{k}}\mathcal{E}_{\mathbf{k}}V^{*}_{\mathbf{k}},$
(2)

where $\mathcal{E}_{\mathbf{k}}$ and $V_{\mathbf{k}}$ are the eigenvalues and
the eigenvectors, respectively; $*$ denotes conjugate transpose. Usually the
eigenvalues and eigenvectors returned from linear algebra computer programs
are already sorted in ascending order of eigenvalues; if not, we sort them in
ascending order, so that the matrices are partitioned into two blocks,

$\displaystyle\mathcal{E}_{\mathbf{k}}$
$\displaystyle=\begin{bmatrix}\mathcal{E}^{1}_{\mathbf{k}}&0\\\
0&\mathcal{E}^{2}_{\mathbf{k}}\end{bmatrix},$ (3) $\displaystyle
V_{\mathbf{k}}$
$\displaystyle=\begin{bmatrix}V^{1}_{\mathbf{k}}&V^{2}_{\mathbf{k}}\end{bmatrix},$

where $V^{1}_{\mathbf{k}}\in\mathbb{C}^{N\times P}$
($V^{2}_{\mathbf{k}}\in\mathbb{C}^{N\times Q}$) corresponds to states whose
eigenvalues $\mathcal{E}^{1}_{\mathbf{k}}$ ($\mathcal{E}^{2}_{\mathbf{k}}$)
are below (above) the band gap, and $0$ represents a zero matrix. Here, $N$ is
the number of WFs of the VCB manifold, $P$ and $Q$ are the number of WFs in
the valence (below band gap) and the conduction (above band gap) submanifolds,
respectively, such that $N=P+Q$. Next, all the Wannier-gauge operators are
rotated according to $V^{1}_{\mathbf{k}}$ for the valence submanifold: for
instance, the overlap matrices $M^{W}_{\mathbf{k},\mathbf{b}}$ (for computing
the spread functional) is rotated by

$M_{\mathbf{k},\mathbf{b}}^{W,1}=V^{1*}_{\mathbf{k}}M^{W}_{\mathbf{k},\mathbf{b}}V^{1}_{\mathbf{k+b}},$
(4)

where

$\displaystyle M^{W}_{\mathbf{k},\mathbf{b}}$
$\displaystyle=U^{*}_{\mathbf{k}}M_{\mathbf{k},\mathbf{b}}U_{\mathbf{k}+\mathbf{b}},$
(5) $\displaystyle M_{\mathbf{k}\mathbf{b}}$
$\displaystyle=\braket{u_{m,\mathbf{k}}}{u_{n,\mathbf{k+b}}},$

$U_{\mathbf{k}}$ are the unitary transformations from the VCB manifold
Wannierization, and $\ket{u_{m,\mathbf{k}}}$ is the periodic part of Bloch
wavefunction
$\ket{\psi_{m,\mathbf{k}}}=e^{i\mathbf{k}\mathbf{r}}\ket{u_{m,\mathbf{k}}}$.
For more details on the notations of $M_{\mathbf{k},\mathbf{b}}$ and
$\mathbf{b}$-vectors, see Marzari and Vanderbilt [1]. Consistently, the
$\mathcal{E}^{1}_{\mathbf{k}}$ is used as the new eigenvalues. Now the problem
is reformulated into a Wannierization of an isolated submanifold with $P$ WFs
for VB. Similarly, the conduction manifold operators are rotated by
$V^{2}_{\mathbf{k}}$, and an analogous Wannierization of an isolated
submanifold with $Q$ WFs. Indeed, the first-step Wannierization of VCB has
already disentangled the MLWFs from all the remaining higher-energy bands, so
that these MLWFs span all the target submanifolds we are interested in. To
achieve our goal, we are left with Wannierizations of two isolated
submanifolds and thus the subsequent steps do not need any disentanglement.
Such a two-step procedure makes the whole algorithm quite robust, especially
when Wannierizing the CB, for which it is difficult to provide good initial
projections of the corresponding anti-bonding orbitals.

The remaining difficulty of the Wannierization of the two isolated
submanifolds is caused by the diagonalization in Eq. 2. Indeed, since the
Hamiltonians $H^{W}_{\mathbf{k}}$ are independently diagonalized at each
$\mathbf{k}$, the resulting eigenvectors will have different gauges at
different $k$-points, requiring additional Wannierizations in each
submanifold. Since these Wannierizations are carried out on submanifolds that
have isolated bands, the minimization algorithm is typically more robust to
the choice of initial projections compared to the case of disentanglement. One
could simply resort to random Gaussian initial projections followed by maximal
localization to reach the MLWFs for the two submanifolds, respectively; or
even brute-force maximal localization starting from the random gauge after the
Wannier Hamiltonian diagonalization. However, a direct maximal localization
starting from a random gauge is not robust—we observe that, in many cases, the
maximal localization fails due to zeros in the diagonal of the overlap
matrices $M_{\mathbf{k},\mathbf{b}}$; and, even if it converges, it displays
the same issues of random Gaussian projections: a large number of iterations,
and oscillatory evolution of spread and the sum of MLWF centers during the
minimization process (see LABEL:supp-fig:si2_val_4gauge in supplementary
information (SI) for discussions on the convergences of these choices).
Moreover, when the number of $k$-points $N_{\mathbf{k}}$ is large, the maximal
localization is much harder to converge. A better solution is finding good
starting gauge for the two submanifolds in an automated fashion.

To tackle this challenge, we adopt the parallel transport algorithm [18] to
construct smooth gauges for the two submanifolds. For an isolated manifold,
the existence of a smooth gauge is determined by its topological obstructions,
which are characterized by the Chern numbers (one in 2D and three in 3D). If
the Chern numbers are 0 (as it is the case for systems with time-reversal
symmetry), then it can be proven [20] that it is possible to construct a
continuous gauge explicitly by the following procedure [18]: (a) Suppose
$k_{i}\in[0,1]$ (in fractional coordinates) for $i=x,y,z$: propagate (using
singular value decomposition of overlap matrices
$M_{\mathbf{k},\mathbf{b}}^{W}$ to maximally align the gauge between
neighboring $k$-points) the Bloch wavefunctions $\ket{u_{n\mathbf{0}}}$ at
$\Gamma$ along $k_{x}$ from $k$-point $\mathbf{0}=(0,0,0)$ to
$\mathbf{1}=(1,0,0)$, to construct a continuous gauge across these $k$-points.
The new gauge is not necessarily quasi-periodic, i.e., satisfying the Bloch
theorem imposed on $\ket{u_{n\mathbf{k}}}$ by
$\ket{u_{n\mathbf{k+K}}}=\tau_{\mathbf{K}}\ket{u_{n\mathbf{k}}}$ where
$\tau_{\mathbf{K}}=e^{-i\mathbf{K}\cdot\mathbf{r}}$ is the translation
operator in reciprocal space, and $\mathbf{K}$ is a reciprocal lattice vector.
In general, instead, the two states are related by
$\ket{u_{n\mathbf{1}}}=(\tau_{\mathbf{1}}\ket{u_{n\mathbf{0}}})V_{\mathrm{obs}}$,
where this expression defines the obstruction matrix $V_{\mathrm{obs}}$
quantifying the misalignment of the propagated gauge and the gauge required by
Bloch theorem at $(1,0,0)$. To fulfill the quasi-periodic boundary condition,
we can therefore multiply each $\ket{u_{n\mathbf{k}}}$ by $e^{-k_{x}L}$ (note
$k_{x}\in[0,1]$) where $V_{\mathrm{obs}}=\exp(L)$: in this way, we obtain a
continuous gauge that also satisfies Bloch theorem, i.e., the new obstruction
matrix in this modified gauge is the identity matrix. (b) For each $k_{x}$,
propagate along $k_{y}$ from $(k_{x},0,0)$ to $(k_{x},1,0)$. Now we obtain a
series of obstruction matrices $V_{\mathrm{obs}}(k_{x})$ along $k_{x}$. If the
winding number [21, 22, 18] of the determinants of $V_{\mathrm{obs}}(k_{x})$
vanishes (i.e., the Chern number is 0), then there is a continuous function
that maps $V_{\mathrm{obs}}(k_{x})$ to identity [18]. We then multiply the
gauge by this mapping, so that the new gauge satisfies the quasi-periodic
boundary condition in the $k_{x}-k_{y}$ plane. Ref. [18] explicitly constructs
the continuous mapping by their “column interpolation” method for the Kane-
Mele model, which is a 2D fermionic time-reversal-symmetric model (i.e.,
having a vanishing Chern number) but can present a non-zero $\mathbb{Z}_{2}$
number; as a comparison, previous methods had difficulties in handling
$\mathbb{Z}_{2}$ systems [23, 24, 18], sometimes requiring model-specific
information [25, 26]. (c) For each ($k_{x},k_{y}$), propagate along $k_{z}$
from $(k_{x},k_{y},0)$ to $(k_{x},k_{y},1)$. Now the obstruction matrices
$V_{\mathrm{obs}}(k_{x},k_{y})$ depend both on $k_{x}$ and $k_{y}$. Similar to
point (b), if the two winding numbers of the determinants of
$V_{\mathrm{obs}}(k_{x},0)$ and $V_{\mathrm{obs}}(0,k_{y})$ vanish, then there
is a continuous function that maps $V_{\mathrm{obs}}(k_{x},k_{y})$ to
identity. We then multiply the gauge with this mapping and obtain the final
gauge satisfying the quasi-periodic boundary condition in 3D. Ref. [18]
demonstrates this constructive algorithm to obtain a continuous gauge for a 3D
system (silicon). The results also show that the continuous gauge can be
further smoothened by the standard maximal localization procedure [1] to
construct MLWFs. We stress that the algorithm is non-iterative and fast, thus
solving the problem of finding good initial WFs for isolated manifolds in an
efficient and robust way.

As shown in LABEL:supp-fig:si2_val_4gauge, parallel transport generates a much
better starting point than random Gaussian projections or random gauges: the
convergence of maximal localization is much faster, and the spread and the sum
of MLWF centers smoothly evolve during minimization. We note that since the
propagation of gauge requires overlap matrices between a particular set of
nearest-neighboring $k$-points, in the SI LABEL:supp-sec:bvec_pt we present a
procedure so that parallel transport can be applied to any arbitrary crystal
structure.

In summary, the sequential parallel transports move the obstructions to the
Brillouin zone edges, and the “column interpolation” method fixes the quasi-
periodicity. Our tests on a set of 77 insulators (see discussion later in
Section II.5) show that this algorithm is able to construct a good initial
gauge, and maximal localization is able to construct MLWFs without issue.

We now mention that since we propagate the gauge starting from the first
$k$-point $(0,0,0)$, there is still one gauge arbitrariness at this $\Gamma$
point. Here, we suggest to insert an additional step that first minimizes the
spread functional w.r.t. a single rotation matrix $W$ for the first $k$-point,
before performing the standard maximal localization w.r.t. all $k$-points to
obtain MLWFs. Indeed, thanks to the small size of $W$, this first preliminary
step is computationally efficient, and can help in further improving the
overall robustness of the full algorithm that we are presenting here. To
achieve this, we derive the expression of the gradient of the spread $\Omega$
w.r.t. the rotation matrix $W$ in SI LABEL:supp-sec:opt_rot. We then use this
gradient with a manifold optimization algorithm [27] to minimize $\Omega$
w.r.t. $W$, where $W$ is constrained on the unitary matrix manifold
$\left\\{W\in\mathbb{C}^{K\times K}|W^{*}W=I\right\\}$, where $K=P$ for the
valence manifold, or $K=Q$ for the conduction manifold. This minimization
provides us with a single rotation matrix $W$ that further improves the
localization, while still preserving the parallel transport gauge: we stress
that, in addition to increasing the robustness of the algorithm as mentioned
earlier, this additional step can thus be beneficial for cases where the
parallel transport gauge is implicitly assumed during the derivation of
equations (for instance, Wannier interpolation of Berry curvature [28], or
Wannier interpolation of nonlinear optical responses [29, 30]).

After the parallel transport and the single rotation, the resulting WFs are
close to the ideal MLWFs. However, since parallel transport only generates a
continuous quasi-periodic gauge, it typically does not provide the smallest
possible spread. It is therefore helpful to perform a final smoothing of the
gauge [18] by running a final maximal localization step (see examples in
Sections II.2 and II.3). This can be achieved using either the original
Marzari-Vanderbilt localization [1] or a matrix manifold optimization w.r.t.
gauge matrices at all the $k$-points, i.e., optimization on a product manifold
of a series of unitary matrices
$\prod_{\mathbf{k}}\left\\{U_{\mathbf{k}}\in\mathbb{C}^{K\times
K}|U^{*}_{\mathbf{k}}U_{\mathbf{k}}=I\right\\}$, where $K=P$ for valence and
$K=Q$ for conduction manifolds. As already mentioned, the multi-step procedure
that we propose here aims at making the whole algorithm more robust, since
every step produce a better starting point for the final iterative
localization algorithm.

In summary, we start from an initial manifold that has been already singled
out from the remaining high-energy states using standard procedure such as
disentanglement and maximal localization (e.g., very accurately using
projectability disentanglement [17] to extract as much as possible the bonding
and anti-bonding characters from all the bands). The subsequent
diagonalizations of Wannier-gauge Hamiltonians separate the manifold into
(two) orthogonal submanifolds (for VB & CB, respectively). The (two) parallel-
transport steps (for the relevant submanifolds) construct continuous gauges,
fixing the randomness caused by the independent Hamiltonian diagonalization at
each $k$-point. The rotation w.r.t. a single unitary matrix removes the gauge
arbitrariness of parallel transport at the first $k$-point. The final maximal
localizations ultimately smoothen the gauge, leading to two sets of MLWFs,
each of which spans the submanifold for VB or CB. In SI LABEL:supp-
sec:MRWF_gauge, we prove that the final gauge transformation has block
diagonal structure, i.e., the MRWFs are transformed according to

$U(\mathbf{k})=\begin{bmatrix}U_{\mathrm{VB}}(\mathbf{k})&0\\\
0&U_{\mathrm{CB}}(\mathbf{k})\end{bmatrix},$ (6)

where $U_{\mathrm{VB}}(\mathbf{k})$ and $U_{\mathrm{CB}}(\mathbf{k})$ are
unitary matrices for VB and semi-unitary matrices for CB, respectively.

###  II.2 Silicon

To test the validity of the present method, we first disentangle and maximally
localize the VCB of silicon into 8 WFs, using the standard hydrogenic $s$ and
$p$ projections with energy window disentanglement (we use hydrogenic
projectors and energy window disentanglement here to demonstrate that the
present approach works well as long as the entire VCB are accurately
described; one can also use PDWF to construct MLWFs spanning the entire VCB).
The resulting WFs, two of which have $s$-character and six of which have
$p$-character, are shown in the VCB column of Fig. 1.

Figure 1: WF shapes and band interpolation of silicon. 1 WF shape of various
Wannierization methods. In the header, the left, center, and right columns
correspond to the valence plus conduction bands (VCB), valence bands (VB), and
conduction bands (CB), respectively. The PT column shows WFs after running
parallel transport (PT); The PT+SR column shows WFs after running PT and
single rotation (SR) of $W$ matrix; the PT(+SR)+ML column shows WFs after
running PT (and optionally SR) and maximal localization (ML), since PT+SR+ML
and PT+ML fall to the same minimum. The numbers under each WF are the spread
and its multiplicity: e.g., 1.668 (2) means there are two WFs having spread
$1.668\text{\,}\mathrm{\text{Å}}^{2}$. The last row $\sum\Omega_{\mathrm{WF}}$
is the total spread of each calculation. The blue spheres are the silicon
atoms. 1 Band structure comparison of DFT, VCB, VB, and CB. The horizontal
dashed blue line is the upper limit of the frozen window for disentanglement
(DIS). The VCB are computed by Wannier interpolation from DIS+ML; the VB and
CB are computed by Wannier interpolation from PT+ML.

For the valence manifold, after running the Hamiltonian diagonalization and
parallel transport, we obtain four WFs with spreads around
$4\text{\,}\mathrm{\text{Å}}^{2}$, but their shapes do not have clear physical
meaning, since the gauge of the first $k$-point is still arbitrary (see PT
column inside the VB column of Fig. 1); after running the single rotation of
the $W$ matrix, their spreads are further minimized to around
$3\text{\,}\mathrm{\text{Å}}^{2}$, and the shapes now resemble the bonding
orbitals between neighboring silicon atoms (PT+SR column inside the VB column
of Fig. 1); a final maximal localization further reduces the spreads to around
$2\text{\,}\mathrm{\text{Å}}^{2}$ and the four spreads end up becoming
identical (PT(+SR)+ML column inside the VB column of Fig. 1) thus respecting
the symmetry of the full system. For the conduction manifold (CB columns of
Fig. 1), again the shapes of WFs after parallel transport have no clear
meaning. However, even after the single rotation, the shapes of WFs still do
not resemble the expected anti-bonding orbitals, and only after the final
maximal localization the anti-bonding shape is recovered. Note that, in this
simple case of silicon, for both valence and conduction manifolds we reach the
same set of WFs whether we run a maximal localization directly after parallel
transport, or a maximal localization after parallel transport + single
rotation, so the two cases are merged into one column in Fig. 1 under the
header PT(+SR)+ML. The total spreads for the VB and the CB manifolds after
parallel transport + maximal localization are
$8.584\text{\,}\mathrm{\text{Å}}^{2}$ and
$31.899\text{\,}\mathrm{\text{Å}}^{2}$, respectively. As expected, their sum
is larger than the value for the VCB manifold
($29.619\text{\,}\mathrm{\text{Å}}^{2}$) after disentanglement and maximal
localization, since in the VCB case there is additional freedom to further
minimize the spread by remixing bonding and anti-bonding WFs into pure $s$ and
$p$ orbitals (we highlight that using atom-centered $s,p$ projections does not
lead to the most localized orbitals for VCB in silicon; with a choice of atom-
centered $sp^{3}$ projections, the total spreads can further decrease to
$26.761\text{\,}\mathrm{\text{Å}}^{2}$, where four WFs have spreads
$3.522\text{\,}\mathrm{\text{Å}}^{2}$ and another four
$3.168\text{\,}\mathrm{\text{Å}}^{2}$). In addition, we note that since the
Hamiltonian diagonalization returns a random gauge, the spreads for parallel
transport and parallel transport + single rotation are different in each run,
but the spreads of PT(+SR)+ML should always be the same, since the algorithm
should always manage to find the maximally-localized gauge in this simple
case. To quantify how our multi-step procedure increases the overall
robustness of the algorithm while at the same time reducing its computational
cost, we show in SI LABEL:supp-fig:si2_val_4gauge the evolution of WF spreads
and centers during maximal localizations. Starting from the random gauge
directly after Hamiltonian diagonalization (LABEL:supp-
fig:si2_val_4gauge&rand_gauge), it takes $28,430$ iterations to converge;
using random Gaussians as initial guesses (LABEL:supp-
fig:si2_val_4gauge&rand_proj), the number of iterations decreases
significantly to $812$; with the parallel transport gauge (LABEL:supp-
fig:si2_val_4gauge&pt), the number of iterations further decreases to $228$,
and the evolution of spreads and centers is much smoother; the best starting
gauge is the one after single rotation (LABEL:supp-fig:si2_val_4gauge&pt_or),
which only takes $40$ iterations to converge, without any oscillations in the
evolution of spreads and centers. Note that the spreads of valence MRWFs from
PT(+SR)+ML are the same as MLWFs obtained from a direct Wannierization of the
valence bands, i.e., the valence MRWFs after separation span the original DFT
valence manifold, thus the initial VCB Wannierization does not cause any
delocalization of the valence MRWFs (see SI LABEL:supp-sec:MRWF_gauge for a
proof).

We now discuss the quality of the band interpolation. The WFs for VB & CB are
constructed from the initial VCB manifold obtained from a preliminary
disentanglement and maximal localization. Therefore, if VB/CB are properly
Wannierized, their band interpolation quality should be similar to that of VCB
MLWFs. Thus, in the following paragraphs, we compare the band interpolations
of VB/CB MLWFs, VCB MLWFs, and DFT bands. Once the starting VCB manifold was
properly disentangled and could well reproduce the DFT band structure, as
shown in Fig. 1, the WFs after parallel transport + maximal localization for
the VB and the CB manifolds can reproduce the corresponding part of the VCB
Wannier interpolated bands with high accuracy, being visually
indistinguishable. To quantitatively evaluate the band interpolation quality,
we compute the average band distance, $\eta_{\mathrm{isolated}}$, between the
VCB and VB/CB bands [31, 12, 17]. The $\eta_{\mathrm{isolated}}$ is defined as

$\eta_{\mathrm{isolated}}^{\mathrm{A-B}}=\sqrt{\frac{\sum_{n\mathbf{k}}(\epsilon_{n\mathbf{k}}^{\mathrm{A}}-\epsilon_{n\mathbf{k}}^{\mathrm{B}})^{2}}{N_{b}N_{\mathbf{k}}}},$
(7)

where $\epsilon_{n\mathbf{k}}$ are the eigenvalues of a band structure, and
its superscript A or B refers to the eigenvalues of two different bands, A or
B, which can be DFT bands, or Wannier interpolated bands of VCB, VB, or CB;
$N_{b}$ and $N_{\mathbf{k}}$ are the number of bands and $k$-points,
respectively. For silicon, we obtain: $\eta_{\mathrm{isolated}}^{\mathrm{VCB-
VB}}=$ $6.6\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ and
$\eta_{\mathrm{isolated}}^{\mathrm{VCB-CB}}=$
$15.5\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$. In general, the VB
interpolation is more accurate than CB since the VB MLWFs usually have smaller
spreads. To improve the CB interpolation quality, one might need to increase
$k$-point sampling, as we discuss in Section II.3.2.

###  II.3 MoS2

Next, we test the method on a two-dimensional (2D) MoS2 monolayer. For VCB
Wannierization, we use the standard hydrogenic Mo $d$ and S $s,p$ projections
(the semicore states are excluded, so in total 9 VBs and the lowest 4 CBs are
Wannierized). Since the VB and the lowest four CBs of MoS2 are isolated, 13
WFs are maximally localized from 13 bands without disentanglement. The 13
MLWFs can be well characterized into 4 groups by their angular momentum: as
shown in the VCB column of Fig. 2, from top to bottom, 3 resemble $d_{z^{2}}$,
2 resemble $d_{xy}$, and the remaining 8 resembles $sp^{3}$ hybridized
orbitals.

####  II.3.1 Valence and conduction bands

For the valence manifold, both after PT and after PT+SR, the WFs still do not
have a clear resemblance to bonding orbitals; after PT+ML or PT+SR+ML, the WFs
can be well grouped into six hybrids of Mo $d_{z^{2}}$ \+ S $p$, two $s$-like
WFs near sulfur atoms, and one WF floating inside the hexagonal cage and
having $C_{3h}$ symmetry, originating from the hybridization of three properly
oriented Mo $d_{z^{2}}$ orbitals from the three nearest Mo atom. For the
conduction manifold, WFs after PT are already close to the anti-bonding hybrid
orbitals, and the further SR or ML steps help to slightly reduce the spreads
and result in more symmetrized WF shapes. We notice that in contrast to
intuition, in this case it took more iterations to converge starting from
PT+SR gauge than directly from the PT gauge, as shown in the SI LABEL:supp-
fig:MoS2_evolution. Although the PT+SR cases start from a smaller total
spread, the maximal localizations got stuck longer in plateaus in the final
stages, leading to longer iterations. However, in all cases, the evolutions
are smooth and converge in less than $500$ iterations, since both valence and
conduction are isolated manifolds themselves, and PT is able to construct good
starting gauges.

Figure 2: WF shapes and band interpolation of MoS2. 2 WF shape of various
Wannierization methods. The notations are the same as Fig. 1, except that in
some cases the additional star sign (*) indicates that the WFs are grouped
together if their spreads are roughly similar, and only one of the shapes is
shown: e.g., 4.560 (3)* means there are three WFs having similar spreads, and
their average spread is $4.560\text{\,}\mathrm{\text{Å}}^{2}$. The yellow and
the silver spheres are the S and Mo atoms, respectively. 2 Band structure
comparison of DFT, VCB, VB, and CB. The VCB are computed by Wannier
interpolation after maximal localization (ML); the VB and CB are computed by
Wannier interpolation from PT+ML. The DFT bands are almost indistinguishable
since the VCB, VB, and CB Wannier-interpolated bands overlap essentially
exactly with the DFT bands.

In terms of band interpolation, again the Wannier interpolated VB and CB
overlap essentially exactly with the DFT bands as well as the Wannier
interpolated VCB, as shown in Fig. 2, and demonstrated quantitatively by the
excellent band-distance values: $\eta_{\mathrm{isolated}}^{\mathrm{VCB-VB}}=$
$0.19\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ and
$\eta_{\mathrm{isolated}}^{\mathrm{VCB-CB}}=$
$0.51\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$.

####  II.3.2 Single top valence band

In practical applications, the highest valence and lowest conduction bands are
of high interest since they are critical for electronic transport properties.
However, the Wannierization of a single band remains elusive since it is
difficult to write down a proper initial projection, resulting from a complete
hybridization of many different atomic orbitals. However, in the MoS2 case,
since the top valence band is isolated with respect to all other bands, we can
use our algorithm to construct a smooth gauge for that single band,
demonstrating the more general applicability of our method, beyond the
separation of VB & CB.

Figure 3 shows the band interpolation of this single-band WF, and the inset
shows the shape of this highly symmetric WF in real space. As usual, since the
separate Wannierizations in each submanifold have less degrees of freedom
compared with the Wannierization of the initial manifold, the WF spreads for
separate Wannierization are usually larger. Indeed, the single WF has a
relatively large spread ($9.288\text{\,}\mathrm{\text{Å}}^{2}$). For such a
large spread, artificial interactions between periodic copies of the same WF
in different supercells (where the supercell size is determined by the
$k$-point sampling) may start to become non-negligible. Indeed, we observe
some small oscillations at the minimum of the band along $\Gamma$ to M and
along M to K, whose zoom-ins are shown in Figs. 3 and 3. By increasing the
$k$-point sampling from the $12\times 12\times 1$
($0.2\text{\,}\mathrm{\text{Å}}^{-1}$ density, same as the VCB Wannierization)
to $18\times 18\times 1$, the interpolation quality improves significantly
(see Figs. 3 and 3). This means that the band interpolation error is not
caused by our separation method, but by the insufficient $k$-point sampling.
Therefore, if one targets a very high band interpolation quality, the
$k$-point sampling might need to be increased when considering a submanifold
only.

Figure 3: Band structure comparison of the single MoS2 top valence band. 3
Comparison of DFT and Wannier-interpolated top VB. The bottom-left inset shows
the shape of the resulting single MLWF. 3 and 3 are zoom-in comparisons of the
Wannier-interpolated bands using different $k$-point samplings: red dashed
line for $12\times 12\times 1$, blue dotted dashed line for $14\times 14\times
1$, and green dotted line for $18\times 18\times 1$. The $16\times 16\times 1$
is not shown since it has similar quality to the $14\times 14\times 1$.

###  II.4 SrVO3

Here we test the present method on the metallic perovskite SrVO3, which is a
correlated material on the $t_{2g}$ Hubbard manifold. At the DFT level, there
is a small gap ($37.8\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ at R point)
separating the $3d$ manifold from bands above, thus satisfying the requirement
of isolated manifold for the present method. Starting from 24 MLWFs for VCB,
we generate two sets of MLWFs: 5 for $3d$ and 19 for the remaining manifold.
For brevity, we only show the comparison of DFT and Wannier-interpolated bands
in Fig. 4, while the respective WF shapes and spreads are shown in SI
LABEL:supp-fig:SrVO3_wf and LABEL:supp-fig:SrVO3_evolution. For the $3d$
manifold, the PT gauge is already quite close to the maximally localized
gauge: maximal localization only slightly decreases the total spread from
$9.815\text{\,}\mathrm{\text{Å}}^{2}$ of PT to
$9.629\text{\,}\mathrm{\text{Å}}^{2}$, and symmetrizes the shapes of WFs (the
$3d$ columns of LABEL:supp-fig:SrVO3_wf). For the remaining manifold, it is
quite hard to converge: only after $13,027$ iterations (LABEL:supp-
fig:SrVO3_evolution) the maximal localization can converge to real-valued,
spatially-symmetrized MLWFs (LABEL:supp-fig:SrVO3_wf). In this case, the
single rotation greatly helps in improving the convergence: only $1544$
iterations are needed to converge to the same MLWFs starting from PT+SR gauge,
also removing the oscillations in spreads and centers during maximal
localization (LABEL:supp-fig:SrVO3_evolution). For the band interpolation,
again the respective manifolds are accurately reproduced, as demonstrated by
the bands in Fig. 4 and the band distances:
$\eta_{\mathrm{isolated}}^{\mathrm{VCB-}3d}=$
$6.74\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ and
$\eta_{\mathrm{isolated}}^{\mathrm{VCB-others}}=$
$3.00\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$. Furthermore, since the $t_{2g}$
and $e_{g}$ manifolds are gapped in energy, we can also separate them into two
submanifolds. As shown in Fig. 4, the $t_{2g}$ and $e_{g}$ bands are again
reproduced very well.

Figure 4: Comparison of DFT and Wannier-interpolated bands for SrVO3. 4 The
VCB are computed by Wannier interpolation using PDWF; the $3d$ bands ($3d$)
and the remaining bands (others) are computed by Wannier interpolations from
PT+ML. At DFT level, the $3d$ bands have a small gap
($37.8\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ at R point) separated from the
bands above; the present algorithm successfully separates the $3d$ submanifold
from the remaining manifold. 4 The $t_{2g}$ and $e_{g}$ bands are further
separated starting from the $3d$ bands of 4.

###  II.5 Results on 77 insulators

Finally, we test our method on a set of 77 insulators with number of atoms in
the unit cell ranging from 1 to 45. This is the same as the insulator set of
Ref. [12], except that 4 (He, Ne, Ar2, Kr2) of the 81 materials are excluded
since they consist of closed-shell noble-gas atoms, where the valence
electrons are fully occupied (and there is thus no need for separate
Wannierizations). This comprehensive test set not only validates the
correctness of the present method, but also helps improve its generality to
cover edge cases (e.g., the additional treatment of $\mathbf{b}$-vectors in
the SI LABEL:supp-sec:bvec_pt) that would be difficult to discover with only a
few test cases. The separate Wannierization is implemented as a fully
automated AiiDA [32, 33, 34] workflow, which first runs the Wannierization of
VCB using PDWF [17], then splits the VCB manifold with the method discussed
here (see Section VI for the Wannier.jl code implementation), and then runs
two separate maximal localizations using Wannier90 for the VB and the CB
manifolds, respectively.

All the Wannierizations finish successfully and have excellent band
interpolation quality, which we measure by the band distance [31, 12, 17]
$\eta_{\mathrm{isolated}}$ for comparisons between isolated bands (VB of
separate Wannierization w.r.t. VB of DFT, VB of separate Wannierization w.r.t.
VB of VCB Wannierization, CB of separate Wannierization w.r.t. CB of VCB
Wannierization), and $\eta_{2}$ for comparisons involving CB of DFT (VCB
Wannierization w.r.t. VCB of DFT, CB of separate Wannierization w.r.t. CB of
DFT),

$\eta_{2}^{\mathrm{A-B}}=\sqrt{\frac{\sum_{n\mathbf{k}}\tilde{f}_{n\mathbf{k}}(\epsilon_{n\mathbf{k}}^{\mathrm{A}}-\epsilon_{n\mathbf{k}}^{\mathrm{B}})^{2}}{\sum_{n\mathbf{k}}\tilde{f}_{n\mathbf{k}}}},$
(8)

where
$\tilde{f}_{n\mathbf{k}}=\sqrt{f^{\mathrm{A}}_{n\mathbf{k}}(E_{F}+2,\sigma)f^{\mathrm{B}}_{n\mathbf{k}}(E_{F}+2,\sigma)}$
and $f(E_{F}+2,\sigma)$ is the Fermi-Dirac distribution with the Fermi energy
set to $2\text{\,}\mathrm{e}\mathrm{V}$ above the real Fermi energy $E_{F}$ to
compare also part of the conduction bands; the smearing width is set to
$\sigma=$ $0.1\text{\,}\mathrm{e}\mathrm{V}$. The statistics of $\eta$ are
shown in Fig. 5. It is worth noting that the VB interpolation
($\eta_{\mathrm{isolated}}^{\mathrm{VB-DFT}}$ =
$0.859\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$) is even more accurate than the
VCB interpolation ($\eta_{2}^{\mathrm{VCB-DFT}}$ =
$2.609\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$): this is partly because
$\eta_{2}^{\mathrm{VCB-DFT}}$ is averaged over all bands, including the larger
error of CB interpolation; to exclude the effect of averaging, we also compute
the band distance of the VB of VCB Wannierization w.r.t. the VB of DFT,
$\eta_{\mathrm{isolated}}^{\mathrm{VB(VCB)-DFT}}$ =
$1.721\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$, which is larger than
$\eta_{\mathrm{isolated}}^{\mathrm{VB-DFT}}$ =
$0.859\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$, showing that the VB
interpolation is indeed more accurate than the VCB interpolation—this can be
explained by two facts: (1) the valence MRWFs are constructed by unitary
transformations of Bloch states, thus the valence manifold is exactly
preserved (see SI LABEL:supp-sec:MRWF_gauge for a proof); (2) the valence
MRWFs are more localized than the VCB MLWFs (will be discussed in the next
paragraph), leading to a better Fourier interpolation quality. The CB
interpolation ($\eta_{2}^{\mathrm{CB-DFT}}$ =
$7.619\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$) is slightly worse than the CB
of VCB Wannierization ($\eta_{2}^{\mathrm{CB(VCB)-DFT}}$ =
$6.616\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$), since the CB MLWFs are more
delocalized than the VCB MLWFs (will be discussed in the next paragraph);
moreover, it appears much larger than the $\eta_{2}^{\mathrm{VCB-DFT}}$ =
$2.609\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ since $\eta$ is defined as an
average over all bands—the (accurate) VB interpolations are not taken into
account in the computation of $\eta_{2}^{\mathrm{CB-DFT}}$. In addition, and
most importantly, the VB/CB to VCB distances are
$\eta_{\mathrm{isolated}}^{\mathrm{VB-VCB}}$ =
$2.219\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ and
$\eta_{\mathrm{isolated}}^{\mathrm{CB-VCB}}$ =
$3.835\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$, showing that the submanifolds
are well separated with little loss of interpolation accuracy compared with
the starting-point VCB Wannierization. For completeness, we also show the
statistics of max band distance, which is a stricter measure of band
interpolation quality, in the SI LABEL:supp-fig:maxbandsdist.

Now we discuss the localization of MLWFs by comparing the average spread from
the 77 materials. For VCB Wannierization, the average spread
$\Omega^{\mathrm{VCB}}$ = $1.178\text{\,}\mathrm{\text{Å}}^{2}$; after
separation (followed by maximal localization of VB and CB, respectively), the
average spread of VB MLWFs is slightly more localized ($\Omega^{\mathrm{VB}}$
= $1.079\text{\,}\mathrm{\text{Å}}^{2}$); while that of CB MLWFs are more
delocalized ($\Omega^{\mathrm{CB}}$ = $2.919\text{\,}\mathrm{\text{Å}}^{2}$).
This is consistent with the intuition that the VB MLWFs are the more localized
bonding orbitals whereas the CB MLWFs are the more delocalized anti-bonding
orbitals. Finally, as discussed in the previous section, the separated MLWFs
have less degrees of freedom compared with the VCB MLWFs, thus the sum of the
spreads of VB and CB ($\Omega^{\mathrm{VB}}+\Omega^{\mathrm{CB}}$) is in
general larger than that of VCB ($\Omega^{\mathrm{VCB}}$). Figure 6 shows the
percentage increase of $\Omega^{\mathrm{VB}}+\Omega^{\mathrm{CB}}$ over
$\Omega^{\mathrm{VCB}}$. On average, there is a 52.9% increase of the spread.
Note, however, that in three cases (Na2Se, Ca4O14V4, and HK) there are a
39.2%, 56.0%, and 5.1% decreases of spreads, respectively. For the first two
cases, there are few large-spread VCB MLWFs, probably because the VCB
Wannierizations are trapped in local minima; during the separation
Wannierizations, the parallel transport algorithm is able to find a smoother
gauge, thus reaching more localized MLWFs for both VB and CB, leading to
smaller $\Omega^{\mathrm{VB}}+\Omega^{\mathrm{CB}}$ than
$\Omega^{\mathrm{VCB}}$.

For completeness, we show the band-structure comparisons, the band distances,
the evolution of spreads and the evolution of the sum of WF centers during
maximal localization in SI [35], for each of the 77 materials. The smooth
evolution of spreads and sum of WF centers during maximal localization for CB
and VB demonstrate that parallel transport is able to construct continuous
gauge, thus maximal localization has no difficulty in further smoothening the
gauge.

Figure 5: Band distances of 77 insulators. From top to bottom: histograms of
band distances for valence plus conduction bands (VCB) w.r.t. DFT, valence
bands (VB) of VCB Wannierization w.r.t. DFT, VB Wannierization w.r.t. DFT, VB
Wannierization w.r.t. VCB Wannierization, conduction bands (CB) of VCB
Wannierization w.r.t. DFT, CB Wannierization w.r.t. DFT, and CB Wannierization
w.r.t. VCB Wannierization. The red and blue bars are the histograms and
cumulative histograms, respectively. The vertical lines indicate the mean and
median values.

Figure 6: WF spreads of 77 insulators. 6 From top to bottom: histograms of
VCB, VB, and CB MLWF spreads. The upper limit of the right y-axis of each
panel is the total number of WFs, respectively. 6 Histogram of the percentage
increase of the sum of VB and CB spreads w.r.t. VCB spreads.

##  III Conclusions

We introduce an automated method (manifold-remixed Wannier functions (MRWF))
to separate band manifolds by constructing MLWFs for the respective
submanifolds that have finite energy gaps between them. The method naturally
extends to the case of valence and conduction manifolds, but also to any other
case of isolated groups of bands. First, we start with a properly Wannierized
valence plus conduction manifold, obtained using any manual or automated
method (the recently introduced projectability disentangled WF [17] is
particularly suitable for this application since it is able to robustly and
reliably construct MLWFs that preserve as much as possible the anti-bonding
characters). Then, we split the manifolds by diagonalizing the Wannier-gauge
Hamiltonian into submanifolds for target energy ranges, respectively. Next,
using parallel transport, we construct smooth gauges for each submanifold to
fix the randomness caused by the independent Hamiltonian diagonalization at
every $k$-point. Finally, we maximally localize the parallel-transport gauge
to obtain smooth MLWFs for the desired manifolds. Before the final maximal
localization, we can optionally run a preliminary rotation w.r.t. a single
unitary matrix to fix the remaining gauge randomness intrinsic to parallel
transport. Often the final maximal localization is able to find the maximally-
localized gauge directly; however, the single rotation step helps to improve
the robustness of the final maximal localization, and has the additional
benefit of improving localization while still preserving the parallel
transport gauge, which might be relevant in some applications.

Results on silicon and MoS2 show that the final valence (conduction) MLWFs
restore faithfully chemical intuition for bonding/anti-bonding orbitals, and
accurately reproduce the valence/conduction part of the band structure of the
valence plus conduction manifold. Moreover, we demonstrate that the method is
not limited to the separation of valence and conduction manifolds, but also
applicable to any system with band groups separated by a finite gap: for
instance, the single top valence band of MoS2; or the $3d$, $t_{2g}$, and
$e_{g}$ manifolds of SrVO3. Furthermore, we implement fully automated AiiDA
[32, 33, 34] workflows to carry out the whole separate Wannierization process,
and test the present method on a set of 77 insulators. Statistics show that
the band interpolation achieves excellent accuracy at the
$\mathrm{m}\mathrm{e}\mathrm{V}$ scale, and on average the sum of VB and CB
MLWF spreads increase around 50% w.r.t. the VCB MLWF spreads. Thus, we
highlight that to ensure accurate band interpolation quality, the $k$-point
sampling density for separated Wannierization might need to be increased, as
demonstrated in the Wannierization of the MoS2 top valence band.

As an outlook, we envision several applications that one may find useful with
the preset approach: the analysis of bonding/anti-bonding orbitals based on
MLWFs; material properties that rely solely on the occupied manifold, such as
the electric polarization; spectral theories that require separate sets of
localized orbitals for both occupied and unoccupied states (for instance, the
Koopmans functionals to predict accurately the electronic band gap [6]); and
the dynamical mean field theory for correlated electrons.

##  IV Methods

The DFT calculations are carried out by Quantum ESPRESSO [36], using the SSSP
efficiency (version 1.1, PBE functional) library[31] for pseudopotentials and
its recommended energy cutoffs. The high-throughput calculations for 77
insulators are managed by the AiiDA [32, 33, 34] infrastructure which submits
Quantum ESPRESSO and Wannier90 [4] calculations to remote clusters, parses,
and stores the results into a database, while also orchestrating all sequences
of simulations and workflows. The automated AiiDA workflows are open-source
and hosted on GitHub [37]. Semicore states from pseudopotentials are excluded
from Wannierizations, except for a few cases where the semicore states overlap
with valence states; in such cases, all the semicore states are Wannierized. A
regular $k-$point mesh is used for the Wannier calculations, with a $k-$point
spacing of $0.2\text{\,}{\mathrm{\text{Å}}}^{-1}$, as selected by the protocol
in Vitale  _et al._ [12]. Figures are generated by matplotlib [38].

##  V Data Availability

All data generated for this work can be obtained from the Materials Cloud
Archive (https://doi.org/10.24435/materialscloud:2f-hs).

##  VI Code Availability

All codes used for this work are open-source; the latest stable versions can
be downloaded at http://www.wannier.org/ for Wannier90, https://www.quantum-
espresso.org/ for Quantum ESPRESSO, https://www.aiida.net/ for AiiDA, and
https://github.com/aiidateam/aiida-wannier90-workflows for aiida-
wannier90-workflows.

The MRWF method is implemented in an open-source Julia [39] package named
Wannier.jl, which is available at https://github.com/qiaojunfeng/Wannier.jl,
and https://www.wannierjl.org/ for the accompanying documentation/tutorials.

##  VII Acknowledgements

We thank Antoine Levitt and Michael F. Herbst for helpful discussions and
feedback on the Wannier.jl implementation. We acknowledge financial support
from the NCCR MARVEL (a National Centre of Competence in Research, funded by
the Swiss National Science Foundation, grant No. 205602), the Swiss National
Science Foundation (SNSF) Project Funding (grant 200021E_206190 “FISH4DIET”).
The work is also supported by a pilot access grant from the Swiss National
Supercomputing Centre (CSCS) on the Swiss share of the LUMI system under
project ID “PILOT MC EPFL-NM 01”, a CHRONOS grant from the CSCS on the Swiss
share of the LUMI system under project ID “REGULAR MC EPFL-NM 02”, and a grant
from the CSCS under project ID s0178.

##  VIII Author Contributions

J. Q. implemented and tested the method. G. P. and N. M. supervised the
project. All authors analyzed the results and contributed to writing the
manuscript.

##  IX Competing Interests

The authors declare that there are no competing interests.

## References

  * Marzari and Vanderbilt [1997] N. Marzari and D. Vanderbilt, Maximally localized generalized Wannier functions for composite energy bands, Phys. Rev. B 56, 12847 (1997).
  * Souza  _et al._ [2001] I. Souza, N. Marzari, and D. Vanderbilt, Maximally localized Wannier functions for entangled energy bands, Phys. Rev. B 65, 035109 (2001).
  * Marzari  _et al._ [2012] N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, Maximally localized Wannier functions: Theory and applications, Rev. Mod. Phys. 84, 1419 (2012).
  * 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, Wannier90 as a community code: new features and applications, J. Phys.: Condens. Matter 32, 165902 (2020).
  * [5] G. S. H. Pau,  _Reduced basis method for quantum models of crystalline solids_ , Thesis.
  * De Gennaro  _et al._ [2022] R. De Gennaro, N. Colonna, E. Linscott, and N. Marzari, Bloch’s theorem in orbital-density-dependent functionals: Band structures from Koopmans spectral functionals, Phys. Rev. B 106, 035106 (2022).
  * Georges  _et al._ [1996] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions, Rev. Mod. Phys. 68, 13 (1996).
  * Kotliar  _et al._ [2006] 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, 865 (2006).
  * Maier  _et al._ [2005] T. Maier, M. Jarrell, T. Pruschke, and M. H. Hettler, Quantum cluster theories, Rev. Mod. Phys. 77, 1027 (2005).
  * Gull  _et al._ [2011] E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, and P. Werner, Continuous-time Monte Carlo methods for quantum impurity models, Rev. Mod. Phys. 83, 349 (2011).
  * Damle  _et al._ [2015] A. Damle, L. Lin, and L. Ying, Compressed Representation of Kohn–Sham Orbitals via Selected Columns of the Density Matrix, J. Chem. Theory Comput. 11, 1463 (2015).
  * Vitale  _et al._ [2020] V. Vitale, G. Pizzi, A. Marrazzo, J. R. Yates, N. Marzari, and A. A. Mostofi, Automated high-throughput Wannierisation, npj Comput. Mater. 6, 66 (2020).
  * Mustafa  _et al._ [2015] J. I. Mustafa, S. Coh, M. L. Cohen, and S. G. Louie, Automated construction of maximally localized Wannier functions: Optimized projection functions method, Phys. Rev. B 92, 165134 (2015).
  * Mahler  _et al._ [2022] A. Mahler, J. Z. Williams, N. Q. Su, and W. Yang, Wannier Functions Dually Localized in Space and Energy, arXiv (2022), arXiv:2201.07751 [cond-mat.mtrl-sci] .
  * Thygesen  _et al._ [2005a] K. S. Thygesen, L. B. Hansen, and K. W. Jacobsen, Partly Occupied Wannier Functions, Phys. Rev. Lett. 94, 026405 (2005a).
  * Thygesen  _et al._ [2005b] K. S. Thygesen, L. B. Hansen, and K. W. Jacobsen, Partly occupied Wannier functions: Construction and applications, Phys. Rev. B 72, 125119 (2005b).
  * Qiao  _et al._ [2023] J. Qiao, G. Pizzi, and N. Marzari, Projectability disentanglement for accurate and automated electronic-structure Hamiltonians (2023), arXiv:2303.07877 [physics.comp-ph] .
  * Gontier  _et al._ [2019] D. Gontier, A. Levitt, and S. Siraj-dine, Numerical construction of Wannier functions through homotopy, J. Math. Phys. 60, 031901 (2019).
  * Damle  _et al._ [2019] A. Damle, A. Levitt, and L. Lin, Variational Formulation for Wannier Functions with Entangled Band Structure, Multiscale Model. Simul. 17, 167 (2019).
  * Brouder  _et al._ [2007] C. Brouder, G. Panati, M. Calandra, C. Mourougane, and N. Marzari, Exponential Localization of Wannier Functions in Insulators, Phys. Rev. Lett. 98, 046402 (2007).
  * Cancès  _et al._ [2017] É. Cancès, A. Levitt, G. Panati, and G. Stoltz, Robust determination of maximally localized Wannier functions, Phys. Rev. B 95, 075114 (2017).
  * Cornean  _et al._ [2017] H. D. Cornean, D. Monaco, and S. Teufel, Wannier functions and $\mathbb{Z}_{2}$ invariants in time-reversal symmetric topological insulators, Rev. Math. Phys. 29, 1730001 (2017).
  * Damle  _et al._ [2017] A. Damle, L. Lin, and L. Ying, SCDM-k: Localized orbitals for solids via selected columns of the density matrix, J. Comput. Phys. 334, 1 (2017).
  * Cornean  _et al._ [2016] H. D. Cornean, I. Herbst, and G. Nenciu, On the Construction of Composite Wannier Functions, Ann. Henri Poincaré 17, 3361 (2016).
  * Mustafa  _et al._ [2016] J. I. Mustafa, S. Coh, M. L. Cohen, and S. G. Louie, Automated construction of maximally localized Wannier functions for bands with nontrivial topology, Phys. Rev. B 94, 125151 (2016).
  * Winkler  _et al._ [2016] G. W. Winkler, A. A. Soluyanov, and M. Troyer, Smooth gauge and Wannier functions for topological band structures in arbitrary dimensions, Phys. Rev. B 93, 035453 (2016).
  * Mogensen and Riseth [2018] P. K. Mogensen and A. N. Riseth, Optim: A mathematical optimization package for Julia, J. Open Source Software 3, 615 (2018).
  * Wang  _et al._ [2006] X. Wang, J. R. Yates, I. Souza, and D. Vanderbilt, Ab initio calculation of the anomalous Hall conductivity by Wannier interpolation, Phys. Rev. B 74, 195118 (2006).
  * Wang  _et al._ [2017] C. Wang, X. Liu, L. Kang, B.-L. Gu, Y. Xu, and W. Duan, First-principles calculation of nonlinear optical responses by Wannier interpolation, Phys. Rev. B 96, 115147 (2017).
  * Ibañez-Azpiroz  _et al._ [2018] J. Ibañez-Azpiroz, S. S. Tsirkin, and I. Souza, Ab initio calculation of the shift photocurrent by Wannier interpolation, Phys. Rev. B 97, 245143 (2018).
  * Prandini  _et al._ [2018] G. Prandini, A. Marrazzo, I. E. Castelli, N. Mounet, and N. Marzari, Precision and efficiency in solid-state pseudopotential calculations, npj Comput. Mater. 4, 72 (2018).
  * Pizzi  _et al._ [2016] G. Pizzi, A. Cepellotti, R. Sabatini, N. Marzari, and B. Kozinsky, AiiDA: automated interactive infrastructure and database for computational science, Comput. Mater. Sci. 111, 218 (2016).
  * Huber  _et al._ [2020] S. P. Huber, S. Zoupanos, M. Uhrin, L. Talirz, L. Kahle, R. Häuselmann, D. Gresch, T. Müller, A. V. Yakutovich, C. W. Andersen, F. F. Ramirez, C. S. Adorf, F. Gargiulo, S. Kumbhar, E. Passaro, C. Johnston, A. Merkys, A. Cepellotti, N. Mounet, N. Marzari, B. Kozinsky, and G. Pizzi, AiiDA 1.0, a scalable computational infrastructure for automated reproducible workflows and data provenance, Sci. Data 7, 300 (2020).
  * Uhrin  _et al._ [2021] M. Uhrin, S. P. Huber, J. Yu, N. Marzari, and G. Pizzi, Workflows in AiiDA: Engineering a high-throughput, event-based engine for robust and modular computational workflows, Comput. Mater. Sci. 187, 110086 (2021).
  * [35] See Supplemental Material at [URL will be inserted by publisher], xxx.
  * Giannozzi  _et al._ [2020] P. Giannozzi, O. Baseggio, P. Bonfà, D. Brunato, R. Car, I. Carnimeo, C. Cavazzoni, S. de Gironcoli, P. Delugas, F. F. Ruffino, A. Ferretti, N. Marzari, I. Timrov, A. Urru, and S. Baroni, Quantum ESPRESSO toward the exascale, J. Chem. Phys. 152, 154105 (2020).
  * [37] aiida-wannier90-workflows: A collection of advanced automated workflows to compute Wannier functions using AiiDA and the Wannier90 code, https://github.com/aiidateam/aiida-wannier90-workflows, [Online; accessed 2023-05-24].
  * Hunter [2007] J. D. Hunter, Matplotlib: A 2D Graphics Environment, Comput. Sci. Eng. 9, 90 (2007).
  * Bezanson  _et al._ [2017] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, Julia: A Fresh Approach to Numerical Computing, SIAM Rev. 59, 65 (2017).