# Dissipation bounds precision of current response to kinetic perturbations

Krzysztof Ptaszyński<EMAIL_ADDRESS>Department of
Physics and Materials Science, University of Luxembourg, L-1511 Luxembourg
City, Luxembourg  Institute of Molecular Physics, Polish Academy of Sciences,
Mariana Smoluchowskiego 17, 60-179 Poznań, Poland  Timur Aslyamov
<EMAIL_ADDRESS>Department of Physics and Materials Science, University
of Luxembourg, L-1511 Luxembourg City, Luxembourg  Massimiliano Esposito
<EMAIL_ADDRESS>Department of Physics and Materials Science,
University of Luxembourg, L-1511 Luxembourg City, Luxembourg

###### Abstract

The precision of currents in Markov networks is bounded by dissipation via the
so-called thermodynamic uncertainty relation (TUR). In our work, we
demonstrate a similar inequality that bounds the precision of the static
current response to perturbations of kinetic barriers. Perturbations of such
type, which affect only the system kinetics but not the thermodynamic forces,
are highly important in biochemistry and nanoelectronics. We prove that our
inequality cannot be derived from the standard TUR. Instead, it implies the
standard TUR and provides an even tighter bound for dissipation. We also
provide a procedure for obtaining the optimal response precision for a given
model.

Introduction.—Among the most fundamental results of statistical physics are
the relations between the system response to external perturbations and
stationary thermodynamic observables. Close to equilibrium, such a relation is
given by the seminal fluctuation-dissipation theorem (FDT) linking the linear
response to external forces and equilibrium fluctuations [1]. The link between
dissipation and fluctuations has recently been generalized to a far-from-
equilibrium regime [2, 3, 4, 5, 6, 7, 8, 9, 10, 11]. In particular, for Markov
jump processes, the theory of stochastic thermodynamics [12] gave rise to the
thermodynamics uncertainty relation (TUR) linking the entropy production rate
$\dot{\sigma}$ which measures dissipation to the average of any current
$\mathcal{J}$ and its variance $\langle\langle\mathcal{J}\rangle\rangle$ as
[13, 14, 15, 16, 17, 18, 19, 20]

$\displaystyle\frac{\mathcal{J}^{2}}{\langle\langle\mathcal{J}\rangle\rangle}\leq\frac{\dot{\sigma}}{2}\,.$
(1)

This relation shows that to reduce the fluctuations, one needs to pay by
increased dissipation. In other words, it describes the thermodynamic cost of
current precision $\mathcal{J}^{2}/\langle\langle\mathcal{J}\rangle\rangle$.
It can also be used to infer the value of entropy production, which is often
directly inaccessible, using measurable currents and their fluctuations [21].

In our work, we get closer to the original formulation of FDT by providing a
thermodynamic trade-off between the current fluctuations and the static linear
current response, rather than the current itself. Recently, static responses
for Markov jump processes and chemical reaction networks have attracted
notable attention [22, 23, 24, 25, 26, 27, 28, 29, 30] due to their importance
in biophysical applications (proofreading, sensing [31, 22]). Specifically, we
focus on the response to kinetic perturbations, i.e., perturbations affecting
only the kinetics of transitions between the system states but not the
thermodynamic forces. Perturbations of this type are crucial in chemistry
(including biochemistry), as they correspond to the control of kinetic rates
by changing the concentration of a catalyst (e.g., enzyme) [32]. Among others,
kinetics may play a crucial role in determining the direction of motion of
molecular motors [33] and providing the driving mechanism of Maxwell demons
[34]. Chemical kinetics can also be controlled by magnetic fields, e.g., via
the radical pair mechanism [35], which is hypothesized to be the basis for
magnetoreception and other magnetic field effects in biology [36, 37, 38, 39].
Another important example of kinetic perturbations appears in the field of
nanoelectronics, where it corresponds to the adjustment of tunnel barriers in
quantum dots [40, 41, 42] or potential barriers in CMOS devices [43, 44] by
gate voltages. Crucially, the response of currents to kinetic perturbations
vanishes at equilibrium, where the state of the system depends only on its
thermodynamics, and not on its kinetics. This leads to the intuitive
expectation that the response to kinetic perturbations requires some
thermodynamic cost. Indeed, certain thermodynamics bounds on the response to
kinetic perturbations have been obtained for Markov jump processes [22, 24],
chemical reaction networks [25], and nonequilibrium diffusion [45, 46].
However, they apply to static system observables and their fluctuations rather
than to currents.

In this Letter, we provide numerical evidence to conjecture a new bound on the
current response to kinetic perturbations. We call it the response
thermodynamic uncertainty relation (R-TUR). We can prove it in special cases
detailed below, e.g., for unicyclic networks, linear response regime close to
equilibrium, and for single edge perturbations. To formulate that bound, we
parameterize the transition rates of the opposite Markov jumps $\pm e$ as

$\displaystyle W_{\pm e}=\exp\big{[}B_{e}(\varepsilon)\pm S_{e}/2\big{]}\,,$
(2)

where $B_{e}$ and $S_{e}$ parameterize the symmetric and asymmetric part of
the transition rate, respectively. The former term characterizes the kinetic
barriers, while the latter is the entropy change in the reservoir due to the
jump $e$. Only the kinetic part $B_{e}$ is assumed to depend on the control
parameter $\varepsilon$ (which can correspond, e.g., to the enzyme
concentration [32] or the gate voltage [40, 41, 42, 43, 44]). The bound reads

$\displaystyle\frac{(d_{\varepsilon}\mathcal{J})^{2}}{\langle\langle\mathcal{J}\rangle\rangle}\leq\frac{b_{\text{max}}^{2}\dot{\sigma}}{2}\,,$
(3)

where $d_{\varepsilon}\mathcal{J}\equiv d\mathcal{J}/d\varepsilon$ is the
static response of any current $\mathcal{J}$ to the parameter $\varepsilon$
and $b_{\text{max}}=\max_{e}|\partial_{\varepsilon}B_{e}|$ is the maximal rate
of change of kinetic barriers in the system. If we perturb a single barrier
corresponding to a given edge $e$, our bound (3) simplifies to

$\displaystyle\frac{(d_{B_{e}}\mathcal{J})^{2}}{\langle\langle\mathcal{J}\rangle\rangle}\leq\frac{\dot{\sigma}}{2}\,.$
(4)

Although the above bounds look similar to TUR (1), we prove that they cannot
be derived from TUR. Instead, TUR is a consequence of our result. Moreover, we
develop the optimization procedure that transforms R-TUR into a bound for
entropy production that is tighter than the optimized TUR [20]. However, in
some specific cases R-TUR can be derived from TUR. In particular, Eq. 3 can be
proven for unicyclic networks and the linear response regime close to
equilibrium. Furthermore, Eq. 4 can be derived when the current response at
the perturbed edge is considered, i.e., when $\mathcal{J}=j_{e}$.

Framework.—We consider a continuous-time Markov process describing stochastic
jumps between $N$ discrete states of the system (corresponding, e.g., to
chemical configurations, charge states of nanoelectronic devices, or molecule
conformations). It is described by the directed graph $\mathcal{G}$, where the
nodes of $\mathcal{G}$ correspond to the system states and the edges $\pm
e\in\mathcal{E}$ to transitions between states. We focus on the steady state
$\bm{\pi}$ given by the condition

$\displaystyle\mathbb{W}\cdot\bm{\pi}=0\,,$ (5)

where $\bm{\pi}=(\dots,\pi_{i},\dots)^{\intercal}$ is the vector of state
probabilities $\pi_{i}$ (with $\sum_{i}\pi_{i}=1$). $\mathbb{W}$ is here the
rate matrix with non-diagonal elements $W_{t(\pm e)s(\pm e)}=W_{\pm
e}(\varepsilon)$ describing the transition rate from the state $s(\pm e)$ (the
source of the edge $\pm e$) to the state $t(\pm e)$ (the tip of the edge $\pm
e$) and with the diagonal elements $W_{ii}=-\sum_{j\neq i}W_{ji}$.

Figure 1: (a): The colored dots are $m(\bm{\beta},\bm{x})$ numerically
calculated for random values of the vectors $\bm{\beta}$ and $\bm{x}$. (b):
The colored dots are the optimized $m(\bm{\beta},\bm{x}^{*})$ with
$\bm{x}^{*}$ from Eq. 9. The dashed lines in (a,b) denote the bound (3).
Simulations: we generate graphs with random topology with $N=5$ nodes and
various numbers of edges (from $6$ to $10$ highlighted by colors from purple
to yellow). We calculate Eq. 2 with randomly distributed $-2\leq B_{e}\leq 2$
and $0\leq\exp(S_{e}/2)\leq 10$. For the vectors $\bm{\beta}$ and $\bm{x}$ we
use $-1\leq\beta_{e}\leq 1$ and $-2\leq x_{e}\leq 2$. For the simulation
methods, see Appendix A

The bound (3) can be applied to an arbitrary current observable, namely,
$\mathcal{J}=\sum_{e}x_{e}j_{e}$, where
$j_{e}=W_{+e}\pi_{s(+e)}-W_{-e}\pi_{s(-e)}$ are the edge currents and
$\bm{x}=(\dots,x_{e},\dots)^{\intercal}$ is an arbitrary vector. Using the
chain rule for the static response of $\mathcal{J}$ we arrive at

$\displaystyle d_{\varepsilon}\mathcal{J}$
$\displaystyle=\sum_{e}x_{e}d_{\varepsilon}j_{e}=\sum_{e}x_{e}\sum_{e^{\prime}}(d_{B_{e^{\prime}}}j_{e})\partial_{\varepsilon}B_{e^{\prime}}$
$\displaystyle=\sum_{e^{\prime}}b_{e^{\prime}}d_{B_{e^{\prime}}}\mathcal{J}=\bm{b}^{\intercal}\nabla\mathcal{J}\,,$
(6)

where we introduce
$\bm{b}=(\dots,\partial_{\varepsilon}B_{e},\dots)^{\intercal}$ known from the
model parameterization and
$\nabla\mathcal{J}=(\dots,d_{B_{e}}\mathcal{J},\dots)^{\intercal}$ where
$d_{B_{e^{\prime}}}\mathcal{J}=\sum_{e}x_{e}d_{B_{e^{\prime}}}j_{e}$ is the
static response to the edge parameters $B_{e^{\prime}}$. The entropy
production rate can be calculated as $\dot{\sigma}=\sum_{e}j_{e}S_{e}$ [12].

Numerical evidence.—We start by providing numerical evidence for our bound.
Let us first rewrite Eq. 3 by expressing $d_{\varepsilon}\mathcal{J}$ via
Dissipation bounds precision of current response to kinetic perturbations and
using the expression
$\langle\langle\mathcal{J}\rangle\rangle=\bm{x}^{\intercal}\mathbb{C}\bm{x}$
[20], with the covariance matrix $\mathbb{C}$ of the edge currents defined in
Appendix A. We obtain

$\displaystyle\forall_{\bm{\beta},\bm{x}}:m(\bm{\beta},\bm{x})\equiv\frac{(d_{\varepsilon}\mathcal{J})^{2}}{\langle\langle\mathcal{J}\rangle\rangle
b_{\text{max}}^{2}}=\frac{\big{(}\bm{\beta}\nabla\mathcal{J}\big{)}^{2}}{\bm{x}^{\intercal}\mathbb{C}\bm{x}}$
$\displaystyle\leq\frac{\dot{\sigma}}{2}\,,$ (7)

where $\bm{\beta}=\bm{b}/b_{\text{max}}$ with $b_{\text{max}}=\max_{e}|b_{e}|$
and $|\beta_{e}|\leq 1$. From Fig. 1a one can see that the bound (7) with
random vectors $\bm{\beta}$ and $\bm{x}$ indeed holds, but is very loose.
Thus, we now look for $\bm{\beta}$ and $\bm{x}$ that maximize
$m(\bm{\beta},\bm{x})$ for a given model.

Figure 2: (a) Numerical results (color dots) for the optimized bound $\mu$.
Dashed lines show Eq. 8. (b) Numerical results (color dots) for the ratio
$\ell/\mu$. For the details of the simulations see the caption of Fig. 1.

Bound optimization.—We notice that Eq. 7 holds for an arbitrary linear
combination of $j_{e}$ and an arbitrary dependence on $\varepsilon$, which
implies that both $\bm{x}$ and $\bm{b}$ are arbitrary vectors. Therefore, for
a given graph $\mathcal{G}$ and rate matrix $\mathbb{W}$, the optimized
version of Eq. 7 reads

$\displaystyle\mu\equiv\max_{\bm{\beta}}\max_{\bm{x}}m(\bm{\beta},\bm{x})$
$\displaystyle\leq\frac{\dot{\sigma}}{2}.$ (8)

The problem $\max_{\bm{x}}m(\bm{\beta},\bm{x})$ is a standard convex
optimization of a quadratic form, and has been previously applied to TUR (1)
in Ref. [20]. The optimal vector $\bm{x}^{*}$ reads

$\displaystyle\bm{x}^{*}(\bm{\beta})=(d_{\varepsilon}\bm{j}^{\intercal}\mathbb{C}^{-1}d_{\varepsilon}\bm{j})^{-1}\mathbb{C}^{-1}d_{\varepsilon}\bm{j}\,,$
(9)

which results in $m(\bm{\beta},\bm{x}^{*})$ for the inner maximum in Eq. 8. To
account for possible zero eigenvalues, the matrix $\mathbb{C}^{-1}$ is here
the Drazin inverse [47], which for real symmetric (and generally Hermitian)
matrices is equivalent to Moore–Penrose pseudoinverse [48]. However, as shown
in Fig. 1b, the optimization Eq. 9 does not greatly improve the tightness of
the bound Eq. 7. To do this, we also need to optimize $\bm{\beta}$.

The problem $\max_{\bm{\beta}}m(\bm{\beta},\bm{x})$ for a fixed $\bm{x}$ is
solved simply by $\beta_{e}^{*}=\text{sign}(d_{B_{e}}\mathcal{J})$ due to the
constraint $|\beta_{e}|\leq 1$. Therefore, the optimal vector $\bm{\beta}^{*}$
is a certain combination of $\pm 1$ elements. Thus, joint optimization over
$\bm{\beta}$ and $\bm{x}$ can be performed by maximizing
$m(\bm{\beta},\bm{x}^{*})$ over all possible vectors $\bm{\beta}$ being
combinations of $\pm 1$ elements. As shown in Fig. 2a in this case the bound
(8) can become tight, especially close to equilibrium, i.e., for small values
of the entropy production.

Relation between R-TUR and TUR.—We now prove that the conjectured R-TUR (3)
cannot be derived from the standard TUR (1). Instead, it implies the validity
of TUR. To that end, we recall the formulation of R-TUR via Eq. 7 (which is
equivalent to Eq. 3) and show that the standard TUR is equivalent to a weaker
version of that bound,

$\displaystyle\forall_{\bm{x}}:\quad
m(\bm{1},\bm{x})=\frac{\big{(}\bm{1}\nabla\mathcal{J}\big{)}^{2}}{\langle\langle\mathcal{J}\rangle\rangle}$
$\displaystyle\leq\frac{\dot{\sigma}}{2}\,,$ (10)

where $\bm{\beta}=\bm{1}$ (vector of ones) corresponds to homogeneous
perturbation of the kinetic rates. This can be proven by noting that

$\displaystyle\bm{1}\nabla\mathcal{J}=\sum_{e}d_{B_{e}}\mathcal{J}=\sum_{e^{\prime}}x_{e^{\prime}}\sum_{e}d_{B_{e}}j_{e^{\prime}}=\sum_{e^{\prime}}x_{e^{\prime}}j_{e^{\prime}}=\mathcal{J}.$
(11)

Here, in the second-last step, we used the summation response relation from
Ref. [27] applied to Eq. 2:

$\displaystyle\sum_{e}d_{B_{e}}j_{e^{\prime}}=j_{e^{\prime}}\,.$ (12)

Inserting $\bm{1}\nabla\mathcal{J}=\mathcal{J}$ into Eq. 10 we obtain the
standard TUR (1). Since Eq. 10 is a weaker version of Eq. 7, R-TUR implies
TUR, while TUR does not necessarily imply R-TUR. Consequently, in the generic
case, R-TUR cannot be derived from TUR.

Quantitative comparison with TUR.—We now show that for most Markov networks
Eq. 10 is indeed a weaker condition than Eq. 7. Consequently, for a given
Markov network, the optimal response precision
$(d_{\varepsilon}\mathcal{J})^{2}/\langle\langle\mathcal{J}\rangle\rangle$ can
be greater than the optimal current precision
$\mathcal{J}^{2}/\langle\langle\mathcal{J}\rangle\rangle$. As shown in Ref.
[20], the latter can be maximized as

$\displaystyle\ell\equiv\max_{\bm{x}}\frac{\mathcal{J}^{2}}{\langle\langle\mathcal{J}\rangle\rangle}=m(\bm{1},\bm{x}^{*})\geq\frac{\mathcal{J}^{2}}{\langle\langle\mathcal{J}\rangle\rangle}\,,$
(13)

where Eq. 1 also implies $\ell\leq\dot{\sigma}/2$; here we use the above
proven relation
$m(\bm{1},\bm{x})=\mathcal{J}^{2}/\langle\langle\mathcal{J}\rangle\rangle$. We
thus obtain the chain of inequalities

$\displaystyle\frac{\mathcal{J}^{2}}{\langle\langle\mathcal{J}\rangle\rangle}\leq\ell\leq\mu\leq\frac{\dot{\sigma}}{2}\,,$
(14)

with the equality $\ell=\mu$ for $\bm{\beta}^{*}=\bm{1}$ (i.e., when the
optimal $\bm{\beta}$ corresponds to homogeneous perturbation). This implies
that the optimized R-TUR (8) provides a tighter bound for the entropy
production than the standard TUR. We demonstrate this in Fig. 2b. Indeed, for
certain models the ratio $\ell/\mu$ is close to 0, which implies that the
optimized R-TUR (8) is far tighter than the standard TUR. Only in less than
10% of the instances $\ell=\mu$, i.e., R-TUR (8) is optimized for homogeneous
perturbation, which makes it equivalent to the optimized TUR (13).

Proof for local responses.—We have proven that in general R-TUR cannot be
derived from TUR. However, this is possible for some specific cases. The first
example is the response of a single edge current $\mathcal{J}=j_{e}$ to the
perturbation of that edge. Applying the results of Ref. [27] to the rate
parameterization (2) (which implies $\partial_{B_{e}}j_{e}=j_{e}$), we arrive
at

$\displaystyle 0\leq j_{e}^{-1}d_{B_{e}}j_{e}\leq 1.$ (15)

In conjunction with TUR (1), this implies Eq. (4):

$\displaystyle\frac{(d_{B_{e}}j_{e})^{2}}{\langle\langle
j_{e}\rangle\rangle}\leq\frac{j_{e}^{2}}{\langle\langle
j_{e}\rangle\rangle}\leq\frac{\dot{\sigma}}{2}\,.$ (16)

Proof for unicyclic networks.—Eq. 3 can also be proven for unicyclic networks,
i.e., systems in which each state $i\in\\{0,\ldots,N-1\\}$ is a source of a
single directed edge pointing to the tip $i+1$ (with $i$ defined modulo $N$).
In such a case, all edge currents $j_{e}=j$ are equal to each other and
$\mathcal{J}=\sum_{e}x_{e}j_{e}=x_{\Sigma}j$, where
$x_{\Sigma}=\sum_{e}x_{e}$. Moreover, we have

$\displaystyle\text{sign}\,d_{B_{e}}\mathcal{J}=\text{sign}\,x_{\Sigma}d_{B_{e}}j=\text{sign}\,x_{\Sigma}j\,,$
(17)

where we use Eq. 15 for $d_{B_{e}}j/j\geq 0$. We can always define the
orientation of $\mathcal{G}$ so that $x_{\Sigma}j$ is positive, which implies
$\beta_{e}^{*}=\text{sign}(d_{B_{e}}\mathcal{J})=1$ (see Eq. 17), i.e., $\mu$
is maximized for homogeneous perturbation. Then, as discussed below Eq. 14,
$\mu=\ell\leq\dot{\sigma}/2$, which proves Eq. 3.

Proof for the linear response regime.—Finally, R-TUR (3) can be proved to hold
close to equilibrium, where the currents are linear in the applied
thermodynamic forces. To define this regime, let us parameterize the
asymmetric transition rate parameters in Eq. 2 as
$S_{e}=F_{e}+E_{s(e)}-E_{t(e)}$ [22, 49], where $s(e)$ and $t(e)$ are the
source and the tip of the edge $e$, $E_{i}$ are node parameters (which may
sometimes be interpreted as state energies), and $F_{e}$ are the
nonconservative thermodynamic forces. In the absence of forces, the system
relaxes to the equilibrium state $\pi_{i}^{\text{eq}}=e^{-E_{i}}/Z$, with
$Z=\sum_{i}e^{-E_{i}}$, which is independent of $B_{e}$. Close to equilibrium
(i.e., for small values of nonconservative forces) the currents respond
linearly to applied forces as $\bm{j}=\mathbb{L}\bm{F}$, where
$\bm{F}=(\ldots,F_{e},\ldots)^{T}$ and $\mathbb{L}$ is the positive
semidefinite Onsager matrix. Consequently, the current observables and the
entropy production can be calculated as
$\mathcal{J}=\bm{x}^{T}\mathbb{L}\bm{F}$ and
$\dot{\sigma}=\bm{F}^{T}\mathbb{L}\bm{F}$. Furthermore, in this regime, the
covariance matrix of edge currents is related to Onsager matrix as
$\mathbb{C}=2\mathbb{L}$ by virtue of the fluctuation-dissipation theorem
[50]. Let us also define the matrix $\mathbb{L}_{\bm{\beta}}\equiv
d_{\varepsilon}\mathbb{L}/b_{\text{max}}=\sum_{e}\beta_{e}d_{B_{e}}\mathbb{L}$,
which describes the kinetic response of Onsager matrix. Inserting the above
expressions for $\mathbb{C}$ and $\dot{\sigma}$ into Eq. 7, using
$\bm{\beta}\nabla\mathcal{J}=\sum_{e}\beta_{e}d_{B_{e}}\mathcal{J}=\bm{x}^{T}\mathbb{L}_{\bm{\beta}}\bm{F}$,
and applying optimization (9), we find that R-TUR (3) holds if and only if

$\displaystyle\forall_{\bm{F}}:\quad\bm{F}^{T}\mathbb{L}_{\bm{\beta}}\mathbb{L}^{-1}\mathbb{L}_{\bm{\beta}}\bm{F}\leq\bm{F}^{T}\mathbb{L}\bm{F}\,.$
(18)

We recall that $\mathbb{L}^{-1}$ here denotes the Moore-Penrose pseudoinverse
of $\mathbb{L}$. This is equivalent to

$\displaystyle\mathbb{L}-\mathbb{L}_{\bm{\beta}}\mathbb{L}^{-1}\mathbb{L}_{\bm{\beta}}\succcurlyeq
0\,,$ (19)

where $\mathbb{A}\succcurlyeq 0$ means that the matrix $\mathbb{A}$ is
positive semidefinite. To prove the above inequality, we note that
$\mathbb{L}-\mathbb{L}_{\bm{\beta}}\mathbb{L}^{-1}\mathbb{L}_{\bm{\beta}}$ is
the Shur complement $\mathbb{M}/\mathbb{L}$ of the matrix

$\displaystyle\mathbb{M}=\begin{pmatrix}\mathbb{L}&\mathbb{L}_{\bm{\beta}}\\\
\mathbb{L}_{\bm{\beta}}&\mathbb{L}\end{pmatrix}\,.$ (20)

Thus, Eq. 19 holds if $\mathbb{M}\succcurlyeq 0$. Since eigenvalues of
$\mathbb{M}$ correspond to that of $\mathbb{L}\pm\mathbb{L}_{\bm{\beta}}$
[51], $\mathbb{M}\succcurlyeq 0$ if and only if
$\mathbb{L}\pm\mathbb{L}_{\bm{\beta}}\succcurlyeq 0$. To prove the latter
occurrence, in Appendix B we show that the responses $d_{B_{e}}\mathbb{L}$ are
positive semidefinite and obey the summation response relation

$\displaystyle\sum_{e}d_{B_{e}}\mathbb{L}=\mathbb{L}\,.$ (21)

As a consequence,
$\mathbb{L}\pm\mathbb{L}_{\bm{\beta}}=\sum_{e}(1\pm\beta_{e})d_{B_{e}}\mathbb{L}$
is positive semidefinite because $1\pm\beta_{e}\in[0,2]\geq 0$ and
$d_{B_{e}}\mathbb{L}\succcurlyeq 0$. This proves Eq. 19, and thus R-TUR (3).

We also note that Eq. 19 is saturated for homogeneous perturbation
$\bm{\beta}=\bm{1}$. Thus, in the linear response regime, the optimized R-TUR
(8) is equivalent to the optimized TUR (13). (We note that, in contrast to the
unicyclic case, this does not imply that homogeneous perturbation will always
yield the optimal response precision when considering a generic rather than
optimized current observable.) Using Eq. 9, both are optimized for the current
observable corresponding to the entropy production:
$\mathcal{J}=\dot{\sigma}$. Indeed, in the linear response regime,
$\langle\langle\dot{\sigma}\rangle\rangle=2\dot{\sigma}$, and thus TUR (1)
becomes saturated:
$\dot{\sigma}^{2}/\langle\langle\dot{\sigma}\rangle\rangle=\dot{\sigma}/2$
[15]. Furthermore, using $d_{B_{e}}\mathbb{L}\succcurlyeq 0$, Eq. 21 and
$\dot{\sigma}=\bm{F}^{T}\mathbb{L}\bm{F}$, we find that close to equilibrium
the entropy production response itself, and not only its precision, is bounded
by dissipation as

$\displaystyle 0\leq
d_{B_{e}}\dot{\sigma}\leq\dot{\sigma}\quad\text{with}\quad\sum_{e}d_{B_{e}}\dot{\sigma}=\dot{\sigma}\,,$
(22) $\displaystyle|d_{\varepsilon}\dot{\sigma}|\leq
b_{\text{max}}\dot{\sigma}\,.$ (23)

Finally, we note that the Onsager matrix can be represented in the basis of
measurable currents (e.g., charge, heat, or chemical currents) conjugated to
fundamental thermodynamic forces (e.g., temperature or chemical potential
gradients) [52]. Such physical Onsager matrix can be expressed as
$\mathbb{L}_{\mathbb{X}}=\mathbb{X}\mathbb{L}\mathbb{X}^{T}$, where
$\mathbb{X}$ is the matrix that connects measurable currents to edge currents
[50]. Using $d_{B_{e}}\mathbb{L}\succcurlyeq 0$ and Eq. 21, its kinetic
response can also be bounded as

$\displaystyle\mathbb{L}_{\mathbb{X}}\succcurlyeq
d_{B_{e}}\mathbb{L}_{\mathbb{X}}\succcurlyeq
0\quad\text{with}\quad\sum_{e}d_{B_{e}}\mathbb{L}_{\mathbb{X}}=\mathbb{L}_{\mathbb{X}}\,,$
(24) $\displaystyle b_{\text{max}}\mathbb{L}_{\mathbb{X}}\succcurlyeq
d_{\varepsilon}\mathbb{L}_{\mathbb{X}}\succcurlyeq-
b_{\text{max}}\mathbb{L}_{\mathbb{X}}\,.$ (25)

This enables one to experimentally bound $b_{\text{max}}$ using only
measurable currents.

Final remarks.—Our bound establishes a fundamental thermodynamic cost for
controlling currents by modulating the system’s kinetics. It may prove useful
in optimizing the energetic cost needed to accurately control nanodevices such
as molecular machines or nanoelectronic devices. It may also be used for
thermodynamic [21] or model [53] inference. In particular, when all parameters
$b_{e}$ are known (e.g., when the control parameter $\varepsilon$ affects only
some known transition rates), the empirical current responses and fluctuations
may provide a better estimate of the entropy production than the standard TUR.
On the other hand, when the model details are unknown but the entropy
production is accessible, the inequality (3) can be used to bound the
transition rate responses $b_{e}$.

###### Acknowledgements.

K.P., T.A. and M.E. acknowledge the financial support from, repectively,
project No. 2023/51/D/ST3/01203 funded by the National Science Centre, Poland,
project ThermoElectroChem (C23/MS/18060819) from Fonds National de la
Recherche-FNR, Luxembourg, project TheCirco (INTER/FNRS/20/15074473) funded by
FRS-FNRS (Belgium) and FNR (Luxembourg).

##  Appendix A Numerical Simulations

To calculate the response $\nabla\mathcal{J}$ we use the results of Ref. [27]
as follows

$\displaystyle\nabla^{\intercal}\mathcal{J}=\Big{[}\sum_{e}x_{e}d_{B_{e^{\prime}}}j_{e}\Big{]}_{e^{\prime}}=\bm{x}^{\intercal}\mathbb{P}\mathbb{J}\,,$
(26)

where the matrix $\mathbb{P}$ is the projection matrix with analytical
expression [27] and
$\mathbb{J}=[\partial_{B_{e^{\prime}}}j_{e}]_{\\{e,e^{\prime}\\}}=\text{diag}(\dots,j_{e},\dots)$
the diagonal Jacobian matrix. Thus, to simulate Eq. 26, we only need the
numerical values of the steady state $\bm{\pi}$ satisfying Eq. 5.

To find the covariance matrix of the edge currents $\mathbb{C}$, we follow the
method from Ref. [54] and define the matrix $\mathbb{W}^{\phi}(\bm{q})$ with
nondiagonal elements $W^{\phi}_{t(\pm e)s(\pm e)}(\bm{q})=W_{t(\pm e)s(\pm
e)}\exp(\pm q_{e})$ and diagonal elements the same as $\mathbb{W}$. The
elements of the covariance matrix are defined as

$\displaystyle C_{ee^{\prime}}=\frac{\partial}{\partial
q_{e}}\frac{\partial}{\partial
q_{e^{\prime}}}\lambda(\bm{q})\Big{|}_{\bm{q}=\bm{0}}\,,$ (27)

where $\lambda$ is eigenvalue of the matrix $\mathbb{W}^{\phi}$ with the
largest real part. In simulations of Eq. 27, we used the method of finite
differences to numerically calculate $\lambda(\bm{q})$.

##  Appendix B Proof of $d_{B_{e}}\mathbb{L}\succcurlyeq 0$ and Eq. (21)

To prove $d_{B_{e}}\mathbb{L}\succcurlyeq 0$ we use the expression for Onsager
matrix from Ref. [50],

$\displaystyle\mathbb{L}=\tfrac{1}{2}\mathbb{K}\mathbb{H}^{-1}\mathbb{K}^{T}\quad\text{with}\quad\mathbb{H}=\mathbb{K}^{T}\mathbb{T}^{-1}\mathbb{H}\,.$
(28)

Here $\mathbb{K}$ is the cycle matrix that depends only on the system topology
and not on transition rates, while
$\mathbb{T}=\text{diag}(\ldots,\tau_{e}^{\text{eq}},\ldots)$ is the matrix
with equilibrium traffics
$\mathcal{\tau}_{e}^{\text{eq}}=W_{+e}\pi^{\text{eq}}_{s(+e)}+W_{-e}\pi^{\text{eq}}_{s(-e)}$
at the diagonal. We note that Ref. [50] actually considered the Onsager matrix
describing current responses to edge affinities $\mathcal{F}_{e}$ rather than
forces $F_{e}$ defined here; however, the Onsager matrix for both force
definitions is the same, as it is related to the covariance matrix via the
fluctuation-dissipation theorem [15]
$\langle\langle\dot{\sigma}\rangle\rangle=2\dot{\sigma}\Leftrightarrow\mathbb{L}=\mathbb{C}/2$.
We now use the expressions
$d_{B_{e}}\mathbb{H}^{-1}=-\mathbb{H}^{-1}(d_{B_{e}}\mathbb{H})\mathbb{H}^{-1}$,
and the equality
$d_{B_{e}}\tau_{e^{\prime}}^{\text{eq}}=\delta_{ee^{\prime}}\tau_{e}^{\text{eq}}$
that holds because probabilities $\pi_{i}^{\text{eq}}$ do not depend on
$B_{e}$; the latter yields $d_{B_{e}}\mathbb{T}^{-1}=-\mathbb{T}_{e}^{-1}$
with $\mathbb{T}_{e}=\text{diag}(\ldots,0,\tau_{e}^{\text{eq}},0,\ldots)$. We
get

$\displaystyle
d_{B_{e}}\mathbb{L}=\tfrac{1}{2}\mathbb{K}\mathbb{H}^{-1}\mathbb{K}^{T}\mathbb{T}_{e}^{-1}\mathbb{K}\mathbb{H}^{-1}\mathbb{K}^{T}=2\mathbb{L}\mathbb{T}_{e}^{-1}\mathbb{L}\,.$
(29)

This is a positive semidefinite matrix, as this is a sandwich of the positive
semidefinite diagonal matrix $\mathbb{T}_{e}^{-1}$ between the matrix
$\mathbb{L}$ and its transpose $\mathbb{L}^{T}=\mathbb{L}$, which concludes
the proof of $d_{B_{e}}\mathbb{L}\succcurlyeq 0$.

Summing Eq. 29 over edges, we further obtain Eq. 21:

$\displaystyle\sum_{e}d_{B_{e}}\mathbb{L}=\tfrac{1}{2}\mathbb{K}\mathbb{H}^{-1}\mathbb{H}\mathbb{H}^{-1}\mathbb{K}^{T}=\mathbb{L}\,.$
(30)

Here we use $\sum_{e}\mathbb{T}_{e}^{-1}=\mathbb{T}^{-1}$ and the relation
$\mathbb{H}^{-1}\mathbb{H}\mathbb{H}^{-1}=\mathbb{H}^{-1}$ that holds for
Moore-Penrose pseudoinverse.

## References

  * Kubo [1966] R. Kubo, Rep. Prog. Phys. 29, 255 (1966).
  * Agarwal [1972] G. S. Agarwal, Zeitschrift für Physik A Hadrons and nuclei 252, 25 (1972).
  * Seifert and Speck [2010] U. Seifert and T. Speck, Europhys. Lett. 89, 10007 (2010).
  * Prost  _et al._ [2009] J. Prost, J.-F. Joanny, and J. M. Parrondo, Phys. Rev. Lett. 103, 090601 (2009).
  * Altaner  _et al._ [2016] B. Altaner, M. Polettini, and M. Esposito, Phys. Rev. Lett. 117, 180601 (2016).
  * Chun  _et al._ [2021] H.-M. Chun, Q. Gao, and J. M. Horowitz, Phys. Rev. Res. 3, 043172 (2021).
  * Baiesi  _et al._ [2009] M. Baiesi, C. Maes, and B. Wynants, Phys. Rev. Lett. 103, 010602 (2009).
  * Dechant and Sasa [2020] A. Dechant and S.-i. Sasa, Proceedings of the National Academy of Sciences 117, 6430 (2020).
  * Di Terlizzi and Baiesi [2018] I. Di Terlizzi and M. Baiesi, Journal of Physics A: Mathematical and Theoretical 52, 02LT03 (2018).
  * Falasco  _et al._ [2022] G. Falasco, M. Esposito, and J.-C. Delvenne, Journal of Physics A: Mathematical and Theoretical 55, 124002 (2022).
  * Hasegawa [2024] Y. Hasegawa, Phys. Rev. Lett. 132, 087102 (2024).
  * Seifert [2012] U. Seifert, Rep. Prog. Phys. 75, 126001 (2012).
  * Barato and Seifert [2015] A. C. Barato and U. Seifert, Phys. Rev. Lett. 114, 158101 (2015).
  * Gingrich  _et al._ [2016] T. R. Gingrich, J. M. Horowitz, N. Perunov, and J. L. England, Phys. Rev. Lett. 116, 120601 (2016).
  * Pietzonka  _et al._ [2016] P. Pietzonka, A. C. Barato, and U. Seifert, Phys. Rev. E 93, 052145 (2016).
  * Pietzonka  _et al._ [2017] P. Pietzonka, F. Ritort, and U. Seifert, Phys. Rev. E 96, 012101 (2017).
  * Horowitz and Gingrich [2017] J. M. Horowitz and T. R. Gingrich, Phys. Rev. E 96, 020103 (2017).
  * Falasco  _et al._ [2020] G. Falasco, M. Esposito, and J.-C. Delvenne, New Journal of Physics 22, 053046 (2020).
  * Horowitz and Gingrich [2020] J. M. Horowitz and T. R. Gingrich, Nature Physics 16, 15 (2020).
  * Van Vu  _et al._ [2020] T. Van Vu, V. T. Vo, and Y. Hasegawa, Phys. Rev. E 101, 042138 (2020).
  * Seifert [2019] U. Seifert, Annual Review of Condensed Matter Physics 10, 171 (2019).
  * Owen  _et al._ [2020] J. A. Owen, T. R. Gingrich, and J. M. Horowitz, Phys. Rev. X 10, 011066 (2020).
  * Owen and Horowitz [2023] J. A. Owen and J. M. Horowitz, Nature Communications 14, 1280 (2023).
  * Fernandes Martins and Horowitz [2023] G. Fernandes Martins and J. M. Horowitz, Phys. Rev. E 108, 044113 (2023).
  * Chun and Horowitz [2023] H.-M. Chun and J. M. Horowitz, The Journal of Chemical Physics 158, 174115 (2023).
  * Aslyamov and Esposito [2024a] T. Aslyamov and M. Esposito, Phys. Rev. Lett. 132, 037101 (2024a).
  * Aslyamov and Esposito [2024b] T. Aslyamov and M. Esposito, arXiv preprint 10.48550/arXiv.2402.13990 (2024b).
  * Harunari  _et al._ [2024] P. E. Harunari, S. D. Cengio, V. Lecomte, and M. Polettini, arXiv preprint 10.48550/arXiv.2402.13193 (2024).
  * Floyd  _et al._ [2024] C. Floyd, A. R. Dinner, and S. Vaikuntanathan, arXiv preprint 10.48550/arXiv.2404.03798 (2024).
  * Zheng and Lu [2024] J. Zheng and Z. Lu, arXiv preprint 10.48550/arXiv.2403.10952 (2024).
  * Murugan  _et al._ [2014] A. Murugan, D. A. Huse, and S. Leibler, Phys. Rev. X 4, 021016 (2014).
  * Wachtel  _et al._ [2018] A. Wachtel, R. Rao, and M. Esposito, New Journal of Physics 20, 042002 (2018).
  * Maes and de Galway [2015] C. Maes and W. O. de Galway, Physica A: Statistical Mechanics and its Applications 436, 678 (2015).
  * Esposito and Schaller [2012] M. Esposito and G. Schaller, Europhys. Lett. 99, 30003 (2012).
  * Steiner and Ulrich [1989] U. E. Steiner and T. Ulrich, Chemical Reviews 89, 51 (1989).
  * Ritz  _et al._ [2000] T. Ritz, S. Adem, and K. Schulten, Biophys. J. 78, 707 (2000).
  * Hore and Mouritsen [2016] P. J. Hore and H. Mouritsen, Annu. Rev. Biophys. 45, 299 (2016).
  * Wiltschko and Wiltschko [2019] R. Wiltschko and W. Wiltschko, Journal of the Royal Society Interface 16, 20190295 (2019).
  * Zadeh-Haghighi and Simon [2022] H. Zadeh-Haghighi and C. Simon, Journal of the Royal Society Interface 19, 20220325 (2022).
  * Gustavsson  _et al._ [2006] S. Gustavsson, R. Leturcq, B. Simovič, R. Schleser, T. Ihn, P. Studerus, K. Ensslin, D. C. Driscoll, and A. C. Gossard, Phys. Rev. Lett. 96, 076605 (2006).
  * Gustavsson  _et al._ [2009] S. Gustavsson, R. Leturcq, M. Studer, I. Shorubalko, T. Ihn, K. Ensslin, D. Driscoll, and A. Gossard, Surface Science Reports 64, 191 (2009).
  * Sánchez  _et al._ [2019] R. Sánchez, P. Samuelsson, and P. P. Potts, Phys. Rev. Res. 1, 033066 (2019).
  * Freitas  _et al._ [2021] N. Freitas, J.-C. Delvenne, and M. Esposito, Phys. Rev. X 11, 031064 (2021).
  * Gopal  _et al._ [2022] A. Gopal, M. Esposito, and N. Freitas, Phys. Rev. B 106, 155303 (2022).
  * Gao  _et al._ [2022] Q. Gao, H.-M. Chun, and J. M. Horowitz, Phys. Rev. E 105, L012102 (2022).
  * Gao  _et al._ [2024] Q. Gao, H.-M. Chun, and J. M. Horowitz, Europhys. Lett. 146, 31001 (2024).
  * Stanimirović  _et al._ [2014] P. S. Stanimirović, D. Pappas, and S. Miljković, Linear and Multilinear Algebra 62, 252 (2014).
  * Landi  _et al._ [2024] G. T. Landi, M. J. Kewming, M. T. Mitchison, and P. P. Potts, PRX Quantum 5, 020201 (2024).
  * Rao and Esposito [2018] R. Rao and M. Esposito, New Journal of Physics 20, 023007 (2018).
  * Vroylandt  _et al._ [2019] H. Vroylandt, D. Lacoste, and G. Verley, J. Stat. Mech.: Theory Exp. 2019 (5), 054002.
  * Nath and Paul [2014] M. Nath and S. Paul, Linear Algebra and its Applications 460, 97 (2014).
  * Forastiere  _et al._ [2022] D. Forastiere, R. Rao, and M. Esposito, New Journal of Physics 24, 083021 (2022).
  * Moffitt  _et al._ [2010] J. R. Moffitt, Y. R. Chemla, and C. Bustamante, in  _Methods Enzymol._, Vol. 475 (Elsevier, 2010) p. 221.
  * Wachtel  _et al._ [2015] A. Wachtel, J. Vollmer, and B. Altaner, Phys. Rev. E 92, 042132 (2015).


# Wannier function localisation
using Bloch intrinsic atomic orbitals

Andrew Zhu and David P. Tew  Physical & Theoretical Chemistry Laboratory
University of Oxford

###### Abstract

We extend the Intrinsic Atomic Orbital (IAO) method for localisation of
molecular orbitals to calculate well-localised generalised Wannier functions
in crystals using the Pipek–Mezey locality metric. We furthermore present a
one-shot diabatic Wannierisation procedure that aligns the phases of the Bloch
functions, providing immediate Wannier localisation, which serves as an
excellent initial guess for optimisation. We test our Wannier localisation
implementation on a number of solid state systems, highlighting the
effectiveness of the diabatic preparation, especially for localising core
bands. Partial charges of Wannier functions generated using Bloch IAOs align
well with chemical intuition, which we demonstrate through the example of
adsorption of CO on a MgO surface.

††preprint: APS/123-QED

##  I Introduction

Mean field theories, such as Hartree-Fock (HF) or Kohn-Sham density functional
theory (DFT)[1, 2], provide a description of the electronic structure of a
system through an one-particle orbital model, which enables interpretation of
the bonding in molecules and band structure for materials. However, the
canonical molecular orbitals (MOs) or Bloch functions, are typically
delocalised across the entire system, and thus do not intuitively map to the
interpretation of bonding in terms of overlap of atomic orbitals (AOs), which
is a local picture. By applying unitary rotations to the occupied orbitals,
one can obtain localised objects, commonly known as Wannier functions (WFs)[3]
for periodic systems. Localisation of occupied orbitals aid in interpretation
of electronic structure, and also provide a basis for reduced scaling quantum
chemistry methods, which exploit this locality to truncate the virtual
space[4, 5, 6, 7].

Methods to evaluate localised molecular orbitals have focused upon defining a
localisation metric or functional; the stationary points of this functional
thus correspond to localised orbitals. The two most commonly employed metrics
are from Foster and Boys (FB)[8, 9], and Pipek and Mezey (PM)[10], both of
which have been adapted for periodic systems[11, 12, 13, 14]. The FB metric,
which minimises the spread of the orbitals, has seen widespread usage, namely
through the Wannier90[12] package, which has now established interfaces with
various periodic, plane-wave based, codes[15, 16, 17, 18]. In contrast, the PM
metric, which uses the Mulliken partial atomic charges, is naturally suited to
codes employing localised basis sets, under a linear-combination-of-atomic-
orbitals (LCAO) framework[19, 20, 21, 22, 23], where AO coefficients are
directly accessible and overlaps are easily computed. In addition, WFs
localised with the PM metric produces orbitals with separate $\sigma$ and
$\pi$ bonding character, giving advantages in chemical interpretation, as
opposed to FB. The Mulliken charges, however, are ill-defined for non-minimal
basis sets, and thus require alternate charge definitions to alleviate this
basis set dependence. This issue arises from the near-redundancy of LCAO
expansion with large basis sets and is exacerbated in crystals. The intrinsic
atomic orbital (IAO) method[24] is one partial charge estimate that has
successfully been applied for molecules.

In this article, we introduce Bloch intrinsic atomic orbitals (Bloch IAOs) as
the natural periodic extension of IAOs, and then present the overall
optimisation scheme to generate localised WFs with Bloch IAOs. The initial
guess is a crucial step in the optimisation and we propose a simple and
effective procedure for generating localised orbitals by defining a natural
gauge and by constructing diabatic Bloch orbitals and diabatic Wannier
functions. We then present and analyse performance and stability of the
optimisation, with particular discussion of the solver’s performance when
separating core and valence bands. Finally, the chemical interpretability of
IBAO WFs is commented upon, using a surface adsorption system as an example.

##  II Theory

###  II.1 Review of Wannier functions

Under Born-von-Karman (BvK) boundary conditions, within a linear-combination-
of-atomic-orbitals (LCAO) framework, the crystal orbitals are expanded in a
basis of Bloch AOs. The Bloch AOs are eigenfunctions of the momentum operator
with crystal momentum wave vector $\mathbf{k}$, and are defined as

$\ket{\mu_{\mathbf{k}}}=\frac{1}{\sqrt{N}}\sum_{\bm{R}}^{N}e^{i\mathbf{k}\cdot\bm{R}}\ket{\tilde{\mu}_{\bm{R}}}.$
(1)

$N$ is the number of unit cells within the BvK ‘supercell’, and $\bm{R}$ is
the lattice vector of the unit cell. $\ket{\tilde{\mu}_{\bm{R}}}$ is the
infinite sum of real-space AOs from each of the supercells throughout the
crystal, and is periodic under the BvK boundary conditions. The crystal
orbitals, also referred to as Bloch functions, are eigenstates of the one-
particle Hamiltonian of a periodic system and are given by

$\ket{\psi_{i,\mathbf{k}}}=\sum_{\mu}C_{\mu,i}^{\mathbf{k}}\ket{\mu_{\mathbf{k}}}.$
(2)

The Bloch functions are delocalised across the entire system. By superimposing
the Bloch functions of a single band across the first Brillouin zone, a
conventional Wannier function[3], centred on a unit cell given by lattice
vector $\bm{R}$, is given by,

$\ket{\phi_{i,\bm{R}}}=\frac{1}{\sqrt{N}}\sum_{\mathbf{k}}^{N}e^{-i\mathbf{k}\cdot\bm{R}}\ket{\psi_{i,\mathbf{k}}}.$
(3)

The WFs span the same space as their Bloch counterparts, with translational
copies found in each unit cell.

Bloch functions are defined up to an arbitrary phase only. However, the
spatial distribution of the resultant WFs are highly dependant on the relative
phases of the contributing Bloch functions. The WFs are thus gauge variant. To
obtain localised conventional WFs, the relative phases of the Bloch functions,
for a given band, must be optimised. By rotating the gauge such that the Bloch
functions appear smooth in reciprocal space, the resulting WFs in real space
are in turn localised, as a property of Fourier transforms,

$\ket{\phi_{i,\bm{R}}}=\frac{1}{\sqrt{N}}\sum_{\mathbf{k}}^{N}e^{-i\mathbf{k}\cdot\bm{R}}e^{i\theta^{\mathbf{k}}_{i}}\ket{\psi_{i,\mathbf{k}}}.$
(4)

A natural gauge for each Bloch function can be defined by requiring that the
scalar product of the coefficients between Bloch functions at $\mathbf{k}$ and
the $\Gamma$-point $\mathbf{0}$ is real valued. By first computing the phase
difference

$\sum_{\mu}C_{\mu,i}^{*\mathbf{k}}C_{\mu,i}^{\mathbf{0}}=Re^{i\theta^{\mathbf{k}}_{i}},$
(5)

the Bloch functions can be rotated into their natural gauge
$\ket{\psi_{i,\mathbf{k}}^{\mathrm{n}}}$, from their original gauge,
$\ket{\psi_{i,\mathbf{k}}^{\mathrm{o}}}$, straightforwardly,

$\ket{\psi_{i,\mathbf{k}}^{\mathrm{n}}}=\ket{\psi_{i,\mathbf{k}}^{\mathrm{o}}}e^{-i\theta^{\mathbf{k}}_{i}}.$
(6)

For bands that have small dispersion, or minimal mixing, imposing the natural
gauge is often sufficient to produce well-localised Wannier functions.

###  II.2 Generalised localised Wannier functions

Generalised Wannier functions are defined by allowing Bloch functions from
several bands to mix[11],

$\ket{\phi_{i,\bm{R}}}=\frac{1}{\sqrt{N}}\sum_{\mathbf{k}}^{N}\sum_{j}e^{-i\mathbf{k}\cdot\bm{R}}U_{j,i}^{\mathbf{k}}\;\ket{\psi_{j,\mathbf{k}}}.$
(7)

By allowing mixing between bands, WFs can be further localised not only to a
unit cell, but also to atomic sites within a unit cell. The resulting
generalised WFs are even more strongly ‘non-unique’ than the conventional
definition, and their locality is highly dependant on the choice of
$U_{ji}^{\mathbf{k}}$. To ensure real-valued WFs, inversion symmetry about the
$\Gamma$-point must be imposed between the Bloch functions,
$\ket{\psi_{i,\mathbf{k}}}=(\ket{\psi_{i,-\mathbf{k}}})^{*}$. The Bloch
functions at the $\Gamma$-point are real, following convention, and the choice
of unitary carries the constraint of
$\bm{U}^{\mathbf{k}}=(\bm{U}^{\mathbf{-k}})^{*}$.

Localisation of WFs is achieved by varying $\bm{U}^{\mathbf{k}}$ to optimise a
chosen locality metric. The Foster-Boys[8, 9] and Pipek-Mezey[10] metrics are
two of the most important examples. The FB method localises orbitals by
defining a metric which minimises the orbital spread, as given by its
variance,

${}\braket{O}_{\textbf{FB}}=\sum_{i}\bra{\phi_{i,\mathbf{0}}}\mathbf{r^{2}}\ket{\phi_{i,\mathbf{0}}}-\bra{\phi_{i,\mathbf{0}}}\mathbf{r}\ket{\phi_{i,\mathbf{0}}}^{2}.$
(8)

Marzari and Vanderbilt[11] generalised the FB approach, originally conceived
for molecules, to evaluate localised WFs, creating so called Foster-Boys
Wannier functions (FBWFs). The Wannier90 package[12], which employs this
method, has been widely used amongst the solid-state community[13].

The original Pipek-Mezey metric was defined as the sum of squares of the
Mulliken partial charges[10]. Jónsson et al.[14] introduced a scheme to
generate localised WFs using the PM metric (PMWFs),

${}\braket{O}_{\textbf{PM}}=\sum_{\bm{R},A,i}|Q_{i}^{A_{\bm{R}}}|^{p}=\sum_{\bm{R},A,i}\bra{\phi_{i,\mathbf{0}}}\hat{P}_{A_{\bm{R}}}\ket{\phi_{i,\mathbf{0}}}^{p},$
(9)

where $Q_{i}^{A_{\bm{R}}}$ is the Mulliken charge[10] associated with WF $i$
on atom $A$, situated in unit cell $\bm{R}$, evaluated using the WFs located
in the reference unit cell. $\hat{P}_{A_{\bm{R}}}$ projects onto a basis of
AOs centred on atom $A_{\bm{R}}$, given by,

$\hat{P}_{A_{\bm{R}}}=\sum_{\mu\in
A}\sum_{\nu,\bm{R^{\prime}}}\ket{\tilde{\mu}_{\bm{R}}}(S^{-1})_{\mu,\bm{R}}^{\nu,\bm{R^{\prime}}}\bra{\tilde{\nu}_{\bm{R^{\prime}}}},$
(10)

where
$S_{\mu,\bm{R}}^{\nu,\bm{R^{\prime}}}=\braket{\tilde{\mu}_{\bm{R}}}{\tilde{\nu}_{\bm{R^{\prime}}}}$
is the overlap. Under a LCAO approach, where atomic orbital coefficients are
naturally available, computing the PM metric is extremely straightforward, as
opposed to the FB metric, which has been used more commonly with plane wave
basis sets. In addition, the strong interpretability of PMWFs, producing
orbitals with $\sigma$ and $\pi$ separation[14], further motivates our choice
to utilise this metric, as opposed to the ‘banana’ bonds found within the FBWF
scheme.

The value of the penalty exponent, $p$, is typically $2$, or $4$. If $p$ is
chosen to equal $1$, then Eq.9 reduces to the normalisation criteria of the
bands,

${}\sum_{\bm{R},A}\sum_{i}^{\mathrm{n_{occ}}}Q_{i}^{A_{\bm{R}}}=\mathrm{n_{occ}},$
(11)

where $\mathrm{n_{occ}}$ are the number of occupied bands.

A key issue with the original definition of the PM metric is that the Mulliken
partial charges do not possess a complete basis set limit, meaning they are
undefined for non-minimal basis sets. Alternative charge definitions for MOs
have been suggested[25, 26, 27], that alleviate this basis set dependence.
Lehtola and Jónsson demonstrated that the localised orbitals obtained were
largely independent of the chosen partial charge estimate[25], providing
significant freedom in choice. Other alternative charge definitions have been
implemented for periodic systems, such as real space partitioning of orbital
charge densities[14, 28, 29], or other types of projection onto a pre-
determined set of minimal basis functions[30].

The intrinsic atomic orbital method (IAO), as proposed by Knizia [24], is one
choice of an alternative partial charge estimate that has been employed
successfully for molecules. Using a free-atom minimal basis as a template,
contraction coefficients from the original basis to IAOs are defined such that
the occupied orbitals are exactly represented, which provides a consistent
assignment of charge to atomic centres. Localised MOs using IAOs align well
with chemical intuition, and quantitative measures such as partial charges and
populations are shown to be resistant to changes in original basis, and are
consistent with chemical understanding, leading to the method being
implemented in many quantum chemistry packages[23, 22, 31, 32]. We thus
propose to adapt IAOs to construct a charge metric suitable to localise
Wannier functions, $Q_{i}^{A_{\bm{R}}^{\mathrm{IAO}}}$, using Bloch intrinsic
atomic orbitals (Bloch IAOs).

###  II.3 Bloch intrinsic atomic orbitals

Given the success of IAOs within molecular schemes, we believe a $\mathbf{k}$
space extension to periodic systems would be desirable. Schäfer et al.[33]
demonstrate the use of IAOs to evaluate localised Wannier functions for a
$\Gamma$-point-only calculation, following the molecular formulation as
described by Janowski[34]. Cui et al.[35] construct crystal IAOs, from which
projected AOs are evaluated. We employ similar principles in our
generalisation to $\mathbf{k}$ space, but crucially outline the additional
augmentations needed to construct localised WFs, optimised using the PM
metric, as a full periodic adaption of the IAO method.

We choose to adapt Knizia’s method[24], such that a set of intrinsic atomic
orbitals are constructed for each $\mathbf{k}$ point within our Monkhorst-Pack
quadrature mesh[36]. The Bloch IAOs are able to exactly describe the original
occupied Bloch functions, providing a basis independent charge metric for WFs.

The original Bloch functions (Eq.2) are expressed in terms of Bloch AOs in the
original basis set, labelled $B_{1}$. Analogous to Knizia’s approach, a
minimal basis, $B_{2}$, of free-atom AOs is first chosen, from which
corresponding Bloch AOs are obtained, $\ket{\rho_{\mathbf{k}}}$ where $\rho\in
B_{2}$, from Eq.1.

The following projection operators are defined,

$\hat{P}_{12}^{\mathbf{k}}=\sum_{\begin{subarray}{c}\mu,\nu\in\mathrm{B_{1}}\end{subarray}}\ket{\mu_{\mathbf{k}}}(S^{-1}_{1})_{\mu,\nu}^{\mathbf{k}}\bra{\nu_{\mathbf{k}}}$
(12)

$\hat{P}_{21}^{\mathbf{k}}=\sum_{\begin{subarray}{c}\rho,\sigma\in\mathrm{B_{2}}\end{subarray}}\ket{\rho_{\mathbf{k}}}(S^{-1}_{2})_{\rho,\sigma}^{\mathbf{k}}\bra{\sigma_{\mathbf{k}}},$
(13)

where
$(S_{1})_{\mu,\nu}^{\mathbf{k}}=\langle\mu_{\mathbf{k}}|\nu_{\mathbf{k}}\rangle$
and
$(S_{2})_{\rho,\sigma}^{\mathbf{k}}=\langle\rho_{\mathbf{k}}|\sigma_{\mathbf{k}}\rangle$
are the Bloch-AO overlap matrices in the original and minimal basis sets,
respectively. Using these operators, depolarisied occupied Bloch functions are
obtained through

$\\{\ket{\tilde{\psi}_{i,\mathbf{k}}}\\}=\mathrm{orth}\\{\hat{P}_{12}^{\mathbf{k}}\hat{P}_{21}^{\mathbf{k}}\ket{\psi_{i,\mathbf{k}}}\\}$
(14)

or in matrix form:

$\mathbf{\tilde{C}^{k}}=\mathrm{orth}\\{\mathbf{P_{12}^{k}P_{21}^{k}C^{k}}\\},$
(15)

Here $\mathrm{orth}\\{\\}$ denotes symmetric orthogonalisation and the
transfer matrices are
$\mathbf{P_{12}^{k}}=(\mathbf{S}_{1}^{-1})^{\mathbf{k}}\mathbf{S}_{12}^{\mathbf{k}}$
and
$\mathbf{P_{21}^{k}}=(\mathbf{S}_{2}^{-1})^{\mathbf{k}}\mathbf{S}_{21}^{\mathbf{k}}$,
where
$(S_{12})^{\mathbf{k}}_{\mu\rho}=\langle\mu_{\mathbf{k}}|\rho_{\mathbf{k}}\rangle$
and $\mathbf{S}_{21}^{\mathbf{k}}=\mathbf{S}_{12}^{\mathbf{k}\dagger}$. The
projector onto the depolarised occupied Bloch functions is
$\tilde{O}^{\mathbf{k}}=\sum_{i}\ket{\tilde{\psi}_{i,\mathbf{k}}}\bra{\tilde{\psi}_{i,\mathbf{k}}}$.

The Bloch IAOs are the minimal Bloch AO basis that contains both the
depolarised and polarisation contributions and are defined through

$\ket{\rho_{\mathbf{k}}^{\mathrm{IAO}}}=(O^{\mathbf{k}}\tilde{O}^{\mathbf{k}}+(1-O^{\mathbf{k}})(1-\tilde{O}^{\mathbf{k}}))\hat{P}^{\mathbf{k}}_{12}\ket{\rho_{\mathbf{k}}}.$
(16)

In matrix notation, Eq.16 is given as:

$\begin{split}\mathbf{A^{k}}=&\mathbf{C^{\mathbf{k}}C^{k\dagger}S_{1}^{k}\tilde{C}^{k}\tilde{C}^{k\dagger}S_{12}^{k}}+\\\
&\mathbf{(1-C^{\mathbf{k}}C^{k\dagger}S_{1}^{k})(P_{12}^{k}-\tilde{C}^{k}\tilde{C}^{k\dagger}S_{12}^{k})},\end{split}$
(17)

where $\mathbf{A^{k}}$ are the contraction coefficients from Bloch AOs to
Bloch IAOs at each $\mathbf{k}$ point, and $\mathbf{1}$ is the identity in the
space of $B_{1}$. Janowski[34] and Knizia[24] both outline a simpler
definition for the IAOs, which is equivalent under the assumption that $B_{2}$
can be directly expressed in $B_{1}$. Having implemented both schemes, we note
that the output Bloch IAOs are very similar, with no significant difference in
localisation performance. Finally, the coefficients of the occupied Bloch
functions in the Bloch IAO basis are given by

$\mathbf{C^{k(\mathrm{IAO})}=(S^{k(\mathrm{IAO})})^{-1}A^{k\dagger}S_{1}^{k}C^{k}},$
(18) $\mathbf{S^{k(\mathrm{IAO})}=A^{k\dagger}S_{1}^{k}A^{k}}.$ (19)

In the original molecular implementation, the output IAO coefficients are
symmetrically orthogonalised. However, in the periodic case, we choose not to
do so. The orthogonalisation procedure introduces arbitrary phases to the
Bloch functions in the IAO basis, specifically when obtaining the eigenvectors
of the IAO coefficient matrix. The relative phase differences between
$\mathbf{k}$ points are thus altered compared to the original Bloch functions
expressed in $B_{1}$, leading to issues when optimising the set of unitary
matrices, across the Brillouin zone, in the IAO basis, since they do not
correspond to the original Bloch functions. The simplest solution to remove
this additional gauge problem is to leave the IAOs un-orthogonalised. The
‘depolarised’ Bloch functions, given by Eqs.14 and 15, which are
orthogonalised, avoid this issue, because any phase augmentation is cancelled
in the projector $\tilde{O}^{k}$.

In summary, obtaining the Bloch IAOs is numerically straightforward, requiring
only a free atom basis, and its corresponding Bloch AO overlaps to perform the
matrix multiplication steps. Computation of inverse overlap matrices can be
avoided by solving instead with a Cholesky decomposition. The Bloch
coefficients in the IAO basis and the IAO overlap matrix can then be used in a
Pipek-Mezey style optimisation to obtain optimally local Wannier functions.
The PM projector (Eq.10), in the Bloch IAO basis, is now defined as

$\hat{P}_{A_{\bm{R}}}^{\mathrm{IAO}}=\sum_{\rho\in
A}\sum_{\sigma,\bm{R^{\prime}}}\ket{\tilde{\rho}^{\mathrm{IAO}}_{\bm{R}}}(S^{-1})_{\rho,\bm{R}}^{\sigma,\bm{R^{\prime}}(\mathrm{IAO})}\bra{\tilde{\sigma}^{\mathrm{IAO}}_{\bm{R^{\prime}}}},$
(20)

where
$S_{\rho,\bm{R}}^{\sigma,\bm{R^{\prime}}(\mathrm{IAO})}=\braket{\tilde{\rho}^{\mathrm{IAO}}_{\bm{R}}}{\tilde{\sigma}_{\bm{R^{\prime}}}^{\mathrm{IAO}}}$
is the IAO overlap in real space, obtained from Fourier transforming
$\mathbf{S^{k(\mathrm{IAO})}}$.

###  II.4 Localisation procedure

Recent work obtaining PM localised WFs and MOs have involved global
optimisation algorithms to determine the stationary points of the
functional[30, 14, 37, 25]. Clement et al.[30] recently demonstrated that a
solver using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm leads to
significantly faster convergence compared to previous steepest ascent (SA) or
conjugate gradient implementations. Our localisation procedure uses a BFGS
based algorithm, that we employ in conjunction with our Bloch IAO charges. To
generate an effective initial guess for the optimisation, a novel procedure
generates approximately localised WFs, which we call diabatic Wannier
functions.

####  II.4.1 Diabatic Wannierisation

The initial guess for the WFs is an important step in the localisation
procedure in order to avoid encountering local maxima. Methods which project
Bloch functions onto a set of trial functions have been outlined[11, 13],
whilst other implementations ensure the unitary space is probed fully by
running multiple calculations using randomly sampled unitary matrices[14, 30].
Clement et al.[30] combine random unitary sampling with a procedure to remove
the gauge freedom of the Bloch functions.

As mentioned prior, Bloch functions are defined with an arbitrary gauge
(Eq.4). By fixing the gauge such that the variations between Bloch functions
in $\mathbf{k}$-space are gradual, the Fourier transform produces WFs which
are largely localised to a single cell, serving as an excellent starting guess
for further optimisation. We have defined the natural gauge to be where the
scalar product of the coefficients between Bloch functions within a band at
$\mathbf{k}$ and the $\Gamma$-point $\mathbf{0}$ is real. For generalised WFs
(Eq.7), where the gauge uncertainty is increased by mixing bands, we extend
the intuition of the natural gauge, to construct diabatic Bloch orbitals and
diabatic Wannier functions. First, the Bloch orbitals of the $\Gamma$-point
are localised by orthogonal transformation,

$\displaystyle|\psi_{i,\mathbf{0}}\rangle=\sum_{j}|\psi_{j,\mathbf{0}}\rangle
O_{ji}.$ (21)

The Bloch orbitals of the remaining $\mathbf{k}$-points are then chosen to be
those with maximal similarity with the $\Gamma$-point. The locality of the
orbitals of the $\Gamma$-point is thus transferred diabatically across the
first Brillouin zone. This is obtained by calculating the unitary matrices,
outside the $\Gamma$-point, which give the minimal least squares difference to
the Bloch coefficients of the $\Gamma$-point,

$\underset{\bm{U}_{\mathbf{k}}}{\mathrm{min}}||C_{\mu,p}^{0}-\sum_{j}C_{\mu,j}^{\mathbf{k}}U^{\mathbf{k}}_{j,p}||^{2}.$
(22)

where $||\dots||$ is the Frobenius[38] norm. As this is an example of an
Orthogonal Procrustes problem[39], a solution can be easily obtained via the
singular value decomposition of the product of Bloch coefficients
$\bm{C^{\mathbf{k}\dagger}C^{0}}$,

$\displaystyle|\psi_{i,\mathbf{k}}\rangle$
$\displaystyle=\sum_{\mu,j,j^{\prime}}C_{\mu,j}^{\mathbf{k}}U_{j,j^{\prime}}V_{j^{\prime},i}|\mu_{\mathbf{k}}\rangle$
(23) $\displaystyle\bm{C^{\mathbf{k}\dagger}C^{0}}$
$\displaystyle=\mathbf{U}\mathbf{\Sigma}\mathbf{V}.$ (24)

In this work, we employ a convenient approximate localisation procedure for
the Bloch orbitals of the $\Gamma$-point, where we simply replace them with
the Cholesky vectors of the $\Gamma$-point density, ensuring computation of
the diabatic WFs is fast.

####  II.4.2 Optimisation of the PM metric

We implement a gradient based optimisation method to obtain the stationary
points of the PM functional. Similar to Clement et al.[30] and Lehtola and
Jónsson[37], a Riemannian geometry approach is adopted to maintain the unitary
constraint, as outlined in ref. [40, 41]. This method has proved successful
since the unitary constraint is maintained implicitly, whilst other methods,
including Lagrange multipliers[42], may suffer from slow convergence or only
obtain a solution which only approximately maintains orthonormality.

Given the extensive discussion of the unitary optimisation algorithm in ref.
[41, 40], we only briefly outline our procedure here. Crucially, the PM charge
metric and all associated expressions are evaluated in the Bloch IAO basis,
using $\mathbf{C^{k(\mathrm{IAO})}}$ (Eq.18). The real space IAO overlaps,
$S_{\sigma,\bm{R^{\prime}}}^{\rho,\bm{R}(\mathrm{IAO})}$, are also used.

In the following expressions, the IAO labels are omitted for clarity. The
Bloch IAO charges are defined as

$\begin{split}&Q_{i}^{A_{\bm{R}}^{\mathrm{IAO}}}=\bra{\phi_{i,\mathbf{0}}}\hat{P}_{A_{\bm{R}}}^{\mathrm{IAO}}\ket{\phi_{i,\mathbf{0}}}=\\\
&\frac{1}{N^{2}}\sum_{\rho\in
A}\Bigr{[}\sum_{\begin{subarray}{c}j,\mathbf{k}\end{subarray}}\bar{C}_{\rho,j}^{*\mathbf{k},\bm{R}}\;U_{j,i}^{*\mathbf{k}}\;\Bigr{]}\Bigr{[}\sum_{j^{\prime},\mathbf{k^{\prime}}}C_{\rho,j^{\prime}}^{\mathbf{k^{\prime}},\bm{R}}\;U_{j^{\prime},i}^{\mathbf{k^{\prime}}}\Bigr{]},\end{split}$
(25)

where we have introduced
$C_{\rho,j}^{\mathbf{k},\bm{R}}=C_{\rho,j}^{\mathbf{k}}e^{i\mathbf{k}\bm{R}}$
and
$\bar{C}_{\rho,j}^{\mathbf{k},\bm{R}}=\sum_{\sigma,\bm{R^{\prime}}}C_{\sigma,j}^{\mathbf{k},\bm{R^{\prime}}}S_{\sigma,\bm{R^{\prime}}}^{\rho,\bm{R}}$.
The Euclidean derivative of the PM functional, $\braket{O}_{\textbf{PM}}$,
with respect to the unitary at $\mathbf{k}$, is given by:

$\begin{split}&\frac{\partial\braket{O}_{\textbf{PM}}}{\partial
U_{j,i}^{*\mathbf{k}}}=\frac{\partial\braket{O}_{\textbf{PM}}}{\partial
Q_{i}^{A_{\bm{R}}^{\mathrm{IAO}}}}\frac{\partial
Q_{i}^{A_{\bm{R}}^{\mathrm{IAO}}}}{\partial U_{j,i}^{*\mathbf{k}}}\\\
&=\frac{p}{N^{2}}\sum_{\bm{R},A}|Q_{i}^{A_{\bm{R}}^{\mathrm{IAO}}}|^{p-1}\sum_{\rho\in
A}\bar{C}_{\rho,j}^{*\mathbf{k},\bm{R}}\Bigr{[}\sum_{j^{\prime},\mathbf{k^{\prime}}}C_{\rho,j^{\prime}}^{\mathbf{k^{\prime}},\bm{R}}U_{j^{\prime},i}^{\mathbf{k^{\prime}}}\Bigr{]}\\\
&+C_{\rho,j}^{*\mathbf{k},\bm{R}}\Bigr{[}\sum_{\begin{subarray}{c}j^{\prime\prime},\mathbf{k^{\prime\prime}}\end{subarray}}\bar{C}_{\rho,j^{\prime\prime}}^{\mathbf{k^{\prime\prime}},\bm{R}}\;U_{j^{\prime\prime},i}^{\mathbf{k^{\prime\prime}}}\Bigr{]}.\end{split}$
(26)

The Riemannian gradient, $\bm{G}_{\mathbf{k}}$, can then be transformed from
the Euclidean gradient,
$\bm{\Gamma}_{\mathbf{k}}=\frac{\partial\braket{O}_{\textbf{PM}}}{\partial\bm{U}^{*}_{\mathbf{k}}}|_{\bm{U}_{\mathbf{k}}}$,
by

$\bm{G}_{\mathbf{k}}=\bm{\Gamma}_{\mathbf{k}}(\bm{U}_{\mathbf{k}})^{{\dagger}}-\bm{U}_{\mathbf{k}}(\bm{\Gamma}_{\mathbf{k}})^{{\dagger}}.$
(27)

The ‘two loop recursion’ version of the limited-memory BFGS algorithm
(l-BFGS)[43] is used, as first implemented for WFs by Clement et al.[30], to
obtain a search direction, $\\{\bm{H}_{\mathbf{k}}\\}$. The matrix elements in
the upper triangle of the anti-Hermitian matrices for the Riemannian gradient,
$\\{\bm{G}_{\mathbf{k}}\\}$, form the gradient vector for the l-BFGS
algorithm, ensuring the output search direction is located on the unitary
manifold. Given the requirement for the output WFs to be real, as mentioned
earlier, our gradient vector is comprised of a total
$\frac{N-1}{2}o^{2}+\frac{o(o-1)}{2}$ real numbers, where $o$ is the number of
Bloch functions being localised.

To obtain a suitable step size, $\mu_{\mathrm{opt}}$, an Armijo[41, 44] and
Wolfe line search[43] were both implemented. If the line search along the
l-BFGS direction fails, the search direction is reset to the steepest ascent
vector, and the line search is then repeated.

The unitary matrices at each $\mathbf{k}$ point are updated,

$\bm{U}_{\mathbf{k},\mathrm{new}}=\;e^{\mu_{\mathrm{opt}}\bm{H}_{\mathbf{k}}}\bm{U}_{\mathbf{k},\mathrm{old}},$
(28)

until the norm of the Riemannian gradient decreases below a threshold. The
output unitary can then be applied directly to the original Bloch functions,
to obtain localised WFs in the $B_{1}$ basis, using Eq.7.

##  III Results

###  III.1 Computational Details

The Bloch IAO procedure and the Pipek-Mezey Wannier function localisation have
been implemented in a developmental version of the TURBOMOLE[23] package. The
initial mean-field Bloch functions were obtained through the periodic Hartree-
Fock procedure, within the `riper` module of the TURBOMOLE package[45, 46, 47,
48, 49]. To generate the IAOs, a minimal basis was constructed from HF
calculations of isolated atoms in the cc-pVTZ[50, 50, 51] basis, as already
implemented within TURBOMOLE to construct molecular IAOs. The PM functional
was evaluated with a penalty exponent of $p=4$, rather than $2$, as shown in
Eq.9, due to better localisation for $\pi$ character orbitals, as discussed in
prior works[24, 25].

An Armijo step size method and a Wolfe line search were tested in the
localisation procedure. It is known that a line search fulfilling the Wolfe
conditions ensures stability of the BFGS updates, by ensuring the approximate
Hessian, within our maximisation problem, is negative definite[43]. However,
computation of the line search is costly, since multiple gradient evaluations
are required along the trial direction. By contrast, the Armijo line search
only requires PM metric values to be evaluated, and we observed its
convergence performance to be similar to the Wolfe line search, with shorter
wall times. The Armijo search was employed in all calculations subsequently
discussed in this article.

Table 1:  Insulating and semiconducting systems used to test the Bloch IAO and PM localisation procedure  System | Basis $B_{1}$ | Monkhorst-Pack mesh size |  |  |  |
---|---|---|---|---|---|---
diamond | pob-TZVP[52] | 11,11,11 |  |  |  |
silicon | pob-TZVP | 11,11,11 |  |  |  |
boron nitride | pob-TZVP | 15,15 |  |  |  |
graphene | pob-TZVP | 15,15 |  |  |  |
MgO | def-SVP | 11,11,11 |  |  |  |
SiO2 | pob-TZVP | 5,5,5 |  |  |  |
_trans_ -$\mathrm{(C_{2}H_{2})}_{\infty}$ | pob-TZVP | 101 |  |  |  |
(4,4) C-nanotube | def2-SVP | 11 |  |  |  |

Table 1 details the insulating and semi-conducting systems used to probe the
performance of the IAO PM localisation scheme. The original basis set,
$B_{1}$, used in the mean field calculation, and the Monkhorst-Pack[36] mesh,
are shown. All test systems utilised the universal Coulomb-fitting auxiliary
basis sets[53] with the exception of the magnesium oxide and carbon nanotube
systems, which employed auxiliary functions optimised for the def-SVP and
def2-SVP basis sets, respectively[53]. Unit cell parameters and geometries are
provided in the supplementary information.

###  III.2 Overall performance

Figure 1: Bloch IAO localised Wannier functions of (4,4) nanotube, showing $\sigma$ (left) and $\pi$ (right) bonding character. An isosurface value of 0.05 was used. Table 2:  PM metric values of WFs from the initial SCF calculation, after rotation into the natural gauge, after diabatic Wannierisation and after PM optimisation. Number of iterations to converge PM metric, after initial diabatic preparation, using l-BFGS or steepest ascent. System |  $\braket{O}_{\textbf{PM}}$(SCF) |  $\braket{O}_{\textbf{PM}}$(Nat.) |  $\braket{O}_{\textbf{PM}}$(Dia.) |  $\braket{O}_{\textbf{PM}}$ (Opt.) | l-BFGS | SA
---|---|---|---|---|---|---
diamond | $1.49\times 10^{-8}$ | $6.13\times 10^{-2}$ | 1.92 | 2.66 | 35 | 100
silicon | $2.40\times 10^{-2}$ | $6.64\times 10^{-2}$ | 9.79 | 10.66 | 42 | 114
boron nitride | $1.20\times 10^{-5}$ | 2.36 | 3.19 | 3.37 | 55 | 731
graphene | $3.37\times 10^{-6}$ | $8.37\times 10^{-2}$ | 1.94 | 2.56 | 42 | 60
MgO | $2.88\times 10^{-1}$ | 3.87 | 9.56 | 9.61 | 8 | 744
SiO2 | $5.96\times 10^{-4}$ | 9.55 | 30.80 | 37.29 | 338 | 1794
_trans_ -$\mathrm{(C_{2}H_{2})}_{\infty}$ | $4.92\times 10^{-5}$ | $4.25\times 10^{-2}$ | 2.17 | 2.76 | 43 | 1517
(4,4) C-nanotube | $1.65\times 10^{-3}$ | $4.34\times 10^{-3}$ | 33.70 | 40.60 | 62 | 261

Table 2 reports the performance of the Bloch IAO localisation scheme. The
values of the PM metric are presented for the WFs of the SCF calculation,
after rotation into the natural gauge, after diabatic Wannierisation and after
PM optimisation. The number of iterations required to localise the PM
objective function with l-BFGS, compared to steepest ascent (SA), using
diabatic WFs as the initial guess, are also given. The convergence threshold
for the PM gradient norm was set to $10^{-5}$, with the exception of boron
nitride (BN), which was set to $10^{-6}$. All occupied orbitals are included
in the optimisation.

In all cases, the final PM values from the l-BFGS optimisation and steepest
ascent were equal, to a threshold of $10^{-5}$ (or $10^{-6}$ for BN),
confirming that the output WFs were equally localised. A tighter threshold was
required for BN in order for the output WFs from l-BFGS and SA to agree to
target precision. As demonstrated by Clement et al.[30], we confirm that
utilising l-BFGS, compared to SA, markedly improves convergence performance.
In some systems, a tenfold reduction in iterations to converge is observed, if
not greater. The PM optimiser performs robustly across the range of insulating
and semiconducting materials explored, successfully localising every test
system. With the exception of silicon dioxide, all test systems converge
within 100 iterations with l-BFGS. A larger number of iterations was required
for silicon dioxide. We observed that the step size was very small, which
indicates that the Hessian description of the landscape in this case may be
poor. Despite this, l-BFGS still converges five times faster compared to SA,
for this example, showing the robustness of the scheme.

We stress that using diabatic WFs as the initial preparation has an important
role in the robustness and quality of the final localised WFs. As seen in
Table 2, the PM values after diabatic Wannierisation are remarkably close to
the final optimised values, showing that a significant degree of locality has
been captured through the diabatisation. Our experiments using randomly
generated unitary matrices as the initial guess led to final WFs with PM
values that consistently were smaller than that obtained from the diabatic
preparation, and never greater. The choice of objective functional, and the
parameterisation employed for the gradient, gives an optimisation landscape
with many local maxima and use of an appropriate initial guess, such as
diabatic WFs, is required to ensure the global maxima is located. Direct
comparison of the number of iterations required to converge, using a random
guess and the diabatic preparation, is often not possible since the final WFs
are usually inequivalent. Although we have not verified it in this work, we
predict that localising the Bloch functions of the $\Gamma$-point with an IAO
procedure instead of via Choleski decomposition would further increase
$\braket{O}_{\textbf{PM}}$(Dia.) and would reduce the number of iterations
required for full optimisation.

The initial PM values from straight Wannierisation of the SCF Bloch functions
are very small. It should be noted that due to the different gauges of the
Bloch functions in separate calculations, the values for
$\braket{O}_{\textbf{PM}}$(SCF) can vary arbitrarily. The values reported in
Table 2 are in fact the largest value of $\braket{O}_{\textbf{PM}}$(SCF) taken
from 10 separate calculations. Throughout our testing, we have never observed
any example where $\braket{O}_{\textbf{PM}}$(SCF) has been greater or similar
in magnitude to $\braket{O}_{\textbf{PM}}$(Dia.). Applying the natural gauge
to the Bloch functions increases the PM values in all cases, confirming that
the natural gauge smooths the Bloch functions in reciprocal space. For
crystals with well-separated bands formed from weakly interacting AOs, simply
applying the natural gauge results in well-localised WFs. In all cases,
$\braket{O}_{\textbf{PM}}$(Nat.)) sit in between the SCF and diabatic values.

Figure 1 presents example localised WFs of the (4,4) nanotube system,
described in Table 1. All figures were plotted using Avogadro[54]. The left
subfigure clearly shows a carbon-carbon $\sigma$ bond, whilst the right
subfigure illustrates a $\pi$ character WF. This $\pi$ and $\sigma$ separation
is observed in the other test systems, demonstrating that the Bloch IAO
localisation procedure retains the advantages of the original PM metric.

###  III.3 Separate optimisation of valence and core bands

It is often desirable to localise WFs for the core and valence bands
separately. For example, within local correlation theories for molecules,
obtaining localised occupied orbitals with no contribution from uncorrelated
core orbitals is necessary to compute accurate correlation energies[55]. In
view of this, we report that the Bloch IAO localisation scheme performs
exceptionally well in this scenario, due to the quality of the diabatic WFs,
which almost immediately generates localised core orbitals.

Table 3:

Number of iterations to converge PM metric with initial diabatic preparation,
with core and valence band separation  System Initial core PM gradient core
valence diamond $1.08\times 10^{-4}$ 2 20 silicon $8.73\times 10^{-4}$ 3 24
boron nitride $9.70\times 10^{-4}$ 3 27 graphene $1.91\times 10^{-4}$ 2 22 MgO
$5.17\times 10^{-5}$ 2 4 SiO2 $6.55\times 10^{-3}$ 2 113 _trans_
-$\mathrm{(C_{2}H_{2})}_{\infty}$ $1.67\times 10^{-4}$ 3 22 (4,4) C-nanotube
$6.81\times 10^{-2}$ 3 29

Table 3 shows the norm of the PM gradient for the core bands after the
diabatic initial preparation, and the subsequent iterations required to
localise with the PM metric, using l-BFGS, for both the core and valence
bands, across all the test systems. The PM gradient norms of our initial core
WFs are all already close to the convergence threshold ($10^{-5}$, or
$10^{-6}$ for BN), and localise within 3 l-BFGS iterations, demonstrating that
the diabatic preparation yields nearly optimally local core WFs. The number of
iterations required to localise the core increases typically only by 1 or 2
iterations when the convergence threshold is increased to $10^{-8}$. Across
all test systems, the difference between the final localised PM values for the
core orbitals and $\mathrm{n_{core}}$ was within $3\times 10^{-2}$, where
$\mathrm{n_{core}}$ is given by Eq.11, summing only over core bands. This
indicates that the core orbitals are localised to the global maxima, given
that $\mathrm{n_{core}}$ is the upper bound for the PM metric for the core.
Localisation of the valence bands also occurs rapidly, although less markedly
than in the core case, and the total number of iterations required to localise
the core and valence bands separately is less than that required to localise
the full occupied space (Table 2), across all the test systems. This is to be
expected, given that the dimensionality of the optimisation problem is reduced
by separating the core and valence bands.

Diabatic preparation is particularly effective in localising WFs from bands
composed of weakly interacting AOs. The core bands and valence bands with
strong ionic character, for example in MgO, require only a few optimisation
steps for optimisation. Since the localised valence WFs are usually bonding in
character, they are typically centred between atoms, and inherently less local
than their core counterparts and require more steps for optimisaiton. Even for
the core bands, using diabatic WFs as an initial preparation is vital for the
robustness and accuracy of the Bloch IAO scheme. Using a random initial
unitary for core bands frequently leads to output WFs with PM values
significantly smaller than those obtained with the diabatic guess, showing the
encountering local maxima is a common problem without correct preparation of
the WFs.

Figure 2: Bloch IAO localised Wannier functions of the MgO(001) CO adsorption
system. A magnesium 2p-like orbital is shown (left), as well a $\pi$ bonding
orbital on CO (centre). A WF demonstrating oxygen 3p orbital donation, on the
MgO surface, to the CO $\pi^{*}$ orbital, is also presented (right). Only
atoms in the reference cell are shown. An isosurface value of 0.01 was used.

###  III.4 Chemical intuition of Bloch IAO generated WFs

One of the key strengths of the original IAO scheme is the clear
interpretation of these molecular orbitals, and the direct connection to
chemical intuition and concepts. Knizia[24] demonstrated that IAOs allow
robust, basis set independent, partial charges and orbital populations to be
computed. This enables quantitative measures for electronegativities and
oxidation states for molecules, which align with empirical understanding, to
be evaluated. In analogous fashion, we demonstrate that Bloch IAO localised
WFs provide chemical understanding in periodic systems. Bloch IAO partial
charges can be computed in similar fashion to moleuclar IAOs, by summing the
atomic contributions across all the cells within the supercell,

${}q_{A}=Z_{A}-\sum_{\bm{R},i}Q_{i}^{A_{\bm{R}}^{\mathrm{IAO}}},$ (29)

where $Z_{A}$ is the atom’s nuclear charge. IAO partial charges can be
computed for any periodic system, such as those listed in Table 1. However, it
is also worth stressing that IAOs can be used for both periodic and molecular
systems, on equal footing, meaning that they are a robust and consistent
partial charge estimate for probing both material and molecular systems. This
opens the possibility to investigate interesting chemical scenarios, such as
systems containing interactions between material and molecules.

The adsorption of CO onto the MgO(001) surface has been heralded as the
‘hydrogen molecule of surface science’[56], and an important case study for
the theoretical understanding of heterogenous catalysis. Obtaining an accurate
adsorption energy for this system remains a highly discussed topic, in which
many quantum chemical and many body methods have been utilised[57, 58, 59, 60,
61], in order to achieve consensus with experimental data. This adsorption
example is an ideal case to demonstrate the ability that Bloch IAO localised
WFs have to provide insight into the underlying chemistry of the system.

To model the system, a unit cell consisting of a $4\times 4\times 1$ slab of
MgO was constructed. CO, orientated perpendicular to the surface, was
positioned with a C-Mg equilibrium distance of $2.479$Å[58]. To obtain mean-
field Bloch functions, a periodic DFT (PBE[62]) calculation was conducted in
the `riper` module, using the pob-TZVP[52] basis set, on a (3,3) Monkhorst-
Pack mesh to sample the Brillouin zone. Bloch IAO localised WFs were then
obtained.

Although higher level quantum chemical methods have been used elsewhere to
attempt to accurately model the weak van der Waals interactions dominating the
adsorption, we stress that our motivation is to use this system to exhibit the
use of IAO WFs for chemical intuition. Since orbitals provide a zeroth order
description to the motion of the electrons, and are a result from mean-field,
effective one electron theories, using solely DFT to model this picture serves
our purposes.

Figure 2 presents three Bloch IAO localised WFs of the system, all three of
which align with chemical understanding. The left and centre subfigures show a
localised 2p orbital, centred on Mg, and a $\pi$ bonding orbital, centred on
CO, respectively. Similar WFs are observed when WFs for the surface and
adsorbate are computed separately. A more interesting WF is presented on the
right, where back-bonding from the nearest neighbour oxygen atom, on the MgO
surface, to the $\pi^{*}$ orbital of carbon monoxide, is shown. Although the
role of back-bonding within metal oxide adsorption has been scrutinised[61,
63] , we wish to make the point that IAO localised WFs provide direct and
intuitive chemical understanding, consistent with the level of theory used to
generate the original Bloch functions.

Table 4:

Bloch IAO partial charges of the non-interacting and equilibrium MgO(001) CO
adsorption system  System C (CO) O (CO) Mg (MgO) O (MgO) non-interacting 0.42
-0.42 1.68 -1.68 equilibrium 0.32 -0.41 1.70111The partial charge of the
closest Mg atom to CO is given -1.65222The average partial charge of the four
nearest neighbours to CO is given

The Bloch IAO partial charges of the MgO(001) CO adsorption system, in the
equilibrium geometry, were also calculated. In Table 4, we compare these
partial charge values to charges obtained from separate periodic calculations
of the MgO surface slab, and the CO molecule, representing a non-interacting
scenario. Most significantly, a reduction in positive partial charge on the
carbon of CO is observed moving from the non-interacting to equilibrium
geometry, as well as a decrease in negative charge on the oxygen nearest
neighbours of MgO. This can be rationalised from the back-bonding process
presented in Figure 2, and once again shows that quantitative measures, such
as partial charges, derived from IAO WFs, are consistent with chemical
intuition, at the level of theory employed in the mean field picture.

##  IV Conclusions

We have generalised the intrinsic atomic orbital method to periodic solids.
Bloch IAOs form a minimal basis which exactly represent the occupied bands,
and thus alleviate the well known issue of using Mulliken charges for non-
minimal basis sets. They thus enable localised Wannier functions, optimised
using the Pipek-Mezey metric, to be robustly evaluated, as first introduced by
Jónsson et al.[14] We outline a localisation scheme, which prepares the
initial Bloch functions by diabatically transferring locality imposed at the
gamma point through the Brillouin zone, before localising according to the PM
metric. Clement et al.[30] demonstrated the improved performance using l-BFGS
compared to other gradient based solvers, and we confirm this. This scheme
works efficiently across a range of semiconducting and insulating solids, and
in particular we highlight the ability of the diabatic WFs to localise atom-
centred WFs almost immediately. Using the example of CO adsorption onto
MgO(001), we demonstrate that Bloch IAO localised WFs can provide chemical
insight into systems, through visualisation of the WFs, and through computing
measures such as partial charges. We expect that Bloch IAOs will provide a
bridge for understanding chemical phenomena within periodic systems. Bloch
IAOs are not solely restricted to LCAO methods, but can also be applied with
plane wave basis sets by computing the overlap between plane waves and the
minimal Gaussian AO basis[33]. In particular, we note that in the molecular
setting, localised MOs using IAOs have proved popular in constructing
localised occupied orbitals and domains, for use within local correlation
theories[55]. Bloch IAOs may provide an analogous route for similar
implementations within periodic systems.

###### Acknowledgements.

Financial support from the University of Oxford and Turbomole GmbH is
gratefully acknowledged.

## References

  * Hohenberg and Kohn [1964] P. Hohenberg and W. Kohn, Inhomogeneous electron gas, Phys. Rev. 136, B864 (1964).
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, Self-consistent equations including exchange and correlation effects, Phys. Rev. 140, A1133 (1965).
  * Wannier [1937] G. H. Wannier, The structure of electronic excitation levels in insulating crystals, Phys. Rev. 52, 191 (1937).
  * Saebo and Pulay [1993] S. Saebo and P. Pulay, Local treatment of electron correlation, Annual Review of Physical Chemistry 44, 213 (1993).
  * Hampel and Werner [1996] C. Hampel and H. Werner, Local treatment of electron correlation in coupled cluster theory, The Journal of Chemical Physics 104, 6286 (1996).
  * Maschio  _et al._ [2007] L. Maschio, D. Usvyat, F. R. Manby, S. Casassa, C. Pisani, and M. Schütz, Fast local-mp2 method with density-fitting for crystals. i. theory and algorithms, Phys. Rev. B 76, 075101 (2007).
  * Usvyat  _et al._ [2011] D. Usvyat, B. Civalleri, L. Maschio, R. Dovesi, C. Pisani, and M. Schütz, Approaching the theoretical limit in periodic local MP2 calculations with atomic-orbital basis sets: The case of LiH, The Journal of Chemical Physics 134, 214105 (2011).
  * Boys [1960] S. F. Boys, Construction of some molecular orbitals to be approximately invariant for changes from one molecule to another, Rev. Mod. Phys. 32, 296 (1960).
  * Foster and Boys [1960] J. M. Foster and S. F. Boys, Canonical configurational interaction procedure, Rev. Mod. Phys. 32, 300 (1960).
  * Pipek and Mezey [1989] J. Pipek and P. G. Mezey, A fast intrinsic localization procedure applicable for ab initio and semiempirical linear combination of atomic orbital wave functions, The Journal of Chemical Physics 90, 4916 (1989).
  * Marzari and Vanderbilt [1997] N. Marzari and D. Vanderbilt, Maximally localized generalized wannier functions for composite energy bands, Phys. Rev. B 56, 12847 (1997).
  * 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, Journal of Physics: Condensed Matter 32, 165902 (2020).
  * 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).
  * Jónsson  _et al._ [2017] E. O. Jónsson, S. Lehtola, M. Puska, and H. Jónsson, Theory and applications of generalized pipek–mezey wannier functions, Journal of Chemical Theory and Computation 13, 460 (2017), pMID: 28099002.
  * Kresse and Hafner [1993] G. Kresse and J. Hafner, Ab initio molecular dynamics for liquid metals, Phys. Rev. B 47, 558 (1993).
  * Giannozzi  _et al._ [2017] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. A. D. Jr, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, Advanced capabilities for materials modelling with quantum espresso, Journal of Physics: Condensed Matter 29, 465901 (2017).
  * Clark  _et al._ [2005] S. J. Clark, M. D. Segall, C. J. Pickard, P. J. Hasnip, M. I. J. Probert, K. Refson, and M. C. Payne, First principles methods using castep, Zeitschrift für Kristallographie - Crystalline Materials 220, 567 (2005).
  * Gonze  _et al._ [2016] X. Gonze, F. Jollet, F. Abreu Araujo, D. Adams, B. Amadon, T. Applencourt, C. Audouze, J.-M. Beuken, J. Bieder, A. Bokhanchuk, E. Bousquet, F. Bruneval, D. Caliste, M. Côté, F. Dahm, F. Da Pieve, M. Delaveau, M. Di Gennaro, B. Dorado, C. Espejo, G. Geneste, L. Genovese, A. Gerossier, M. Giantomassi, Y. Gillet, D. Hamann, L. He, G. Jomard, J. Laflamme Janssen, S. Le Roux, A. Levitt, A. Lherbier, F. Liu, I. Lukačević, A. Martin, C. Martins, M. Oliveira, S. Poncé, Y. Pouillon, T. Rangel, G.-M. Rignanese, A. Romero, B. Rousseau, O. Rubel, A. Shukri, M. Stankovski, M. Torrent, M. Van Setten, B. Van Troeye, M. Verstraete, D. Waroquiers, J. Wiktor, B. Xu, A. Zhou, and J. Zwanziger, Recent developments in the abinit software package, Computer Physics Communications 205, 106 (2016).
  * Erba  _et al._ [2023] A. Erba, J. K. Desmarais, S. Casassa, B. Civalleri, L. Donà, I. J. Bush, B. Searle, L. Maschio, L. Edith-Daga, A. Cossard, C. Ribaldone, E. Ascrizzi, N. L. Marana, J.-P. Flament, and B. Kirtman, Crystal23: A program for computational solid state physics and chemistry, Journal of Chemical Theory and Computation 19, 6891 (2023), pMID: 36502394.
  * [20] C. Pisani and R. Dovesi, Exact-exchange hartree–fock calculations for periodic systems. i. illustration of the method, International Journal of Quantum Chemistry 17, 501.
  * [21] R. Dovesi, C. Pisani, and C. Roetti, Exact-exchange hartree–fock calculations for periodic systems. ii. results for graphite and hexagonal boron nitride, International Journal of Quantum Chemistry 17, 517.
  * Sun  _et al._ [2020] Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt, N. A. Bogdanov, G. H. Booth, J. Chen, Z.-H. Cui, _et al._ , Recent developments in the pyscf program package, The Journal of chemical physics 153 (2020).
  * Balasubramani  _et al._ [2020] S. G. Balasubramani, G. P. Chen, S. Coriani, M. Diedenhofen, M. S. Frank, Y. J. Franzke, F. Furche, R. Grotjahn, M. E. Harding, C. Hättig, A. Hellweg, B. Helmich-Paris, C. Holzer, U. Huniar, M. Kaupp, A. Marefat Khah, S. Karbalaei Khani, T. Müller, F. Mack, B. D. Nguyen, S. M. Parker, E. Perlt, D. Rappoport, K. Reiter, S. Roy, M. Rückert, G. Schmitz, M. Sierka, E. Tapavicza, D. P. Tew, C. van Wüllen, V. K. Voora, F. Weigend, A. Wodyński, and J. M. Yu, Turbomole: Modular program suite for ab initio quantum-chemical and condensed-matter simulations, The Journal of Chemical Physics 152, 184107 (2020).
  * Knizia [2013] G. Knizia, Intrinsic atomic orbitals: An unbiased bridge between quantum theory and chemical concepts, Journal of chemical theory and computation 9, 4834 (2013).
  * Lehtola and Jónsson [2014] S. Lehtola and H. Jónsson, Pipek–mezey orbital localization using various partial charge estimates, Journal of Chemical Theory and Computation 10, 642 (2014), pMID: 26580041.
  * Cioslowski [1991] J. Cioslowski, Partitioning of the orbital overlap matrix and the localization criteria, Journal of Mathematical Chemistry 8, 169 (1991).
  * [27] D. R. Alcoba, L. Lain, A. Torre, and R. C. Bochicchio, An orbital localization criterion based on the theory of “fuzzy” atoms, Journal of Computational Chemistry 27, 596.
  * Yu and Trinkle [2011] M. Yu and D. R. Trinkle, Accurate and efficient algorithm for Bader charge integration, The Journal of Chemical Physics 134, 064111 (2011).
  * Zicovich-Wilson  _et al._ [2016] C. Zicovich-Wilson, M. Hô, A. Navarrete-López, and S. Casassa, Hirshfeld-i charges in linear combination of atomic orbitals periodic calculations, Theoretical Chemistry Accounts 135 (2016).
  * Clement  _et al._ [2021] M. C. Clement, X. Wang, and E. F. Valeev, Robust pipek–mezey orbital localization in periodic solids, Journal of Chemical Theory and Computation 17, 7406 (2021).
  * Werner  _et al._ [2020] H.-J. Werner, P. J. Knowles, F. R. Manby, J. A. Black, K. Doll, A. Heßelmann, D. Kats, A. Köhn, T. Korona, D. A. Kreplin, _et al._ , The molpro quantum chemistry package, The Journal of chemical physics 152 (2020).
  * Saue  _et al._ [2020] T. Saue, R. Bast, A. S. P. Gomes, H. J. A. Jensen, L. Visscher, I. A. Aucar, R. Di Remigio, K. G. Dyall, E. Eliav, E. Fasshauer, _et al._ , The dirac code for relativistic molecular calculations, The Journal of chemical physics 152 (2020).
  * Schäfer, Tobias and Gallo, Alejandro and Irmler, Andreas and Hummel, Felix and Grüneis, Andreas [2021] Schäfer, Tobias and Gallo, Alejandro and Irmler, Andreas and Hummel, Felix and Grüneis, Andreas, Surface science using coupled cluster theory via local Wannier functions and in-RPA-embedding: The case of water on graphitic carbon nitride, The Journal of Chemical Physics 155, 244103 (2021).
  * Janowski [2014] T. Janowski, Near equivalence of intrinsic atomic orbitals and quasiatomic orbitals, Journal of Chemical Theory and Computation 10, 3085 (2014).
  * Cui  _et al._ [2020] Z.-H. Cui, T. Zhu, and G. K.-L. Chan, Efficient implementation of ab initio quantum embedding in periodic systems: Density matrix embedding theory, Journal of Chemical Theory and Computation 16, 119 (2020).
  * Monkhorst and Pack [1976] H. J. Monkhorst and J. D. Pack, Special points for brillouin-zone integrations, Phys. Rev. B 13, 5188 (1976).
  * Lehtola and Jónsson [2013] S. Lehtola and H. Jónsson, Unitary optimization of localized molecular orbitals, Journal of Chemical Theory and Computation 9, 5365 (2013), pMID: 26592274.
  * Horn and Johnson [2012] R. A. Horn and C. R. Johnson,  _Matrix analysis_ (Cambridge university press, 2012).
  * Gower and Dijksterhuis [2004] J. C. Gower and G. B. Dijksterhuis,  _Procrustes problems_ , Vol. 30 (OUP Oxford, 2004).
  * Abrudan  _et al._ [2009] T. Abrudan, J. Eriksson, and V. Koivunen, Conjugate gradient algorithm for optimization under unitary matrix constraint, Signal Processing 89, 1704 (2009).
  * Abrudan  _et al._ [2008] T. E. Abrudan, J. Eriksson, and V. Koivunen, Steepest descent algorithms for optimization under unitary matrix constraint, IEEE Transactions on Signal Processing 56, 1134 (2008).
  * Hyvärinen  _et al._ [2001] A. Hyvärinen, J. Karhunen, and E. Oja,  _Independent Component Analysis_, Adaptive and Cognitive Dynamic Systems: Signal Processing, Learning, Communications and Control (Wiley, 2001).
  * Nocedal and Wright [2006] J. Nocedal and S. J. Wright,  _Numerical Optimization_ , 2nd ed. (Springer, New York, NY, USA, 2006).
  * Polak [1997] E. Polak,  _Optimization: Algorithms and Consistent Approximations_ (Springer-Verlag, Berlin, Heidelberg, 1997).
  * Łazarski  _et al._ [2015] R. Łazarski, A. M. Burow, and M. Sierka, Density functional theory for molecular and periodic systems using density fitting and continuous fast multipole methods, Journal of Chemical Theory and Computation 11, 3029 (2015).
  * [46] R. Łazarski, A. M. Burow, L. Grajciar, and M. Sierka, Density functional theory for molecular and periodic systems using density fitting and continuous fast multipole method: Analytical gradients, Journal of Computational Chemistry 37, 2518.
  * Burow and Sierka [2011] A. M. Burow and M. Sierka, Linear scaling hierarchical integration scheme for the exchange-correlation term in molecular and periodic systems, Journal of Chemical Theory and Computation 7, 3097 (2011).
  * Burow  _et al._ [2009] A. M. Burow, M. Sierka, and F. Mohamed, Resolution of identity approximation for the Coulomb term in molecular and periodic systems, The Journal of Chemical Physics 131, 214101 (2009).
  * [49] C. Müller, M. Sharma, and M. Sierka, Real-time time-dependent density functional theory using density fitting and the continuous fast multipole method, Journal of Computational Chemistry 41, 2573.
  * Dunning [1989] T. Dunning, Gaussian basis sets for use in correlated molecular calculations. i. the atoms boron through neon and hydrogen, The Journal of Chemical Physics 90, 1007 (1989).
  * Prascher  _et al._ [2011] B. P. Prascher, D. E. Woon, K. A. Peterson, T. H. Dunning, and A. K. Wilson, Gaussian basis sets for use in correlated molecular calculations. vii. valence, core-valence, and scalar relativistic basis sets for li, be, na, and mg, Theoretical Chemistry Accounts 128, 69 (2011).
  * [52] M. F. Peintinger, D. V. Oliveira, and T. Bredow, Consistent gaussian basis sets of triple-zeta valence with polarization quality for solid-state calculations, Journal of Computational Chemistry 34, 451.
  * Weigend [2006] F. Weigend, Accurate coulomb-fitting basis sets for h to rn, Phys. Chem. Chem. Phys. 8, 1057 (2006).
  * Hanwell  _et al._ [2012] M. D. Hanwell, D. E. Curtis, D. C. Lonie, T. Vandermeersch, E. Zurek, and G. R. Hutchison, Avogadro: an advanced semantic chemical editor, visualization, and analysis platform, Journal of Cheminformatics 4, 17 (2012).
  * [55] Q. Ma and H.-J. Werner, Explicitly correlated local coupled-cluster methods using pair natural orbitals, WIREs Computational Molecular Science 8, e1371.
  * Sauer [2019] J. Sauer, Ab initio calculations for molecule–surface interactions with chemical accuracy, Accounts of Chemical Research 52, 3502 (2019), pMID: 31765121.
  * Nygren and Pettersson [1996] M. A. Nygren and L. G. M. Pettersson, Comparing ab initio computed energetics with thermal experiments in surface science: CO/MgO(001), The Journal of Chemical Physics 105, 9339 (1996).
  * Mitra  _et al._ [2022] A. Mitra, M. R. Hermes, M. Cho, V. Agarawal, and L. Gagliardi, Periodic density matrix embedding for co adsorption on the mgo(001) surface, The Journal of Physical Chemistry Letters 13, 7483 (2022).
  * Shi  _et al._ [2023] B. X. Shi, A. Zen, V. Kapil, P. R. Nagy, A. Grüneis, and A. Michaelides, Many-body methods for surface chemistry come of age: Achieving consensus with experiments, Journal of the American Chemical Society 145, 25372 (2023), pMID: 37948071.
  * Ye and Berkelbach [2023] H.-Z. Ye and T. C. Berkelbach, Co adsorption on the surface of mgo from periodic coupled-cluster theory with local natural orbitals: Adding to the consensus (2023), arXiv:2309.14651 [cond-mat.mtrl-sci] .
  * Neyman and Rösch [1992] K. M. Neyman and N. Rösch, Co bonding and vibrational modes on a perfect mgo(001) surface: Lcgto-ldf model cluster investigation, Chemical Physics 168, 267 (1992).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Pacchionia  _et al._ [1994] G. Pacchionia, K. M. Neyman, and N. Rosch, Co adsorption on the (001) surface of mgo: a comparison of hartree-fock and local density functional results, Journal of Electron Spectroscopy and Related Phenomena 69, 13 (1994).