$\\#{\mathcal{B}}\leq\textnormal{OPT}({\mathcal{I}})+1+\\#\textnormal{{FFD}}({\mathcal{I}}\cap(T_{{\mathcal{I}}}\setminus
W))\leq\textnormal{OPT}({\mathcal{I}})+2=6.$

The first inequality holds by Lemma 39, Step 1, and Step 2 of Algorithm 11.
The second inequality holds by Step 4 and Step 5 of Algorithm FFD and that
$s(X_{I}\cap T_{{\mathcal{I}}})\leq 1$ (i.e., no second bin is opened in the
course of FFD as all items fit in a single bin). The third inequality holds
since $\textnormal{OPT}({\mathcal{I}})=4$.

  2. 2.

$s(Y_{I}\cap T_{{\mathcal{I}}})\leq 1$. Then, $\\#{\mathcal{B}}\leq 6$ by
symmetric arguments to Case 1.

  3. 3.

$s(Y_{I}\cap T_{{\mathcal{I}}})>1$ and $s(X_{I}\cap T_{{\mathcal{I}}})>1$. The
proof follows by Claim 44. ∎

We define
$t=\left\lfloor\frac{2\cdot\textnormal{OPT}({\mathcal{I}})}{3}\right\rfloor$
as the extra number of bins allowed to reach our approximation guarantee.

###### Claim 46

If $5\leq\textnormal{OPT}({\mathcal{I}})\leq 100$,
$\\#\textnormal{{FFD}}(X_{I}\cap T_{{\mathcal{I}}})>t-1$, and
$\\#\textnormal{{FFD}}(Y_{I}\cap T_{{\mathcal{I}}})>t-1$, then
$(\\#\textnormal{{FFD}}(X_{I})-1)\cdot(1-{\varepsilon})\leq s(X_{I})$ and
$(\\#\textnormal{{FFD}}(Y_{I})-1)\cdot(1-{\varepsilon})\leq s(Y_{I})$.

###### Proof

We use the following inequality.

$\displaystyle s(B_{{\mathcal{I}}})\leq{}$ $\displaystyle
s(I)-s(T_{{\mathcal{I}}})$ (45) $\displaystyle\leq{}$
$\displaystyle\textnormal{OPT}({\mathcal{I}})-s(Y_{I}\cap
T_{{\mathcal{I}}})-s(X_{I}\cap T_{{\mathcal{I}}})$ $\displaystyle\leq{}$
$\displaystyle\textnormal{OPT}({\mathcal{I}})-2\cdot\left(t-1\right)\cdot(1-{\varepsilon})$
$\displaystyle\leq{}$ $\displaystyle\textnormal{OPT}({\mathcal{I}})-2\cdot
t+2+2{\varepsilon}\cdot(t-1)$ $\displaystyle\leq{}$
$\displaystyle\textnormal{OPT}({\mathcal{I}})-2\cdot\left(\frac{2}{3}\cdot\textnormal{OPT}({\mathcal{I}})-\frac{2}{3}\right)+2+200\cdot{\varepsilon}$
$\displaystyle\leq{}$ $\displaystyle
3+\frac{1}{3}+0.1-\frac{\textnormal{OPT}({\mathcal{I}})}{3}$
$\displaystyle\leq{}$ $\displaystyle
3+\frac{2}{3}-\frac{\textnormal{OPT}({\mathcal{I}})}{3}$ $\displaystyle\leq{}$
$\displaystyle 2$

The third inequality holds because $\\#\textnormal{{FFD}}(X_{I}\cap
T_{{\mathcal{I}}})>t-1$, $\\#\textnormal{{FFD}}(Y_{I}\cap
T_{{\mathcal{I}}})>t-1$, and by Claim 11. The fifth inequality holds since
$\left\lfloor\frac{2\cdot\textnormal{OPT}({\mathcal{I}})}{3}\right\rfloor\geq\frac{2\cdot\textnormal{OPT}({\mathcal{I}})}{3}-\frac{2}{3}$.
The sixth inequality holds since ${\varepsilon}<0.0001$. The last inequality
holds since $\textnormal{OPT}({\mathcal{I}})\geq 5$. Because
$\\#\textnormal{{FFD}}(X_{I}\cap T_{{\mathcal{I}}})>t-1$ and
$\\#\textnormal{{FFD}}(Y_{I}\cap T_{{\mathcal{I}}})>t-1$, then for
$\textnormal{OPT}({\mathcal{I}})\geq 5$ it holds that $s(T_{{\mathcal{I}}}\cap
X_{I})>1$ and $s(T_{{\mathcal{I}}}\cap Y_{I})>1$ (FFD returns a packing with a
single bin for instances with total size at most one). Thus, by (45) and that
$s(T_{{\mathcal{I}}}\cap X_{I})>1$ and $s(T_{{\mathcal{I}}}\cap Y_{I})>1$, by
Claim 11 the proof follows. ∎

###### Claim 47

If $5\leq\textnormal{OPT}({\mathcal{I}})\leq 100$, then Algorithm 11 returns a
packing of at most $\frac{5}{3}\cdot\textnormal{OPT}({\mathcal{I}})$ bins.

###### Proof

The proof follows by the several sub-cases below.

  1. 1.

$\\#\textnormal{{FFD}}(X_{I}\cap T_{{\mathcal{I}}})\leq t-1$. Then,

$\\#{\mathcal{B}}\leq\textnormal{OPT}({\mathcal{I}})+1+\\#\textnormal{{FFD}}({\mathcal{I}}\cap(T_{{\mathcal{I}}}\setminus
W))\leq\textnormal{OPT}({\mathcal{I}})+1+t-1=\textnormal{OPT}({\mathcal{I}})+t\leq\frac{5}{3}\cdot\textnormal{OPT}({\mathcal{I}}).$

The first inequality holds by Lemma 39, Step 1, and Step 2 of Algorithm 11.
The second inequality holds by the definition of Algorithm FFD and that
$\\#\textnormal{{FFD}}(X_{I}\cap T_{{\mathcal{I}}})\leq t-1$ (i.e., no second
bin is opened in the course of FFD).

  2. 2.

$s(Y_{I}\cap T_{{\mathcal{I}}})\leq 1$. Then,
$\\#{\mathcal{B}}\leq\frac{5}{3}\cdot\textnormal{OPT}({\mathcal{I}})$ by
symmetric arguments to Case 1.

  3. 3.

$\\#\textnormal{{FFD}}(X_{I}\cap T_{{\mathcal{I}}})>t-1$ and
$\\#\textnormal{{FFD}}(Y_{I}\cap T_{{\mathcal{I}}})>t-1$. Therefore,

$\displaystyle\\#{\mathcal{B}}\leq{}$
$\displaystyle\textsf{{Color\\_Sets}}({\mathcal{I}})$ $\displaystyle\leq{}$
$\displaystyle\\#\textnormal{{FFD}}(X_{I})\oplus\\#\textnormal{{FFD}}(Y_{I})$
$\displaystyle\leq{}$
$\displaystyle\frac{s(X_{I})}{1-{\varepsilon}}+1+\frac{s(Y_{I})}{1-{\varepsilon}}+1$
$\displaystyle\leq{}$
$\displaystyle(s(X_{I})+s(Y_{I}))\cdot(1+2\cdot{\varepsilon})+2$
$\displaystyle\leq{}$ $\displaystyle s(I)\cdot 1.0002+2$ $\displaystyle\leq{}$
$\displaystyle\textnormal{OPT}({\mathcal{I}})+0.0002\cdot\textnormal{OPT}({\mathcal{I}})+2$
$\displaystyle\leq{}$ $\displaystyle\textnormal{OPT}({\mathcal{I}})+3$
$\displaystyle\leq{}$
$\displaystyle\frac{5}{3}\cdot\textnormal{OPT}({\mathcal{I}}).$

The first inequality holds by Step 2 of Algorithm 11. The second inequality
holds by Step 4, Step 5, and Step 7 of Algorithm 1. The third inequality holds
by Claim 46. The sixth inequality holds since ${\varepsilon}=0.0001$. The last
inequalities hold since $5\leq\textnormal{OPT}({\mathcal{I}})\leq 100$.

∎

In the following proof, we use the next result of [8].

###### Claim 48

For all BPP instance ${\mathcal{I}}=(I,s,E)$, it holds that
$\textnormal{OPT}({\mathcal{I}}\cap X_{I})+\textnormal{OPT}({\mathcal{I}}\cap
Y_{I})\leq\frac{3}{2}\cdot\textnormal{OPT}({\mathcal{I}})+1$.

Now,

###### Claim 49

If $\textnormal{OPT}({\mathcal{I}})\geq 100$, Algorithm 11 returns a packing
of at most $\frac{5}{3}\cdot\textnormal{OPT}({\mathcal{I}})$ bins.

###### Proof

$\displaystyle\\#{\mathcal{B}}\leq{}$
$\displaystyle\textsf{{Color\\_Sets}}({\mathcal{I}})$ $\displaystyle\leq{}$
$\displaystyle\\#\textnormal{{AsymptoticBP}}(X_{I})\oplus\\#\textnormal{{AsymptoticBP}}(Y_{I})$
$\displaystyle\leq{}$ $\displaystyle
1.02\cdot\textnormal{OPT}({\mathcal{I}}\cap
X_{I})+1.02\cdot\textnormal{OPT}({\mathcal{I}}\cap Y_{I})$ $\displaystyle={}$
$\displaystyle 1.02\cdot\left(\textnormal{OPT}({\mathcal{I}}\cap
X_{I})+\textnormal{OPT}({\mathcal{I}}\cap Y_{I})\right)$ $\displaystyle\leq{}$
$\displaystyle 1.02\cdot(1.5\cdot\textnormal{OPT}({\mathcal{I}})+1)$
$\displaystyle\leq{}$ $\displaystyle
1.53\cdot\textnormal{OPT}({\mathcal{I}})+1.02$ $\displaystyle\leq{}$
$\displaystyle
1.53\cdot\textnormal{OPT}({\mathcal{I}})+0.02\cdot\textnormal{OPT}({\mathcal{I}})$
$\displaystyle\leq{}$ $\displaystyle 1.55\cdot\textnormal{OPT}({\mathcal{I}})$
$\displaystyle\leq{}$
$\displaystyle\frac{5}{3}\cdot\textnormal{OPT}({\mathcal{I}}).$

The first inequality holds by Step 1 and Step 2 of Algorithm 11. The second
inequality holds by Step 4, Step 5, and Step 7 of Algorithm 1. The third
inequality holds by Lemma 3. The fourth inequality holds by Claim 48. The
sixth inequality uses $100\leq\textnormal{OPT}({\mathcal{I}})$.

∎

Finally, the proof of Theorem 0.F.1 follows by Claim 43, Claim 45, Claim 47,
and Claim 49. $\blacksquare$

##  Appendix 0.G Complete Multi-partite Conflict Graphs

In this section we give an asymptotic lower bound for Bin Packing with multi-
partite conflict graph (BPM), and show that it matches the (absolute)
approximation guarantee of Algorithm Color_Sets. The main idea is to use the
$\frac{3}{2}$-hardness of approximation result for classic bin packing.
Specifically, we duplicate a BP instance ${\mathcal{I}}$ with
$\textnormal{OPT}({\mathcal{I}})=2$ to generate a BPM instance. Interestingly,
we derive an asymptotic hardness result for the BPM instance using the
absolute hardness of BP. The next results resolve the complexity status of
BPM.

###### Lemma 50

For any $\alpha<\frac{3}{2}$, there is no asymptotic $\alpha$-approximation
for BPM, unless P=NP.

###### Lemma 51

Algorithm 1 is a $\frac{3}{2}$-approximation for BPM.

We now prove Lemmas 50 and 51. Recall that a multipartite graph $G=(V,E)$ has
a partition $V_{1},\ldots,V_{k}$ of $V$ such that for all $i\in[k]$ it holds
that $V_{i}$ is an independent set and for all $u\in V_{i},v\in
V_{j},i,j\in[k],i\neq j$ it holds that $(u,v)\in E$. Therefore, for simplicity
we use $\mathcal{I}=\big{(}I,\\{G_{1},\ldots,G_{n}\\},s\big{)}$ to denote a
BPM instance, where $I$ is a set of items, $\\{G_{1},\ldots,G_{n}\\}$ is the
unique partition of $I$ into independent sets, and $s$ is the size function.
Also, let ${\mathcal{I}}_{j}=(G_{i},s)$ be the corresponding BP instance to
the $j$-th set $G_{j}$. Finally, given a packing ${\mathcal{A}}$ w.r.t.
${\mathcal{I}}$ let ${\mathcal{A}}(G_{j})$ be the packing w.r.t.
${\mathcal{I}}_{j}$ containing all bins in ${\mathcal{A}}$ that contain items
from $G_{j}$ (observe that each bin in ${\mathcal{A}}$ contains item from only
one $G_{j}$ by the definition of complete multipartite graph).

###### Lemma 52

For any BPM instance $\mathcal{I}=\big{(}I,\\{G_{1},\ldots,G_{n}\\},s\big{)}$
it holds that

$\textnormal{OPT}(\mathcal{I})=\sum_{j\in[n]}\textnormal{OPT}({\mathcal{I}}_{j}).$

###### Proof

Let ${\mathcal{A}}=(A_{1},\ldots,A_{m})$ be an optimal packing for
$\mathcal{I}$. Since for each $j\in[n]$ the items of $G_{j}$ must be packed in
bins that contain only items from $G_{j}$, there must be at least
$\textnormal{OPT}({\mathcal{I}}_{j})$ such bins in the packing
${\mathcal{A}}$. If there are more than $\textnormal{OPT}({\mathcal{I}}_{j})$
such bins it follows that ${\mathcal{A}}$ is not optimal for $\mathcal{I}$.
$\blacksquare$

###### Lemma 53

For any BPM instance $\mathcal{I}=\big{(}I,\\{G_{1},\ldots,G_{n}\\},s\big{)}$
such that for all $j\in[n]$ it holds that
$\textnormal{OPT}({\mathcal{I}}_{j})=2$, and a packing
${\mathcal{A}}=(A_{1},\ldots,A_{m})$ of $\mathcal{I}$ with
$\\#{\mathcal{A}}<1.5\cdot\textnormal{OPT}(\mathcal{I})$, there is $j\in[n]$
such that $|{\mathcal{A}}(G_{j})|\leq 2$ and ${\mathcal{A}}(G_{j})$ is a
packing for $G_{j}$.

###### Proof

Assume towards a contradiction that for all $j\in[n]$ it holds that
${\mathcal{A}}(G_{j})\geq 3$. Therefore,

$m=\sum_{j\in[n]}{\mathcal{A}}(G_{j})\geq
3n=\frac{3}{2}\sum_{j\in[n]}OPT(G_{j})=\frac{3}{2}OPT(\mathcal{I}).$ (46)

The first equality is because in BPM each bin contains items only from one
group. The first inequality is by the assumption that
${\mathcal{A}}(G_{j})\geq 3$. The second equality is because for all $j\in[n]$
it holds that $OPT(G_{j})=2$. The last equality is by Lemma 52. By (46) it
holds that $m\geq 1.5\cdot\textnormal{OPT}(\mathcal{I})$ in contradiction.
$\blacksquare$

We define a duplication of a BP instance, which is simply duplicate the set of
items and giving each item the original size. Formally,

###### Definition 6

Let $G=(I,s)$ be a BP instance. A duplication of $G$ is a BP instance
$G^{\prime}=(I^{\prime},s^{\prime})$ where there is a bijection
$\sigma:I^{\prime}\rightarrow I$ such that for all $i\in I^{\prime}$ it holds
that $s(i)=s\left(\sigma(i)\right)$.

###### Lemma 54

Let $G=(I,s)$ be a BP instance and $G^{\prime}=(I^{\prime},s^{\prime})$ be a
duplication of $G$. Then, there is a polynomial time algorithm that given a
packing of $G^{\prime}$ with size $m$ finds a packing of $G$ with size $m$.

###### Proof

Let $(A^{\prime}_{1},\ldots,A^{\prime}_{m})$ be a packing of $G^{\prime}$ and
let $\sigma:I^{\prime}\rightarrow I$ be the bijection promised by Definition
6. Observe that $\sigma$ can be found in polynomial time by sorting the items
according to their sizes. For all $a\in[m]$ define
$A_{a}=\\{\sigma(i)~{}|~{}i\in A^{\prime}_{a}\\}$. Note that for all $i\in I$
there is $a\in[m]$ such that $i\in A_{a}$ since there is $i^{\prime}\in
I^{\prime}$ such that $\sigma(i^{\prime})=i$ because $\sigma$ is a bijection.
In addition, for all $a\in[m]$ it holds that $s(A_{a})=s(A^{\prime}_{a})\leq
1$. The first equality is by the definition of $\sigma$. Hence, it follows
that $(A_{1},\ldots,A_{m})$ is a packing of $G$. $\blacksquare$

###### Lemma 55

For any $n\in\mathbb{N}$, it is NP-hard to find a packing with size $2$ of a
BP instance $G=(I,s)$ with $|I|\geq n$ and $\textnormal{OPT}(G)=2$.

###### Proof

We use a reduction from BP, which is known to be NP-hard even for instances
$G$ with $OPT(G)=2$ [10]. The reduction is defined as follows given a BP
instance $G=(I,s)$ with $OPT(G)=2$. If $|I|\geq n$, then return
$G^{\prime}=G$. Otherwise, let $G^{\prime}=(I\cup I^{\prime},s^{\prime})$
where $|I^{\prime}|=n-|I|$, for all $i\in I$ it holds that
$s^{\prime}(i)=s(i)$, and for all $i\in I^{\prime}$ it holds that
$s^{\prime}(i)=0$. Given a packing of $G^{\prime}$ with size $m$ we can find a
packing to $G$ with size $m$, by removing the items from $I^{\prime}$.
Moreover, given a packing of $G$ with size $m$ we can find a packing to
$G^{\prime}$ with size $m$ by adding each $i\in I^{\prime}$ to one of the bins
in the packing arbitrarily. This results in a feasible packing since the size
of the items in $I^{\prime}$ is $0$. The claim follows. $\blacksquare$

Proof of Lemma 50: Assume towards a contradiction that there is $\alpha<1.5$
such that there is an asymptotic $\alpha$-approximation for BPM, and denote by
$X$ the algorithm achieving the above approximation. Denote by
$X(\mathcal{I})$ the packing returned by $X$ on some instance $\mathcal{I}$
and let $c=\frac{1.5-\alpha}{2}$. By the definition of asymptotic
$\alpha$-approximation, there is $n\in\mathbb{N}$ such that for all
$n^{\prime}\geq n$ and a BPM instance $\mathcal{I}$ with
$\textnormal{OPT}(\mathcal{I})=n^{\prime}$ it holds that

$\\#X(\mathcal{I})\leq\alpha\cdot\textnormal{OPT}(\mathcal{I})+c\cdot\textnormal{OPT}(\mathcal{I})=(\alpha+c)\cdot\textnormal{OPT}(\mathcal{I})<1.5\cdot\textnormal{OPT}(\mathcal{I}).$
(47)

The first inequality is because $X$ is an asymptotic $\alpha$-approximation
for BPM. The second inequality is because $\alpha<1.5$. Let
${\mathcal{J}}=(I,s)$ be a BP instance such that $|I|\geq n$. Let
$(G_{1},s),\ldots,(G_{n},s)$ be $n$ distinct duplications of ${\mathcal{I}}$
and let $\mathcal{I}_{n}=\big{(}I_{n},\\{G_{1},\ldots,G_{n}\\},s\big{)}$ where
$I_{n}=\bigcup_{j\in[n]}G_{j}$. Note that the construction of
$\mathcal{I}_{n}$ can be done in polynomial time because $|I|\geq n$. Let
${\mathcal{A}}=X(\mathcal{I}_{n})$ and let
${\mathcal{A}}=(A_{1},\ldots,A_{m})$. By Lemma 52 it follows that
$\textnormal{OPT}(\mathcal{I}_{n})\geq n$ and by Equation 47 it holds that
$|A|<1.5\cdot\textnormal{OPT}(\mathcal{I})$. Therefore, by Lemma 53 there is
$j\in[n]$ such that $|{\mathcal{A}}(G_{j})|\leq 2$ and that
${\mathcal{A}}(G_{j})$ is a packing for ${\mathcal{I}}_{j}$. Then, using
${\mathcal{A}}(G_{j})$, we can find in polynomial time a packing for
${\mathcal{J}}$ with size at most $2$ by Lemma 54. Unless P=NP, this is a
contradiction by Lemma 55. $\blacksquare$

Proof of Lemma 51: Let $\mathcal{I}=(I,\\{G_{1},\ldots,G_{n}\\},s)$ be a BPM
instance. Then,

$\\#\textsf{{Color\\_Sets}}({\mathcal{I}})=\sum_{j\in[n]}\\#\textsf{FFD}({\mathcal{I}}_{j})\leq\sum_{j\in[n]}\frac{3}{2}\textnormal{OPT}({\mathcal{I}}_{j})=\frac{3}{2}\sum_{j\in[n]}\textnormal{OPT}({\mathcal{I}}_{j})=\frac{3}{2}\textnormal{OPT}(\mathcal{I}).$

The first equality is by Algorithm 1. The first inequality is is since FFD is
a $\frac{3}{2}$-approximation for BP. The last equality is by Lemma 52.
$\blacksquare$

###### Contents

  1. 1 Introduction
  2. 2 Preliminaries
  3. 3 Approximations for Perfect and Bipartite Graphs
  4. 4 Split Graphs
  5. 5 Asymptotic Hardness for Bipartite and Split Graphs
  6. 6 Discussion
  7. 0.A Related Work
  8. 0.B Omitted Proofs from Section 2
  9. 0.C Omitted Proofs from Section 3
  10. 0.D Omitted Proofs from Section 4
  11. 0.E Omitted Proofs from Section 5
  12. 0.F A $\frac{5}{3}$-approximation for BPB
    1. 0.F.1 Deferred Proofs from Section 0.F
  13. 0.G Complete Multi-partite Conflict Graphs


# Spread balanced Wannier functions: Robust and automatable orbital
localization

Pietro F. Fontana  CAMD, Department of Physics, Technical University of
Denmark, DK-2800 Kongens Lyngby, Denmark  Ask H. Larsen  CAMD, Department of
Physics, Technical University of Denmark, DK-2800 Kongens Lyngby, Denmark
Thomas Olsen  CAMD, Department of Physics, Technical University of Denmark,
DK-2800 Kongens Lyngby, Denmark  Kristian S. Thygesen<EMAIL_ADDRESS>CAMD, Department of Physics, Technical University of Denmark, DK-2800 Kongens
Lyngby, Denmark

(August 27, 2024)

###### Abstract

We introduce a new type of Wannier functions (WFs) obtained by minimizing the
conventional spread functional with a penalty term proportional to the
variance of the spread distribution. This modified Wannierisation scheme is
less prone to produce ineffective solutions featuring one or several poorly
localized orbitals, making it well suited for complex systems or high-
throughput applications. Furthermore, we propose an automatable protocol for
selecting the initial guess and determine the optimal number of bands (or
equivalently WFs) for the localization algorithm. The improved performance and
robustness of the approach is demonstrated for a diverse set of test systems
including the NV center in diamond, metal slabs with atomic adsorbates,
spontaneous polarization of ferroelectrics and 30 inorganic monolayer
materials comprising both metals and semiconductors. The methods are
implemented in Python as part of the Atomic Simulation Environment (ASE).

††preprint: APS/123-QED

##  I Introduction

As computational codes are being increasingly automated it becomes possible to
perform complex, high-throughput investigations with minimal human efforts,
creating new advantages and opportunities. Within materials science, these
developments significantly expand the range of materials/properties that can
be examined by a single researcher. Moreover, it increases data quality by
reducing the risk of human errors, and enables researchers to address
materials phenomena or properties outside his/her domain expertise by lowering
barriers related to the technical aspects of the calculation. While some
computational tasks are straightforward to automate, others are more
challenging. An important example of the latter is the generation of localized
representations of the delocalized Bloch states of a crystal, i.e. Wannier
functions [1] (WFs). In a seminal paper, Marzari and Vanderbilt [2] introduced
a practical scheme for calculating maximally localized Wannier functions
(MLWFs) that overcomes the problem of the non-uniqueness (or ”gauge
dependence”). For simple systems, i.e. when the bands of interest form an
isolated group and/or there are few atoms in a unit cell, the standard
algorithms typically yield well localized WFs. In general, however, the
construction of useful WFs requires some hand-holding rendering automation
highly non-trivial.

When the bands of interest (from hereon referred to as the ”target bands”) are
isolated from all higher and lower lying bands by energy gaps, the MLWFs are
obtained by minimizing the sum of the quadratic spread of all the WFs. In the
general case, i.e. when the target bands do not form an isolated group, the
problem of finding a proper localized representation becomes significantly
harder. In this case, extra degrees of freedom (EDF) in the form of states
orthogonal to the target bands, must be included to aid the localization.
Since the target bands typically contain the occupied manifold, the problem of
identifying the optimum EDF can be seen as the process of augmenting the
target bands by their anti-bonding states [3].

The so-called disentanglement procedure [4] identifies the EDF by minimizing
the dispersion of the $\mathbf{k}$-subspaces (the subspace spanned by the
target bands and the EDF at a given $\mathbf{k}$) across the Brillouin zone
(BZ). Having identified the optimal $\mathbf{k}$-subspaces the MLWFs are
obtained following the usual localization procedure. For perfect crystals,
this approach is well suited. On the other hand, for systems where crystal
momentum is not a good quantum number, e.g. molecules, amorphous solids,
crystals with defects, etc. the idea of minimizing the $\mathbf{k}$-dispersion
does not seem a natural strategy. As an alternative, one can determine the EDF
by direct minimization of the spread functional. With this strategy, the
selection of the EDF and the localization into WFs is cast as one global
optimization problem rather the two-step strategy applied by the
disentanglement procedure. This idea leads to the partly occupied Wannier
functions (POWF) developed in Refs. [3, 5]. Being the result of a one- rather
than two-step optimization, the POWFs have smaller spread than the MLWFs.
Moreover, the scheme avoids reference to the $\mathbf{k}$-dispersion, which
seems natural for non-periodic systems. We note that the POWFs were
rediscovered in a different but equivalent form in Ref. [6].

Regardless of how the EDF are selected the standard localization procedure,
i.e. the minimization of the sum of the quadratic spreads, is not always
straightforward. One manifestation of this problem is its starting-point
dependence, i.e. that different WFs are obtained depending on the initial
guess for the orbitals, e.g. orbital type ($s$,$p$,$d$) and position (atom- or
bond centered). This indicates that the conventional spread functional
exhibits several local minima. From a practical point of view, the lack of
robustness/reproducibility arising from the starting point dependence is not a
problem in itself as long as a decent set of WFs is obtained. Unfortunately,
even that can be challenging and sometimes requires tuning of the initial
guess, the target bands, and number of EDF. One specific problem sometimes
encountered is that all WFs become well localized except for one or a few
which remain delocalized; this renders the entire set of WFs useless for many
purposes and significantly complicates automatization.

In this paper, we introduce a new class of spread functionals that explicitly
penalize delocalization of individual WFs. The minimization of these
functionals generally produces WFs with a more balanced spread distribution.
In particular, the problem of ”sacrificing” one WF to improve the total spread
does not occur. We also introduce a specific protocol for automatically
initializing the WFs based on the valence configuration of the involved atoms.
Leveraging the properties of the POWFs, we device a simple trial-and-error,
yet easy to automate, procedure for selecting the optimal number of EDF.
Putting it all together we arrive at a highly robust and fully automatic
scheme for constructing spread balanced WFs for general types of materials. We
demonstrate the method for a number of challenging systems, including atoms
adsorbed on metal slabs, the NV defect in diamond, and a set of 30 two-
dimensional (2D) materials arbitrarily selected from the Computational 2D
Materials Database (C2DB) [7]. The WFs are used to obtain electronic band
structures and spontaneous polarisations within the framework of the modern
theory of polarization [8]. All methods are implemented in the open source
Atomic Simulation Environment (ASE) [9].

##  II Theory

In this section we review the theory and construction of partly occupied WFs.
We then introduce our new spread functionals designed to produce WFs with
narrow size distributions. While we have investigated several different
functionals, we focus on the best performing one, the minimal variance spread
functional, throughout this paper. Finally, we describe our protocols for
initializing WFs and selecting the optimal number of WFs, respectively.

###  II.1 Partly occupied Wannier functions

The partly occupied Wannier functions were introduced in 2005 [3, 5] and were
recently demonstrated [6] to represent the global minimum of the quadratic
spread functional. The POWFs are related to the maximally localized Wannier
functions [2, 4] but avoid explicit reference to the wave vector in the band
disentanglement procedure, and are directly applicable to non-periodic
systems. Instead of maximizing the reciprocal space smoothness the POWF method
is entirely based on the minimization of the real space spread of the WFs.

For systems with periodic boundary conditions and a sufficiently large
supercell, the minimization of the conventional Marzari-Vanderbilt spread
functional [2] for a set of WFs $\\{w_{n}(\mathbf{r})\\}_{n=1}^{N_{w}}$ is
equivalent [10] to the maximization of

$\Omega=\sum_{n=1}^{N_{w}}\sum_{\alpha=1}^{N_{G}}W_{\alpha}|Z_{\alpha,nn}|^{2}$
(1)

where the matrix $Z_{\alpha}$ is defined as

$Z_{\alpha,nm}=\langle
w_{n}|e^{-i\mathbf{G}_{\alpha}\cdot\mathbf{r}}|w_{m}\rangle.$ (2)

The $\\{\mathbf{G}_{\alpha}\\}$ is a set of $N_{G}$ reciprocal lattice vectors
that connect each $k$-point to its neighbors and $W_{\alpha}$ are
corresponding weights accounting for the shape of the unit cell. The value of
$N_{G}$ can range from 3 to 6 depending on the symmetry of the unit cell. For
a discussion about these vectors and weights we refer to Ref. [11, 12]. The
definition of localization we impose with the functional $\Omega$, as in the
case of the Marzari-Vanderbilt spread functional, is equivalent to the Foster-
Boys method [13].

We emphasize that the assumption of a large supercell with periodic boundary
conditions does not represent a limitation. For example, it applies to a
pristine periodic crystal (with the primitive cell repeated a number of times
in all directions), isolated entities like molecules or clusters surrounded by
a sufficiently large vacuum region, a surface slab (possibly with the
primitive cell repeated in the in-plane directions), or a solid with
periodically repeated disorder, e.g. an impurity or point defect.

The goal of this approach is to obtain a set of $N_{w}$ localized WFs that can
reproduce any eigenstate below an energy threshold, $E_{0}$, exactly. Given
$N_{b}$ available eigenstates, a localization subspace is defined as the space
spanned by the $M$ eigenstates with energy below $E_{0}$ and additional $L$
extra degrees of freedom (EDF), where $M+L=N_{w}\leq N_{b}$. Each WF is then
defined as

$w_{n}=\sum_{m=1}^{M}U_{mn}\psi_{m}+\sum_{l=1}^{L}U_{M+l,n}\phi_{l}$ (3)

where the EDFs $\phi_{l}$ are defined as

$\phi_{l}=\sum_{m=1}^{N_{b}-M}c_{ml}\psi_{M+m}.$ (4)

The matrix $c$ has orthonormal columns while the matrix $U$ is unitary.

All the expressions in this section refer to the simple case where the
eigenstates have been obtained in a large supercell for which a $\Gamma$-point
sampling of the Brillouin zone (BZ) is a good approximation. We stress,
however, that for systems exhibiting periodicity on a smaller scale than the
supercell dimensions, e.g. for a perfect crystal, it is possible to formulate
the theory in terms of the eigenstates of the primitive unit cell sampled on a
uniform grid of $k$ points.[5] In this case, the number of fixed states,
$M^{\mathbf{k}}$, and EDF, $L^{\mathbf{k}}$, become $\mathbf{k}$-dependent.

The localization functional $\Omega$ can be maximized with respect to $U$ and
$c$ using any gradient-dependent algorithm under the constraint of
orthonormality of the EDFs, $\phi_{l}$, implemented e.g. via the method of
Lagrange multipliers. This step is referred to as the iterative localization
procedure. As noted above, the maximization of the functional $\Omega$ is
equivalent to the minimization of the Marzari-Vanderbilt spread functional
[2]. The latter often appears in the literature as $\Omega$ but we
nevertheless decided to use this symbol in order to keep the notation
consistent with the previous works on POWF [3, 5].

The POWF method was originally implemented in the open source ASE [9] Python
package. All the methods described in the following were implemented as
extensions/improvements to the existing POWF-ASE code.

###  II.2 Variance reducing spread functionals

The maximization of $\Omega$ in Eq. (1) is equivalent to the minimization of a
cost functional given by the sum of the quadratic spreads (second moments) of
the individual WFs. In our experience, this approach is not robust and can
produce delocalized WFs, in particular, for large numbers of WFs ($>50$). We
hypothetize that this happens because the cost of delocalizing a single
function can be compensated by a small improvement in the localization of a
number of other functions.

To circumvent this problem, we have explored different types of cost
functionals designed to share the spread more evenly across the entire set of
WFs. One approach is to apply a function on top of the $Z$-matrix elements of
the original functional (these are related to the inverse spread of the
corresponding WF)

$\Omega_{f}=\sum_{n=1}^{N_{w}}\sum_{\alpha=1}^{N_{G}}W_{\alpha}f\left(|Z_{\alpha,nn}|^{2}\right).$
(5)

We have tested different functions $f(x)$: square root ($\sqrt{x}$), scaled
error function ($\frac{2}{\sqrt{\pi}}\int_{0}^{2x}e^{-t^{2}}dt$), scaled and
translated sigmoid function ($1/\left[1+e^{-10(x-0.5)}\right]$), see Figure 1.
They span the same range as the original matrix elements, but introduce flat
plateaus for well localized and a steeper slope for delocalized WFs. The
effect of the steeper slope increases the gain of localizing a delocalized WF
relative to the cost of delocalizing an already localized WF. In particular,
the sigmoid function has an additional penalty for delocalized functions
($|Z_{\alpha,nn}|^{2}<0.5$; the threshold can be tuned if needed). We mention
that the modification function, $f$, may alternatively be applied to the
$\alpha$-sum instead of the individual $Z$-matrix elements, but this was not
pursued in the current study.

Figure 1:  Different functions applied to the $|Z|^{2}$ matrix elements
entering the localization functional Eq. (5). The original localization
functional is corresponds to $f(x)=x$ (blue dots).  Figure 2:  Effect of
changing the weight parameter, $w_{\mathrm{var}}$, of the variance penalty
term in $\Omega_{\mathrm{var}}$ for Wannier functions of bulk silicon. For
larger values of $w_{\mathrm{var}}$, the spread of the most delocalised WF
(maximal spread) decreases while the average spread of the WFs increases
slightly. The symbols represent the mean values over 10 independent
Wannierisations with different initial guess and the bars indicate the
standard deviations. The definition of the minimal $N_{w}$ is given in Eq. 7.

Our second approach adds a penalty term to the original functional
proportional to the variance of the spread distribution

$\mathrm{\Omega_{var}}=\Omega-
w_{\mathrm{var}}\mathrm{Var}\left[\sum_{\alpha=1}^{N_{G}}W_{\alpha}|Z_{\alpha,nn}|^{2}\right]$
(6)

where $w_{\mathrm{var}}$ is a parameter setting the weight for the variance
term. In all our calculations we have set $w_{\mathrm{var}}=N_{w}$ (the weight
of the penalty term should grow linearly with $N_{w}$ as the same holds for
$\Omega$). Of all the functionals tested in this study,
$\mathrm{\Omega_{var}}$ has the best overall performance, and we therefore
focus on this functional in the rest of the paper. In particular, the
functional based on the square root only led to minor and inconsistent
improvements in the localization compared to $\Omega$. The functionals based
on the error function and the sigmoid greatly increased the average spread
$\bar{s}$ in order to minimize the maximum spread $s_{\mathrm{max}}$ of the
set of WFs, see App. A for the definition of $\bar{s}$ and $s_{\mathrm{max}}$.
We stress that $\mathrm{\Omega_{var}}$ properly converges to real-valued WFs,
as expected when reaching the global maximum [14].

We now return to $\Omega_{\mathrm{var}}$ and the role of the weight parameter
of the penalty term, $w_{\mathrm{var}}$. Fig. 2 shows the average and maximum
spread of the set of WFs of bulk silicon obtained by maximizing
$\Omega_{\mathrm{var}}$ with different prefactors included in
$w_{\mathrm{var}}$. As expected, increasing the weight of the penalty term
leads to a more narrow spread distribution and thus a smaller spread of the
least localised WF (maximum spread), but at the same time leads to an increase
of the average spread of the WFs. Therefore, this parameter can be tuned as
needed and may even be optimized for specific applications. The initial value
of $w_{\mathrm{var}}=N_{w}$ worked consistently across the set of materials we
tested, hence we did not perform any further optimization.

The appearance of the penalty term in the functional Eq. 6 implies that the
mean of the spread distribution (the quantity minimized by the standard MLWFs)
is balanced with the variance of the spread distribution. For this reason we
refer to the WFs resulting from the maximization of the functional
$\mathrm{\Omega_{var}}$ as _spread balanced WFs_.

Finally, we note that the use of different objective functions may in general
lead to a different characters of the localized orbitals, as described in Sec.
IIIA of Ref. [14]. In particular, the $\Omega_{\mathrm{var}}$ functional may
produce WFs that differ somewhat in shape from those of the original $\Omega$
functional. A detailed investigation of this aspect is, however, beyond the
scope of the current study where we focus solely on the capability of
$\Omega_{\mathrm{var}}$ to produce WFs with narrow spread distributions.

###  II.3 Selecting the number of Wannier functions

Figure 3:  Left: Interpolated band structures for monolayer MoS2 (red dots)
for an energy threshold $E_{0}=\mathrm{CBM}+1$ $\mathrm{eV}$ (dotted black
line). The thin black lines represent the original Kohn-Sham band structure.
The interpolations have been performed with the WFs generated for $N_{w}=16$
and $N_{w}=17$, respectively. The former choice corresponds to the minimum
number of WFs required. Inclusion of one additional WF yields better
localization properties with a maximal spread of $s_{\mathrm{max}}=1.2$
$\mathrm{\SIUnitSymbolAngstrom}^{2}$ as compared to $s_{\mathrm{max}}=5.3$
$\mathrm{\SIUnitSymbolAngstrom}^{2}$. Right: Projection of eigenstates on the
subspace spanned by the WFs obtained for $N_{w}=16$ and $N_{w}=17$,
respectively. The dotted black line represents the energy threshold. The
projection is performed at the $\Gamma$ $k$-point.

For a given set of target bands, defined by the energy threshold $E_{0}$, the
number of WFs, $N_{w}=M^{\mathbf{k}}+L^{\mathbf{k}}$, is a most important
parameter for successful Wannerisation. Depending on the application,
different criteria may be used to quantify the quality of a set of WFs. We
have found that the spread (see App. A) of the most delocalized WF,
$s_{\mathrm{max}}$, is generally a good quality-indicator, and we will use
this measure along with the maximum and average band interpolation errors (see
App. B) throughout the paper.

In simple cases, a natural value of $N_{w}$ may be guessed by analysing the
band structure, considering symmetries, or using chemical intuition. In the
general case, however, the optimal $N_{w}$ cannot be guessed and a more
systematic approach is desired, see Sec. III.2. We note that the minimum
possible value for $N_{w}$ is given by the largest number of bands lying below
$E_{0}$ at any $\mathbf{k}$,

$N_{w}^{\mathrm{min}}(E_{0})=\max_{\mathbf{k}}{\sum_{n}H(E_{0}-\varepsilon_{n\mathbf{k}})}$
(7)

where $H$ is the Heaviside step function. At such $\mathbf{k}$-points we have
$M^{\mathbf{k}}=N_{w}^{\mathrm{min}}$ and thus $L^{\mathbf{k}}=0$, i.e. no
EDF. In the POWF formalism there is no upper limit to $N_{w}$ (apart from the
total number of bands available, $N_{b}$).

To illustrate how the $N_{w}$ parameter may influence the Wannierisation, we
consider the case of monolayer MoS2. In Fig. 3(left) we show the interpolated
band structure obtained from the POWFs obtained by maximising the $\Omega$
spread functional. Results are shown with an energy threshold ($E_{0}$) of 1
eV above the conduction band minimum (CBM) for $N_{w}=16$ and $N_{w}=17$,
respectively. (The choice $N_{w}=16$ corresponds to the minimum number of WFs
because exactly 16 bands fall below $E_{0}$ at approximately 1/3 along the
$\Gamma$-M path.) It is clear that the choice $N_{w}=17$ improves both the
band structure interpolation and $s_{\mathrm{max}}$. In the present case, this
result could have been anticipated by an analysis of the band structure, which
present an energy gap separating the lowest 17 bands from all higher lying
bands. Fig. 3(right) shows resolution of the resulting set of WF over the
eigenstates (at the $\Gamma$ point). From this analysis it is clear, that the
target states, consisting of states below $E_{0}$, can be perfectly completed
by including the lowest 17 bands in the Wannierization. In contrast, a
frustrated solution is obtained for $N_{w}=16$ where one of the EDF becomes a
mix of eigenstate 15 and 22.

In general, the effect of varying $N_{w}$ can be difficult to predict. A
prototype example where this happens is a non-elemental, low symmetry material
with no natural band gaps above $E_{0}$. In such cases, the POWF localization
algorithm, which may be seen as a bonding-antibonding completion procedure,
might select EDF corresponding to high-energy eigenstates well separated from
the target bands, or even mixtures of such [5].

###  II.4 Initial guess

The iterative localization procedure requires an initial guess for the
rotation matrix $U$ and the EDF coefficient matrix $c$. The quality of the
initial guess is essential. This is particularly true for systems with many
WFs where the iterative optimization algorithm is more likely to get trapped
in a local minimum if the initial guess is far from the global minimum. A
natural choice is to start from a set of $N_{\mathrm{AO}}$ atomic orbitals
$\\{g_{i}\\}$, and then project these onto the available eigenstates,
producing the $N_{b}\times N_{\mathrm{AO}}$ matrix
$P_{ni}=\langle\psi_{n}|g_{i}\rangle$. A prescription for extracting $U$ and
$c$ from $P$ can be found in Ref. [5]. We not that $N_{\mathrm{AO}}=N_{w}$ in
this procedure.

##  III Protocol for automated Wannierization

In this section we present a protocol for automated construction of POWF
suitable for high-throughput calculations. We stress that the protocol can be
used with any spread functional, e.g. $\Omega$ or $\mathrm{\Omega_{var}}$.

###  III.1 Initial guess

It is possible, and sometimes useful, to rely on chemical intuition when
selecting a set of atomic orbitals as initial guess for the Wannierisation. On
the other hand, chemical intuition is not easy to schematize in a form valid
for general materials. With automation a key motivation, we therefore propose
a simple, generally applicable protocol for an initial guess, which does not
require any parameters and which proved to be highly effective: For each atom
with d-states in its valence electron configuration, we include a group of 5
atom-centered d orbitals, one for each value of the magnetic quantum number.
This is motivated by our observation that in such systems one typically finds
WFs that closely mimick atom-centered $d$-orbitals. The total set of
$d$-orbitals ($N_{d}$ equals 5 times the number of transition metal atoms in a
unit cell) then sets a lower limit for $N_{w}$, in addition to the one set by
the threshold energy, $E_{0}$. The set of $d$-orbitals is complemented by
$N_{s}=N_{w}-N_{d}$ $s$-orbitals placed at random positions, but always within
a radius of $1.5\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ of an atom. These
$s$-orbitals can act as “nucleation centers” for atom- as well as bond-
centered $s,p$ or $sp$-like WFs. All the aforementioned atomic orbitals are
always set with a Gaussian of half-width
$1\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ as radial dependence. We did not
perform any tests on materials with $f$-electrons. However, due to the highly
localised nature of f-orbitals in general, we propose to treat such states
similarly to $d$-states, i.e. include a group of 7 atom-centered f orbitals in
the initial guess.

We have found that this choice of initial guess is both effective and robust
in the sense that leads to fairly rapid convergence toward solutions
representing either the global minimum or a local minimum close to the global
one. The variations induced by the randomness in the initial guess generally
produce small variations in the resulting WFs. Nevertheless, to reduce the
influence of the randomness we always perform five independent localizations
and pick the best solution (see caption of Figure 7 as example).

Figure 4:  Spread distribution of several calculations with different random
seeds for the initial guess. Each data point represents the result of a single
calculation, for a total of 100 for each functional.

We mention that an alternative method for initializing the Wannierisation is
the selected columns of the density matrix (SCDM) method proposed by Damle et
al. [15], that also does not require human intervention. The SCDM did,
however, not lead to improved performance over the initial guess protocol
described above when applied to a set of materials with isolated groups of
bands. It is not straightforward to combine the SCDM with the disentanglement
method of our POWF scheme and thus we were not able to test the SCDM for more
complex situations.

The element of randomness in our initialization scheme may seem to be a
weakness. However, we see it as a strength because it allows us to repeat the
Wannierisation with slightly different starting points, which often leads to
slightly different outcomes of which the ”best” solution can be selected. The
results in Fig. 4 show the presence of multiple nearby local maxima of both
$\Omega$ and $\Omega_{\mathrm{var}}$ for bulk silicon. In this case, all the
solutions are valid in the sense that the WFs are all well localized and real
valued. Nonetheless, in a given situation one may prefer a specific solution
satisfying certain problem specific requirements. The figure also clearly
reveals the consistent reduction of the maximum spread $s_{\mathrm{max}}$ when
using the $\mathrm{\Omega_{var}}$ functional, in particular when used with a
non-optimal $N_{w}$, such as the minimum value. The improvement becomes less
pronounced in this case when using the optimal $N_{w}$. A more detailed
discussion of these aspects are provided in Sec. IV.4.

###  III.2 Optimal number of Wannier functions

Figure 5:  Top: The spread of the most delocalized WF obtained for GaAs(FCC)
for different values of $N_{w}$. The calculations for each value of $N_{w}$
have been repeated 5 times with different initial orbitals. The circles
indicate the average value and the error bars are the standard deviation over
the 5 independent Wannierisations. Bottom: Interpolated band structure of GaAs
obtained using the minimal $N_{w}=6$ (as defined by the energy threshold of
$E_{0}=\mathrm{CBM}+2$ $\mathrm{eV}$) and the optimal $N_{w}=8$. The thin
black lines represent the Kohn-Sham band structure, the red dots are the
Wannier interpolation.

When the optimal number of WFs cannot be guessed, c.f. Sec. II.3, it must be
computed. We do this by constructing the POWFs for a range of $N_{w}\geq
N_{w}^{\mathrm{min}}$ and selecting the solution presenting the smallest
$s_{\mathrm{max}}$. In other words, we add EDF to the Wannierization space as
long as it reduces the spread of the least localized WF. For the calculations
presented in this work we have varied $N_{w}$ from $N_{w}^{\mathrm{min}}$ to
$N_{w}^{\mathrm{min}}+5$. Depending on the size of the system (number of
bands) and the available computational resources, a higher upper limit may be
chosen. Due to the randomness in the initial guess, we run 5 optimizations for
each value of $N_{w}$ and select the best solution. In total we thus perform
25 Wannierisations for each material. Based on our experience, the resulting
optimal $N_{w}$ is the same for $\Omega$ and $\mathrm{\Omega_{var}}$.

As an example, Fig. 5 shows $s_{\mathrm{max}}$ for bulk GaAs. The vertical
lines represent the variations due to the randomness in the initial guess. The
minimal $\langle s_{\mathrm{max}}\rangle$ appears for $N_{w}=8$. As was found
for MoS2, the optimal $N_{w}$ also produces a better interpolation of the band
structure, in particular around the energy threshold ($E_{0}$=CBM+2 eV), see
Fig. 5. We do stress, however, that the size of $s_{\mathrm{max}}$ is not
always directly correlated with the band interpolation error[16, 6].

##  IV Results

In this section we present the results we obtained with our methods. We start
presenting the differences we observed with the new localization functional on
few specific systems, then we move to a verification of the automated
procedure on a set of 30 two-dimensional materials.

###  IV.1 An illustrative example: WMo3Te8

In order to demonstrate the effect of the variance term in our newly defined
localization functional $\mathrm{\Omega_{var}}$, we construct the Wannier
functions of monolayer WMo3Te8. This material is a 2D semiconductor with 12
atoms per unit cell, 52 occupied bands, and a band gap of
$0.9\text{\,}\mathrm{eV}$ obtained from a DFT calculation with the PBE xc-
functional [17, 7].

In particular, we set an energy threshold $E_{0}=\mathrm{CBM}+2$ $\mathrm{eV}$
and computed 64 WFs, that is the minimal number of WFs needed to describe the
states up to $E_{0}$. The large number of WFs helps in proving our thesis, the
optimization algorithm is in fact more likely to choose a solution with few
delocalized WFs, if the relative contribution to the total localization
functional is lower.

For monolayer WMo3Te8 we obtain an average spread (see App. A) of
$\bar{s}=2.7$ $\mathrm{\SIUnitSymbolAngstrom}^{2}$ and a maximum spread of
$s_{\mathrm{max}}=21.5$ $\mathrm{\SIUnitSymbolAngstrom}^{2}$ with the standard
functional $\Omega$, while the variance minimizing functional
$\mathrm{\Omega_{var}}$ converges to $\bar{s}=2.8$
$\mathrm{\SIUnitSymbolAngstrom}^{2}$ and $s_{\mathrm{max}}=5.1$
$\mathrm{\SIUnitSymbolAngstrom}^{2}$. The most delocalized WFs for both
functionals are shown in Figure 6. For the band structure interpolation error
(see App. B) we obtain $\eta=21$ $\mathrm{meV}$ ($\eta_{\mathrm{max}}=420$
$\mathrm{meV}$) for the standard MLWFs produced by $\Omega$ and $\eta=9$
$\mathrm{meV}$ ($\eta_{\mathrm{max}}=143$ $\mathrm{meV}$) for the spread
balanced WFs generated with $\mathrm{\Omega_{var}}$.

As expected, the average localization of the WFs generated with the two
different spread functionals is almost identical while the spread of the most
delocalized WFs is greatly improved by $\mathrm{\Omega_{var}}$. Moreover, the
variance reducing functional leads to a significant improvement in the tight-
binding interpolation of the band structure. We stress that the significant
improvements found with $\mathrm{\Omega_{var}}$ for this specific material may
not be fully representative. Although we do generally find a significant
improvement, i.e. reduction, of $s_{\mathrm{max}}$, the band interpolation
error is often similar to that obtained with $\Omega$ (see Sec. IV.4).

The effect of the variance penalty term on the optimization of
$\mathrm{\Omega_{var}}$ is plotted in Figure 6 (c). Comparing the value of the
variance term between the two iterative optimizations, with $\Omega$ and
$\mathrm{\Omega_{var}}$, it is clear that the penalty term leads to a
different convergence path. For the vast majority of the materials we tested,
the number of steps required for the optimization of $\mathrm{\Omega_{var}}$
was comparable to or slightly larger than required for $\Omega$. However, due
to the additional terms in the gradient of $\mathrm{\Omega_{var}}$, each step
is roughly twice as expensive to evaluate in terms computational time.

Figure 6:  Isosurface plots of the most delocalized WF of monolayer WMo3Te8
obtained with (a) the standard localization functional $\Omega$ yielding a
spread of $s_{\mathrm{max}}=21.5$ $\mathrm{\SIUnitSymbolAngstrom}^{2}$ and (b)
the variance reducing functional $\mathrm{\Omega_{var}}$ yielding a spread of
$s_{\mathrm{max}}=5.1$ $\mathrm{\SIUnitSymbolAngstrom}^{2}$. The
Wannierisation has been performed for the minimal number of WFs consistent
with an energy threshold of $E_{0}=\mathrm{CBM}+2$ $\mathrm{eV}$. The
isosurface level is $0.7\text{\,}\mathrm{\SIUnitSymbolAngstrom}^{-3/2}$. (c)
The value of the variance penalty term in $\mathrm{\Omega_{var}}$ (see Eq.
(6)) during the iterative optimization. Additional optimization steps with
$\mathrm{\Omega_{var}}$ can further decrease the variance term by a factor of
2.

###  IV.2 Spontaneous polarization

The spontaneous polarization of ferrolectrics comprises a prominent example of
a physical quantity that is easily accessible from Wannier functions. As shown
by King-Smith and Vanderbilt [8] the change in polarization under an adiabatic
deformation can be calculated by a Berry-phase type formula. Typically the
polarization of ferroelectrics is measured with respect to a centrosymmetric
phase which is known to have vanishing polarization and the spontaneous
polarization can then be computed by a single calculation in the polar phase.

To unravel the relation to Wannier functions it is straightforward to show
that the Wannier charge centers can be written as [14]

$\displaystyle\mathbf{r}_{n}\equiv\langle
w_{n}|\mathbf{\hat{r}}|w_{n}\rangle=\frac{V}{(2\pi)^{3}}\int
d\mathbf{k}\langle
u_{n\mathbf{k}}|i\mathbf{\nabla}_{\mathbf{k}}|u_{n\mathbf{k}}\rangle,$ (8)

where $V$ is the unit cell volume and $u_{n\mathbf{k}}$ is the periodic part
of a Bloch function. Except for the factor of $V$ this expression is exactly
the Berry phase formula for the electric contribution to the polarization for
a single band and the full polarization can thus be written as

$\displaystyle\mathbf{P}=\frac{1}{V}\sum_{a}Z_{a}\mathbf{r}_{a}-\frac{1}{V}\sum_{n\in
occ}\mathbf{r}_{n},$ (9)

where $\mathbf{r}_{a}$ is the position of nucleus $a$ with charge $Z_{a}$. Eq.
(9) is formally equivalent to the Berry phase expression for the polarization
and the $2\pi$ ambiguity in the Berry phase is reflected by the fact that the
nuclei positions as well as the Wannier functions can be chosen in an
arbitrary unit cell. This expression for the polarization has the advantage of
providing a clear physical interpretation of the polarization as the dipole
resulting from a diplacement of Wannier charge centers from the nuclei
positions. Moreover, it strongly facilitates a microscopic analysis of how the
polarization is affected by impurities or interfaces since one may monitor the
shift in Wannier charge centers under an applied perturbation.

We have compared the calculated spontaneous polarization of tetragonal BaTiO3
obtained with a direct implementation of the Berry phase method [18] with that
obtained from Eq. (9) using the $\mathrm{\Omega_{var}}$ spread functional Eq.
(6) to construct Wannier functions from the occupied states. In the present
case we have performed a full relaxation with the PBE functional using a
$8\times 8\times 8$ $k$-point mesh and 800 eV plane wave cutoff. The result
from both calculations is 45.4 $\mu\textrm{C}/\textrm{cm}^{2}$, which is in
agreement with previous calculations [19, 20]. As expected, the new type of
spread balanced Wannier functions thus reproduce the result obtained with
standard maximally localized Wannier functions.

###  IV.3 Complex systems

|  | $\Omega$ | $\mathrm{\Omega_{var}}$
---|---|---|---
System | $N_{w}$ |  $\langle s_{\mathrm{max}}\rangle$ [$\mathrm{\SIUnitSymbolAngstrom}^{2}$] |  $\eta_{\mathrm{max}}$ [$\mathrm{meV}$] |  $\langle s_{\mathrm{max}}\rangle$ [$\mathrm{\SIUnitSymbolAngstrom}^{2}$] |  $\eta_{\mathrm{max}}$ [$\mathrm{meV}$]
NV center in diamond | 127 | $6.8\pm 0.3$ | $8$ | $4.46\pm 0.04$ | $8$
Adsorbed H on Ru slab | 133 | $8\pm 1$ | $158$ | $4.6\pm 0.5$ | $172$
Adsorbed N on Ru slab | 134 | $6.8\pm 0.7$ | $89$ | $4.6\pm 0.6$ | $82$
Adsorbed O on Ru slab | 135 | $6.9\pm 0.5$ | $77$ | $4.7\pm 0.4$ | $117$
Table 1:  Comparison of the spread of the most delocalized WF and the maximum
band interpolation error obtained for WFs generated with the standard spread
functional $\Omega$ and the variance reducing functional
$\mathrm{\Omega_{var}}$ for the NV center in diamond and H, N, O adsorbed on a
Ru(111) surface slab. The minimal $N_{w}$ was used for all materials. $\langle
s_{\mathrm{max}}\rangle$ is the average over 5 calculations with different
random seeds (standard deviation also shown). $\eta_{\mathrm{max}}$ is the
maximum band interpolation error for the most localized set of WFs (lowest
maximum spread) obtained in the 5 calculations with different random seeds.

While the generation of well localized Wannier functions is usually relatively
straightforward for simple systems (small number of atoms, isolated group of
bands, etc.), it can be significantly more challenging in the general case. To
test the variance reducing localization functional on more complex systems we
compare its performance to the standard $\Omega$ functional for a nitrogen-
vacancy (NV) defect center in a diamond crystal and a Ru(111) surface slab
with adsorbed H, N, and O atoms, respectively. The results are summarized in
Table 1 and confirm the previous conclusions. In particular, the spread of the
most delocalized WF is significantly reduced when using the
$\mathrm{\Omega_{var}}$ functional.

Additional computational details are provided in App. D.

###  IV.4 Towards high-throughput applications

Figure 7:  Results for the set of 2D materials with different choices for
$N_{w}$ and the localization functional. The energy threshold is set at CBM
($E_{F}$) + $2\text{\,}\mathrm{eV}$ for insulators (metals). Left column:
Spread of WFs for all the materials in the test set averaged over 5 runs with
different initializations. The standard deviation is indicated by lines. The
dashed green line is $s_{\mathrm{max}}=\bar{s}$, the dotted red line is
$s_{\mathrm{max}}=2\bar{s}$. The red circles indicate results for WMo3Te8 and
ZrTi3Te8. Right column: The maximum band interpolation error. Results from the
most localized set of WFs, out of the 5 calculations with different random
seeds.

In the previous sections we have demonstrated our new Wannierisation scheme
for several different types of systems. In this section, we apply the scheme
to a larger set of materials comprising 30 atomically thin two-dimensional
(2D) materials randomly selected from the Computational 2D Materials Database
[7]. The set includes 22 materials with finite band gaps and 8 metals, and
covers several different crystal lattices and a large set of chemical
elements. The complete list of materials is provided in App. E and the
computational details in App. C. Although we focus on 2D materials we expect
our results to be representative for general material types.

We applied our Wannierisation scheme to each material and determined the
optimal $N_{w}$ as the one yielding the lowest maximum spread,
$s_{\mathrm{max}}$, see Sec. II.3. For each value of $N_{w}$ in the range from
$N_{w}^{\mathrm{min}}$ to $N_{w}^{\mathrm{min}}+5$, we performed 5 independent
Wannierisations using different initializations (differing due to the
arbitrariness in the position of the $s$-orbitals). Following the procedure
visualised in Fig. 5(top), we used the average $s_{\mathrm{max}}$ over the 5
sets of WFs to determine the optimal $N_{w}$. In all cases, we include
eigenstates up to $2\text{\,}\mathrm{eV}$ above the conduction band minimum
(Fermi level) for insulating (metallic) materials.

The results of the Wannierisation procedure for the 30 materials are
summarised in Fig. 7. The left column shows the average and maximal spreads of
the WFs generated with the $\Omega$ and $\Omega_{\mathrm{var}}$ functionals
for the minimal and optimal number of WFs, respectively. Note the different
scales on the axes. The symbols indicate the mean values of $\bar{s}$ and
$s_{\mathrm{max}}$ over the 5 initializations, and the lines indicate the
standard deviation. The green and red dashed lines indicate
$\langle\bar{s}\rangle=\langle s_{\mathrm{max}}\rangle$ and
$\langle\bar{s}\rangle=2\langle s_{\mathrm{max}}\rangle$, respectively.

For generating well-localized WFs, the number of WFs, i.e. using $N_{w}^{opt}$
rather than $N_{w}^{min}$, was found to be more critical than the type of
spread functional, i.e. $\Omega$ versus $\Omega_{\mathrm{var}}$. However, the
use of $\Omega_{\mathrm{var}}$ instead of $\Omega$ does lead to a significant
improvement in the localisation for a fixed $N_{w}$. This improvement is most
pronounced for non-optimal values of $N_{w}$, e.g. the minimal $N_{w}$. In
particular, for the two materials WMo3Te8 and ZrTi3Te8 (both with $N_{w}>60$
and indicated by red circles), we were not able to localize all the WFs when
using the $\Omega$-functional and the minimal $N_{w}$. This issue did not
occur with the $\Omega_{\mathrm{var}}$-functional. Even when using the optimal
$N_{w}$, the improvement by the $\Omega_{\mathrm{var}}$-functional is
significant. Not only do we obtain better localisation of the least localised
WF ($s_{\mathrm{max}}$) without sacrificing the average localisation
($\bar{s}$), the standard deviations on both $\bar{s}$ and $s_{\mathrm{max}}$
is also lowered.

The band interpolation errors (again averaged over the 5 initializations) are
shown in the right column of Fig. 7. As previously found and discussed in
Section III.2, the errors decrease significantly when using the optimized
$N_{w}$ as compared to the minimal $N_{w}$. With the optimal $N_{w}$ the
majority of the materials show a maximum error below
$20\text{\,}\mathrm{meV}$. A few materials show higher band errors, which is
related to significant band crossings (band entanglement) with higher energy
bands close to the $E_{0}$ energy cutoff. We note in passing that if accurate
band interpolation in this region is required, one may simply increase
$E_{0}$, which will push the inaccuracies to higher band energies.

While the $\Omega_{\mathrm{var}}$ spread functional improves the localisation
properties of WFs, it does not represent a significant improvement over
$\Omega$ in terms of the band interpolation error. This shows that the band
interpolation error is not directly correlated with the localisation
properties of the WFs. More precisely, it is not directly correlated with the
spread of the most delocalised WF, $s_{\mathrm{max}}$, which is always
significantly and consistently reduced by using $\Omega_{\mathrm{var}}$. We
stress, however, that there are many other applications of WFs where robust
localisation of all WFs, is of key importance. These include the interpolation
of electron-phonon matrix elements[21, 22], calculation of Berry curvatures
and conductivities[23], and basis sets for electron transport calculations
based on non-equilibrium Green’s functions[24, 25].

Functional | $N_{w}$ |  $\langle\langle\bar{s}\rangle_{\mathrm{init}}\rangle_{\mathrm{mat}}$ [$\mathrm{\SIUnitSymbolAngstrom}^{2}$] |  $\langle\sigma_{\mathrm{init}}(\bar{s})\rangle_{\mathrm{mat}}$ [$\mathrm{\SIUnitSymbolAngstrom}^{2}$] |  $\langle\langle s_{\mathrm{max}}\rangle_{\mathrm{init}}\rangle_{\mathrm{mat}}$ [$\mathrm{\SIUnitSymbolAngstrom}^{2}$] |  $\langle\sigma_{\mathrm{init}}(s_{\mathrm{max}})\rangle_{\mathrm{mat}}$ [$\mathrm{\SIUnitSymbolAngstrom}^{2}$]
---|---|---|---|---|---
$\mathrm{\Omega}$ | Minimal | $1.44\pm 0.09$ | $0.06\pm 0.02$ | $4.8\pm 0.8$ | $1.0\pm 0.4$
$\mathrm{\Omega}$ | Optimal | $1.24\pm 0.06$ | $0.05\pm 0.01$ | $2.9\pm 0.2$ | $1.0\pm 0.4$
$\mathrm{\Omega_{var}}$ | Minimal | $1.46\pm 0.09$ | $0.08\pm 0.02$ | $3.2\pm 0.3$ | $0.5\pm 0.1$
$\mathrm{\Omega_{var}}$ | Optimal | $1.23\pm 0.06$ | $0.05\pm 0.01$ | $2.5\pm 0.1$ | $0.4\pm 0.1$
Table 2:  Comparison of the spread distribution over the entire set of 2D
materials. The definitions of $\bar{s}_{\mathrm{tot}}$,
$\bar{\sigma}_{\mathrm{tot}}$, $\bar{s}_{\mathrm{max}}$ and
$\bar{\sigma}_{\mathrm{max}}$ are given in Eqs. 10, 11 and the following text.
For each value we also report the standard deviation of the mean (i.e. the
standard error).

The effect of penalizing the spread variance in $\Omega_{\mathrm{var}}$ can be
illustrated by considering the localisation properties of the resulting sets
of WFs averaged over all 30 materials for a fixed $N_{w}$. Specifically, we
fix $N_{w}$ at either its minimal or optimal value, and consider the mean of
the average spread

$\langle\langle\bar{s}\rangle_{\mathrm{init}}\rangle_{\mathrm{mat}}=\frac{1}{N_{\mathrm{mat}}}\sum^{N_{\mathrm{mat}}}_{m=1}\left(\frac{1}{N_{\mathrm{init}}}\sum^{N_{\mathrm{init}}}_{i=1}\bar{s}_{i}^{m}\right)$
(10)

and the mean standard deviation of the average spread over 5 initialisations

$\langle\sigma_{\mathrm{init}}(\bar{s})\rangle_{\mathrm{mat}}=\frac{1}{N_{\mathrm{mat}}}\sum^{N_{\mathrm{mat}}}_{m=1}\sqrt{\mathrm{Var}_{\textrm{init}}[\bar{s}_{i}^{m}]},$
(11)

where $N_{\mathrm{mat}}$ is the total number of materials (here
$N_{\mathrm{mat}}=30$), $N_{\mathrm{init}}$ is the number of independent
optimizations with different initial guess (here $N_{\mathrm{init}}=5$),
$\bar{s}_{i}^{m}$ is the average spread of the WFs of material $m$ with
initialization $i$, and $\mathrm{Var}_{\mathrm{init}}$ is the variance of the
$N_{\mathrm{init}}$ independent runs. In the same way, we can define
$\langle\langle s_{\mathrm{max}}\rangle_{\mathrm{init}}\rangle_{\mathrm{mat}}$
and $\langle\sigma_{\mathrm{init}}(s_{\mathrm{max}})\rangle_{\mathrm{mat}}$ by
replacing $\bar{s}$ with $s_{\mathrm{max}}$ in the above equations. Note, that
$\langle\sigma_{\mathrm{init}}(\bar{s})\rangle_{\mathrm{mat}}$ and
$\langle\sigma_{\mathrm{init}}(s_{\mathrm{max}})\rangle_{\mathrm{mat}}$
measure the variation in the average and maximal spread in dependence of the
initial guess, i.e. the _robustness_ of the Wannierisation. In particular, it
does not reflect the variation of the spreads within a given set of WFs. As
demonstrated on several places in the paper, the latter is always
significantly and consistently reduced when using $\Omega_{\mathrm{var}}$.

Table 2 shows the results for the four quantities
$\langle\langle\bar{s}\rangle_{\mathrm{init}}\rangle_{\mathrm{mat}}$,
$\langle\langle
s_{\mathrm{max}}\rangle_{\mathrm{init}}\rangle_{\mathrm{mat}}$,
$\langle\sigma_{\mathrm{init}}(\bar{s})\rangle_{\mathrm{mat}}$, and
$\langle\sigma_{\mathrm{init}}(s_{\mathrm{max}})\rangle_{\mathrm{mat}}$ for
each of the spread functionals $\Omega$ and $\Omega_{\mathrm{var}}$. These
numbers summarise the information in Fig. 7. On this basis we conclude that
the penalization of the spread variance as done in the
$\mathrm{\Omega_{var}}$-functional, does not affect the average spread of the
WFs, but consistently decreases the maximum spread and improves the robustness
with respect to the choice of initial orbitals.

###  IV.5 Data availability

The data presented in this article together with information required for its
reproduction, are available in an online open access archive [26]. For more
details see also App. F.

##  V Conclusions

We have introduced a new localization functional for generating Wannier
functions (WFs) with balanced spread distributions. By penalizing the variance
of the spread distribution, the algorithm becomes less prone to produce
individual WFs with large spreads. It thereby resolves a well known problem of
standard Wannierisation schemes, which is particularly important when applied
to complex systems. Furthermore, we have proposed a general and fully
automatic algorithm for selecting the optimal number of WFs and the initial
set of orbitals for the Wannierisation procedure. Application to an extensive
test suite comprising both bulk and monolayer materials, point defects and
atoms on metal slabs, consistently show that these algorithms comprise a
highly robust approach for generating WFs, which should be particularly useful
for applications to complex systems and/or high-throughput studies. The
methods are implemented in Python and are available as part of the open-source
Atomic Simulation Environment (ASE) [9].

##  VI Acknowledgments

The Center for Nanostructured Graphene (CNG) is sponsored by the Danish
National Research Foundation, Project DNRF103. This project has received
funding from the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation program grant agreement No 773122 (LIMA).
K. S. T. is a Villum Investigator supported by VILLUM FONDEN (grant no.
37789).

##  Appendix A Wannier function spread

We define

$s_{n}=-\frac{1}{(2\pi)^{2}}\sum_{\alpha=1}^{N_{G}}W_{\alpha}\log\left(|Z_{\alpha,nn}|^{2}\right).$
(12)

This quantity is an approximation for the spread of a single WF in
$\mathrm{\SIUnitSymbolAngstrom}^{2}$ [11], if the weights $W_{\alpha}$ retain
the physical dimensions. This definition has the feature of exponentially
increasing in case of completely delocalized WFs, which can help with the
identification of localization issues but at the same time it does not retain
any physical meaning for delocalized WFs. We note that most definitions of
spread or localization functional have a convergence with the $k$-point grid
density, as observed and studied in Ref. [27]. In the present work we often
refer to $\bar{s}$, the average of $s_{n}$ over the set of WFs, and $s_{max}$,
the maximum value of $s_{n}$ over the set of WFs.

##  Appendix B Band interpolation error

In order to quantitatively estimate the quality of the band interpolation we
introduce a quantity called band interpolation error, that measures the
difference between the Kohn-Sham and the Wannier-interpolated band structures.
The definition is the same that appears Refs. [16, 6] and allows to easily
compare the results. The band interpolation error is computed up to the energy
threshold $E_{0}$, introduced in Section II.1.

If the band structure has a gap at the energy threshold, then the
interpolation error is computed only on the energy bands below the gap. We
define the average band interpolation error $\eta$ and the maximum
contribution to the band interpolation error $\eta_{\mathrm{max}}$ as

$\eta=\sqrt{\sum_{n\mathbf{k}}\frac{(\varepsilon^{\mathrm{KS}}_{n\mathbf{k}}-\varepsilon^{\mathrm{Wan}}_{n\mathbf{k}})^{2}}{N_{b}N_{\mathbf{k}}}}$
(13)
$\eta_{\mathrm{max}}=\text{max}_{n\mathbf{k}}(|\varepsilon^{\mathrm{KS}}_{n\mathbf{k}}-\varepsilon^{\mathrm{Wan}}_{n\mathbf{k}}|).$
(14)

where $\varepsilon^{KS}$ and $\varepsilon^{Wan}$ are the Kohn-Sham eigenvalues
and their Wannier interpolation respectively, $N_{b}$ is the number of bands
and $N_{\mathbf{k}}$ the number of $k$ points.

In the case of entangled bands we introduce a different definition which uses
a smearing function to weight the contributions from the bands around the
energy threshold

$\eta=\sqrt{\frac{\sum_{n\mathbf{k}}(\varepsilon^{\mathrm{KS}}_{n\mathbf{k}}-\varepsilon^{\mathrm{Wan}}_{n\mathbf{k}})^{2}\tilde{f}_{n\mathbf{k}}}{\sum_{n\mathbf{k}}\tilde{f}_{n\mathbf{k}}}}$
(15)
$\eta_{\mathrm{max}}=\text{max}_{n\mathbf{k}}(\tilde{f}_{n\mathbf{k}}|\varepsilon^{\mathrm{KS}}_{n\mathbf{k}}-\varepsilon^{\mathrm{Wan}}_{n\mathbf{k}}|)$
(16)

where
$\tilde{f}_{n\mathbf{k}}=\sqrt{f^{\mathrm{KS}}_{n\mathbf{k}}(\nu,\tau)f^{\mathrm{Wan}}_{n\mathbf{k}}(\nu,\tau)}$
and $f_{n\mathbf{k}}(\nu,\tau)$ is the Fermi-Dirac distribution for the state
at energy $\varepsilon_{n\mathbf{k}}$, $\nu$ is a fictitious chemical
potential fixed at $E_{0}$ and $\tau$ is a smearing width fixed to
$0.1\text{\,}\mathrm{eV}$.

##  Appendix C Computational details

For each material in this work we performed a self-consistent PBE calculation
with the GPAW code[28] using a Monkhorst-Pack grid [29] with a minimum density
of 5 $k$-points per ${\mathrm{\SIUnitSymbolAngstrom}}^{-1}$ and a real-space
grid with a spacing of $0.2\text{\,}\mathrm{\SIUnitSymbolAngstrom}$. Several
unoccupied states were included in the calculation. The band structure
calculation, used to evaluate the band interpolation accuracy, was performed
by fixing the density to the value from the previous self-consistent
calculation and evaluating the band structure along a path with minimum
density of 50 $k$-points per ${\mathrm{\SIUnitSymbolAngstrom}}^{-1}$. Finally,
we computed the WFs with ASE [9], starting from the self-consistent DFT Bloch
states. For all materials in Sec. IV.4 we used the structure available in the
public C2DB database [7].

##  Appendix D Methods for complex systems

For the NV center in diamond we substituted one carbon atom with a nitrogen
atom in a 64 atom cubic unit cell of bulk diamond. We then proceeded in the
structure relaxation using GPAW [28], with a 2x2x2 Monkhorst-Pack grid of
$k$-points in the Brillouin zone (BZ), a plane waves basis set with an energy
cutoff of $400\text{\,}\mathrm{eV}$ and an initial charge of $-1$. The
structure was optimized with the LBFGS algorithm [30] and the force threshold
was set to
$0.05\text{\,}\mathrm{eV}\text{\,}{\mathrm{\SIUnitSymbolAngstrom}}^{-1}$.
After the structure was relaxed we ran a self-consistent and a non self-
consistent calculation, with a BZ sampling density of $2$ and $5$ $k$-points
per ${\mathrm{\SIUnitSymbolAngstrom}}^{-1}$ respectively. The self-consistent
calculation was converged for every state up to $4\text{\,}\mathrm{eV}$ from
the conduction band minimum (CBM). The rest of the methods were in line with
Section IV.4.

With the adsorption systems we followed a similar workflow with ASE [9]. We
created a 2x2x4 slab of ruthenium with HCP(0001) surface and a vacuum layer of
$7\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ along the $z$ direction. Then the
adsorbate was placed on the HCP site at
$1\text{\,}\mathrm{\SIUnitSymbolAngstrom}$,
$1.108\text{\,}\mathrm{\SIUnitSymbolAngstrom}$,
$1.257\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ for H, N, O respectively. The
structure was then relaxed, fixing the first two layers of Ru, using the
finite-difference method and the PBE [17] xc-functional in GPAW, again using
LBFGS as optimization algorithm and
$0.05\text{\,}\mathrm{eV}\text{\,}{\mathrm{\SIUnitSymbolAngstrom}}^{-1}$ as
forces threshold. After the relaxation, a self-consistent and a non self-
consistent calculations were ran on each system, following the workflow for
the NV center in diamond.

##  Appendix E List of 2D materials

The 2D structures in Table 3 are picked from the C2DB database [7]. The values
for the maximum spread are the mean and the standard deviation over 5
calculations with different random seeds for the initial guess. The standard
deviation for the maximum band interpolation error is computed with the same
method. For the maximum band interpolation and the number of WFs (the latter
varies only in the case of optimal $N_{w}$) we use the value from the most
localized set (lowest maximum spread) over 5 calculations, instead of the mean
value. In the last row we show the average value and the average standard
deviation for each column.

|  | Minimal $N_{w}$ / $\Omega$ | Optimal $N_{w}$ / $\mathrm{\Omega_{var}}$
---|---|---|---
Formula | Crystal type | $N_{w}$ |  $\langle s_{\mathrm{max}}\rangle$ [$\mathrm{\SIUnitSymbolAngstrom}^{2}$] |  $\eta_{\mathrm{max}}$ [$\mathrm{meV}$] | $N_{w}$ |  $\langle s_{\mathrm{max}}\rangle$ [$\mathrm{\SIUnitSymbolAngstrom}^{2}$] |  $\eta_{\mathrm{max}}$ [$\mathrm{meV}$]
Ag2Br2 | AB-129-bc | 26 | $2.3\pm 0.5$ | $11.7\pm 0.2$ | 26 | $1.7\pm 0.7$ | $11.6\pm 0.1$
Al2Cl2O2 | ABC-59-ab | 20 | $4.3\pm 0.2$ | $13\pm 3$ | 23 | $3.6\pm 0.2$ | $5.0\pm 0.7$
Au2S2 | AB-10-fgm | 19 | $3.21\pm 0.08$ | $24\pm 5$ | 22 | $2.6\pm 0.1$ | $5\pm 3$
BrClTi | ABC-59-ab | 18 | $2.43\pm 0.01$ | $5.1\pm 0.1$ | 18 | $2.44\pm 0.09$ | $5.2\pm 0.3$
Br2Hf2S2 | ABC-156-ac | 29 | $6.3\pm 0.3$ | $373\pm 6$ | 32 | $1.57\pm 0.07$ | $105.3\pm 0.4$
C2H2 | AB-164-d | 8 | $5.0\pm 0.5$ | $13.3\pm 0.4$ | 8 | $4.9\pm 0.7$ | $8\pm 2$
CF2Y2 | AB2C2-164-bd | 23 | $6.9\pm 0.5$ | $36\pm 1$ | 27 | $3.2\pm 0.2$ | $18\pm 4$
CH2O2V2 | AB2C2D2-164-bd | 31 | $5.3\pm 1.6$ | $111\pm 58$ | 35 | $3\pm 0.6$ | $254\pm 75$
Cl2Sc2Se2 | ABC-59-ab | 32 | $1.64\pm 0.05$ | $6.3\pm 0.1$ | 33 | $1.9\pm 0.7$ | $6.5\pm 0.6$
Cr2W2S8 | ABC4-28-bcd | 56 | $7\pm 11$ | $8\pm 48$ | 59 | $2\pm 1$ | $5.2\pm 0.9$
CSiH2 | ABC2-156-ab | 9 | $5.763\pm 0.001$ | $6\pm 2$ | 13 | $4.15\pm 0.02$ | $3.0\pm 0.6$
Hf2Cl4 | AB2-11-e | 30 | $4.74\pm 0.01$ | $62\pm 1$ | 34 | $2.67\pm 0.04$ | $30\pm 1$
HgI2 | AB2-115-dg | 17 | $1.806\pm 0.001$ | $2.9\pm 0.1$ | 17 | $1.792\pm 0.001$ | $2.9\pm 0.1$
I2O2Rh2 | ABC-59-ab | 32 | $1.63\pm 0.01$ | $4\pm 0.1$ | 35 | $2\pm 0.9$ | $3\pm 280$
In2S2 | AB-164-cd | 22 | $6.83\pm 0.01$ | $52\pm 2$ | 26 | $2.4\pm 0.2$ | $5\pm 2$
Ir2Br6 | AB3-162-dk | 40 | $1.303\pm 0.001$ | $5.6\pm 0.1$ | 41 | $1.32\pm 0.03$ | $5.1\pm 0.2$
ISSb | ABC-156-abc | 17 | $6\pm 6$ | $4.5\pm 74$ | 17 | $2.5\pm 0.8$ | $5\pm 2$
MoSe2 | AB2-187-bi | 17 | $1.63\pm 0.04$ | $6.6\pm 0.1$ | 17 | $3\pm 2$ | $7\pm 32$
N2O2Zr3 | A2B2C3-187-bghi | 37 | $2.7\pm 0.3$ | $9.5\pm 20$ | 39 | $1.8\pm 0.7$ | $9.2\pm 0.9$
Nb2I4 | AB2-11-e | 35 | $3.12\pm 0.01$ | $100.1\pm 0.2$ | 38 | $3\pm 0.3$ | $97\pm 4$
PdS2 | AB2-164-bd | 13 | $1.26\pm 0.02$ | $2\pm 0.1$ | 13 | $2.5\pm 0.7$ | $2.0\pm 0.6$
PtSe2 | AB2-164-bd | 16 | $1.652\pm 0.001$ | $13.9\pm 0.1$ | 18 | $2.3\pm 0.2$ | $7\pm 5$
Ru2Se4 | AB2-11-e | 34 | $1.55\pm 0.01$ | $9.2\pm 0.4$ | 34 | $2.2\pm 0.7$ | $8.0\pm 0.5$
ScSe2 | AB2-164-bd | 12 | $4\pm 5$ | $8\pm 189$ | 16 | $1.5\pm 0.01$ | $1.8\pm 0.4$
SnTe2 | AB2-164-bd | 14 | $3.2\pm 0.4$ | $16.6\pm 0.1$ | 17 | $2.5\pm 0.9$ | $15\pm 3$
SrCl2 | AB2-164-bd | 17 | $6\pm 0.8$ | $0.8\pm 0.3$ | 21 | $4\pm 0.2$ | $1.2\pm 0.3$
TaS2 | AB2-187-bi | 13 | $5\pm 2$ | $16\pm 470$ | 17 | $1.6\pm 0.6$ | $3.1\pm 0.1$
WMo3Te8 | AB3C8-1-a | 64 | $21.51\pm 0.06$ | $420\pm 50$ | 68 | $2.12\pm 0.05$ | $10.6\pm 0.5$
WO2 | AB2-187-bi | 18 | $1.88\pm 0.01$ | $1.8\pm 0.1$ | 18 | $1.835\pm 0.002$ | $1.8\pm 0.1$
ZrTi3Te8 | AB3C8-1-a | 58 | $18.77\pm 0.03$ | $52\pm 4$ | 62 | $2.29\pm 0.06$ | $25\pm 7$
Average | - | - | $4.8\pm 0.8$ | $47\pm 31$ | - | $2.48\pm 0.15$ | $22\pm 14$
Table 3: List of the 2D materials considered in Sec. IV.4. The ”Crystal type”
in the second column follows the convention of the C2DB database and stands
for stoichiometry-space group-occupied Wyckoff positions.

##  Appendix F Software availability

Most of the software used and developed in this project is open-source and
available for free. The density functional theory code used is GPAW 20.10 [28]
with version 0.9.2 of the atomic setups, together with ASE 3.20 [9]. The
latter also includes the Wannier module in which the new spread functional has
been implemented. The Wannier functions are represented with Vesta [31], which
is available for free. The plots are produced with Matplotlib [32] and the
independent calculations were run in parallel with the help of GNU parallel
[33], both of them are open-source software.

## References

  * Wannier [1937] G. H. Wannier, Phys. Rev. 52, 191 (1937).
  * Marzari and Vanderbilt [1997] N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 (1997).
  * Thygesen  _et al._ [2005a] K. S. Thygesen, L. B. Hansen, and K. W. Jacobsen, Phys. Rev. Lett. 94, 026405 (2005a).
  * Souza  _et al._ [2001] I. Souza, N. Marzari, and D. Vanderbilt, Phys. Rev. B 65, 035109 (2001).
  * Thygesen  _et al._ [2005b] K. S. Thygesen, L. B. Hansen, and K. W. Jacobsen, Phys. Rev. B 72, 125119 (2005b).
  * Damle  _et al._ [2019] A. Damle, A. Levitt, and L. Lin, Multiscale Modeling & Simulation 17, 167 (2019).
  * Haastrup  _et al._ [2018] S. Haastrup, M. Strange, M. Pandey, T. Deilmann, P. S. Schmidt, N. F. Hinsche, M. N. Gjerding, D. Torelli, P. M. Larsen, A. C. Riis-Jensen, J. Gath, K. W. Jacobsen, J. J. Mortensen, T. Olsen, and K. S. Thygesen, 2D Materials 5, 042002 (2018).
  * King-Smith and Vanderbilt [1993] R. D. King-Smith and D. Vanderbilt, Phys. Rev. B 47, 1651 (1993).
  * Larsen  _et al._ [2017] A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dułak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. B. Jensen, J. Kermode, J. R. Kitchin, E. L. Kolsbjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. B. Maronsson, T. Maxson, T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard, J. Schiøtz, O. Schütt, M. Strange, K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Walter, Z. Zeng, and K. W. Jacobsen, Journal of Physics: Condensed Matter 29, 273002 (2017).
  * Resta and Sorella [1999] R. Resta and S. Sorella, Phys. Rev. Lett. 82, 370 (1999).
  * Berghold  _et al._ [2000] G. Berghold, C. J. Mundy, A. H. Romero, J. Hutter, and M. Parrinello, Phys. Rev. B 61, 10040 (2000).
  * Silvestrelli [1999] P. L. Silvestrelli, Phys. Rev. B 59, 9703 (1999).
  * Foster and Boys [1960] J. M. Foster and S. F. Boys, Rev. Mod. Phys. 32, 300 (1960).
  * Marzari  _et al._ [2012] N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, Rev. Mod. Phys. 84, 1419 (2012).
  * Damle and Lin [2018] A. Damle and L. Lin, Multiscale Modeling & Simulation 16, 1392 (2018).
  * Vitale  _et al._ [2020] V. Vitale, G. Pizzi, A. Marrazzo, J. Yates, N. Marzari, and A. Mostofi, npj Computational Materials 6, 66 (2020).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Gjerding  _et al._ [2021] M. N. Gjerding, A. Taghizadeh, A. Rasmussen, S. Ali, F. Bertoldo, T. Deilmann, U. P. Holguin, N. R. Knøsgaard, M. Kruse, S. Manti, T. G. Pedersen, T. Skovhus, M. K. Svendsen, J. J. Mortensen, T. Olsen, and K. S. Thygesen, arxiv.org/abs/2102.03029 , 1 (2021), arXiv:2102.03029 .
  * Castelli  _et al._ [2019] I. E. Castelli, T. Olsen, and Y. Chen, Journal of Physics: Energy 2, 011001 (2019).
  * Petralanda  _et al._ [2020] U. Petralanda, M. Kruse, H. Simons, and T. Olsen, arxiv.org/abs/2012.11254 , 1 (2020), arXiv:2012.11254 .
  * Giustino  _et al._ [2007] F. Giustino, M. L. Cohen, and S. G. Louie, Physical Review B 76, 165108 (2007).
  * Sjakste  _et al._ [2015] J. Sjakste, N. Vast, M. Calandra, and F. Mauri, Physical Review B 92, 054307 (2015).
  * Wang  _et al._ [2006] X. Wang, J. R. Yates, I. Souza, and D. Vanderbilt, Physical Review B 74, 195118 (2006).
  * Calzolari  _et al._ [2004] A. Calzolari, N. Marzari, I. Souza, and M. B. Nardelli, Physical Review B 69, 035108 (2004).
  * Strange  _et al._ [2008] M. Strange, I. S. Kristensen, K. S. Thygesen, and K. W. Jacobsen, The Journal of Chemical Physics 128, 114714 (2008).
  * Fontana and Thygesen [2021] P. F. Fontana and K. S. Thygesen, Robust wannier data (2021), 10.5281/zenodo.5338784.
  * Stengel and Spaldin [2006] M. Stengel and N. A. Spaldin, Phys. Rev. B 73, 075121 (2006).
  * Enkovaara  _et al._ [2010] J. Enkovaara, C. Rostgaard, J. J. Mortensen, J. Chen, M. Dułak, L. Ferrighi, J. Gavnholt, C. Glinsvad, V. Haikola, H. A. Hansen, H. H. Kristoffersen, M. Kuisma, A. H. Larsen, L. Lehtovaara, M. Ljungberg, O. Lopez-Acevedo, P. G. Moses, J. Ojanen, T. Olsen, V. Petzold, N. A. Romero, J. Stausholm-Møller, M. Strange, G. A. Tritsaris, M. Vanin, M. Walter, B. Hammer, H. Häkkinen, G. K. H. Madsen, R. M. Nieminen, J. K. Nørskov, M. Puska, T. T. Rantala, J. Schiøtz, K. S. Thygesen, and K. W. Jacobsen, Journal of Physics: Condensed Matter 22, 253202 (2010).
  * Monkhorst and Pack [1976] H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
  * Byrd  _et al._ [1995] R. H. Byrd, P. Lu, J. Nocedal, and C. Zhu, SIAM Journal on Scientific Computing 16, 1190 (1995).
  * Momma and Izumi [2011] K. Momma and F. Izumi, Journal of Applied Crystallography 44, 1272 (2011).
  * Hunter [2007] J. D. Hunter, Computing in Science & Engineering 9, 90 (2007).
  * Tange [2011] O. Tange, ;login: The USENIX Magazine 36, 42 (2011).