% \section{REVIEW OF BASIC THEORY}
% \label{sec:theory}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Bloch functions and Wannier functions}
\label{sec:bloch-wannier}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Electronic structure calculations are often carried out using periodic
boundary conditions. This is the most natural choice for the study of
perfect crystals, and also applies to the common use of periodic
supercells for the study of non-periodic systems such as liquids,
interfaces, and defects. The one-particle effective Hamiltonian
$H$ 
%
then commutes with the lattice-translation operator
${T}_{\bf R}$, allowing one to choose as common eigenstates the
Bloch orbitals $\vert\,\psi_{n{\bf k}}\,\rangle$:
%
%\nmm{I removed, as discussed, all hats from chap 2 - a few are in 5
%  and 6, and plenty in 7 and 8 - need to chekc with Arash}
%
\begin{equation}
[\,{H},\,{T}_{\bf R}\,]\;=\;0 
\;\;\Rightarrow\;\; \psi_{n{\bf k}}({\bf r})\;
=\;u_{n{\bf k}}({\bf r})\,e^{i{\bf k\cdot r}}\;\;,
\label{eq:bloch}
\end{equation}
%
where $u_{n{\bf k}}({\bf r})$ has the periodicity of the Hamiltonian.

Several Bloch functions are sketched on the left-hand side of
Fig.~\ref{fig:wannier-bloch} for a toy model in which the
band of interest is composed of $p$-like orbitals centered on each atom.
For the moment, we suppose that this band is an isolated band,
i.e., it remains separated by a gap from the bands below and
above at all $\kk$.  Since Bloch functions at different $\kk$ have
different envelope functions $e^{i\kk\cdot\rr}$, one can expect
to be able to build a localized ``wave packet'' by superposing
Bloch functions of different $\kk$.  To get a very localized wave
packet in real space, we need to use a very broad superposition
in $\kk$ space.  But $\kk$ lives in the periodic Brilloin zone,
so the best we can do is to choose equal amplitudes all across
the Brillouin zone.  Thus, we can construct
%
\beq
w_0(\rr) \;=\;\ibz d\kk\; \psi_{n{\bf k}}(\rr)\ \;\;,
\label{eq:wannier-home}
\eeq
%
where $V$ is the real-space primitive cell volume and the integral is
carried over the Brillouin zone (BZ).
(See Sec.~\ref{sec:norm} for normalization conventions.)
\equa{wannier-home}
can be interpreted as the WF located in the
home unit cell, as sketched in the top-right panel of
Fig.~\ref{fig:wannier-bloch}.

\begin{figure}
\begin{center}
\includegraphics[width=7cm]{figs/wannier-bloch}
\end{center}
\caption{\label{fig:wannier-bloch} (Color online)
Left: Bloch functions associated with a single band in 1D, at
three different values of wavevector $\kk$.  Right:
WFs associated with the same band, forming periodic
images of one another.  Blue dots indicate atoms; green curves
indicate envelopes $e^{ikx}$ of the Bloch functions.
Bloch and WFs span the same Hilbert space.}
\end{figure}

More generally, we can insert a phase factor $e^{-i\kk\cdot\RR}$ into
the integrand of \equ{wannier-home}, where $\RR$ is a real-space
lattice vector; this has the effect of translating the real-space
WF by $\RR$, generating additional WFs
such as $w_1$ and $w_2$ sketched in Fig.~\ref{fig:wannier-bloch}.
Switching to the Dirac bra-ket notation and introducing the notation
that $\RR n$ refers to the WF $w_{n\RR}$ in cell $\RR$
associated with band $n$, WFs can be
constructed according to \cite{wannier-pr37}
%
\begin{equation}
\label{eq:wanniertransform}
\ket{\RR n}\;=\; \ibz d{\bf k}\;e^{-i{\bf k\cdot {\bf R}}} \,
\ket{\,\psi\bnk\,} \;\;.
\end{equation}
%
It is easily shown that the $\vert\,{{\bf R}n}\,\rangle$ form
an orthonormal set (see Sec.~\ref{sec:norm})
%
% \ism{\tiny Should we expand the footnote on normalization conventions, since
%   this is a continuous source of confusion? Do I understand correctly
%   that our convention is that the Bloch functions $\psi\bnk(\rr)$ are
%   normalized over a unit cell (so that
%   $\inprod{\psi_{m\kk'}}{\psi_{n\kk}}$ is a unit cell integral), while
%   WFs $\ket{\RR n}$ are instead normalized over a periodic supercell,
%   so that $\inprod{\RR' m}{\RR n}$ is a supercell integral? If so, I
%   think we need to insert by hand factors of $1/N_{\rm cells}$ when
%   expressing the latter in terms of the former via
%   \equ{wanniertransform}.  (The alternative would be to insert
%   explicitly a factor of $\sqrt{N_{\rm cells}}$ in
%   \equ{wanniertransform}, as done, e.g., in Eq.~(9a) of
%   \textcite{king-smith-prb93}.) I believe this would also address the
%   concern raised by Arash in a marginpar around \equ{bloch-w-gauge}.
%   (Arash, do you agree?)\\
% AM: \textcolor{blue}{Yes, this is precisely the issue I refer to later. 
%   Formally, I prefer the latter solution, namely
%   to explicitly include the factor of $\sqrt{N}$ in the definition. In 
%   practice, however, it will be easier to state in a footnote the 
%   normalisation convention that
%   we are using, and I think that this is what we should do.}} 
%
and that two WFs 
$\vert\,{{\bf R}n}\,\rangle$ and $\vert\,{{\bf R}^\prime}n\,\rangle$
transform into each other under translation by the lattice
vector ${\bf R}-{\bf R^\prime}$ \cite{blount-ssp62}.
\equ{wanniertransform} takes the form of a Fourier
transform, and its inverse transform is
%
\begin{equation}
\label{eq:wannier-inverse}
\ket{\,\psi\bnk\,} \;=\;\sum_\RR \,e^{i\kk\cdot\RR}\, \ket{\RR n}
\end{equation}
%
(see Sec.~\ref{sec:norm}).
Any of the Bloch functions on the left side of
Fig.~\ref{fig:wannier-bloch} can thus be built up by linearly superposing
the WFs shown on the right side, when the appropriate phases
$e^{i\kk\cdot\RR}$ are used.

The transformations of \eqs{wanniertransform}{wannier-inverse}
constitute a unitary transformation between Bloch and Wannier
states.  Thus, both sets of states
provide an equally valid description of the band subspace,
even if the WFs are not Hamiltonian eigenstates. 
For example, the charge density obtained by
summing the squares of the Bloch functions $\ket{\psi\bnk}$ or
the WFs $\ket{\RR n}$ is identical;
a similar reasoning applies to the trace of any one-particle operator.
The equivalence between the Bloch and the Wannier representations
can also be made manifest by expressing the band projection
operator ${P}$ in both representations, i.e. as
%
\begin{equation}
{P}=\ibz d{\bf k}\;\ket{\psi\bnk}\bra{\psi\bnk}=\sum_\RR \ket{\RR n}\bra{\RR n}\;.
\label{eq:bandproj}
\end{equation}
%
% Furthermore, a Slater determinant constructed from the
% $\ket{\psi\bnk}$ is guaranteed to be identical to one constructed
% from the $\ket{\RR n}$.

WFs thus provide an attractive option for representing the space
spanned by a Bloch band in a crystal, being localized
while still carrying the same information contained in
the Bloch functions.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Gauge freedom}
\label{sec:gauge-freedom}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

However, the theory of WFs is
made more complex by the presence of a ``gauge freedom'' that exists in
the definition of the $\psi_{n{\bf k}}$.  In fact, we can replace
%
\begin{equation}
\label{eq:gauge-psi}
\vert\,\tilde{\psi}_{n{\bf k}}\,\rangle\;=\;e^{i\varphi_n({\bf k})}\,\vert\,\psi_{n{\bf k}}\,\rangle\;\;,
\end{equation}
or equivalently,
\begin{equation}
\label{eq:gauge-u}
\vert\,\tilde{u}_{n{\bf k}}\,\rangle\;=\;e^{i\varphi_n({\bf k})}\,\vert\,u_{n{\bf k}}\,\rangle\;\;,
\end{equation}
%
without changing the physical description of the system, with
$\varphi_n({\bf k})$ being any real function that is periodic
in reciprocal space.\footnote{
  More precisely, the condition is that
  $\varphi_n({\bf k}+{\bf G})=\varphi_n({\bf k})+{\bf G}\cdot\Delta{\bf R}$
  for any reciprocal-lattice translation ${\bf G}$, where $\Delta{\bf R}$
  is a real-space lattice vector.  This allows for the possibility that
  $\varphi_n$ may shift by $2\pi$ times an integer upon translation by
  $\bf G$; the vector $\Delta{\bf R}$ expresses the corresponding
  shift in the position of the resulting WF.
}
A {\it smooth} gauge could e.g. be defined such that
$\nabla_{\bf k}\vert u_{n{\bf k}}\rangle$ is well defined at all
$\bf k$. Henceforth we shall assume that the Bloch functions on
the right-hand side of \equ{wanniertransform} belong to a smooth
gauge, since we would not get well-localized
WFs on the left-hand side otherwise.
This is typical of Fourier transforms: the smoother the
reciprocal-space object, the more localized the resulting
real-space object, and vice versa (see \ref{sec:wfprop} for more rigorous relations 
between smoothness of a function and asymptotic decay of 
its Fourier transform).

One way to see this explicitly
is to consider the $\RR=0$ home-cell $w_{n0}(\rr)$
evaluated at a distant point $\rr$; using \equ{bloch}
in \equ{wanniertransform}, this is
given by $\int_{\rm BZ}\,\,u\bnk(\rr) e^{i\kk\cdot\rr}\,
d\kk$, which will be small due to cancellations arising from
the rapid variation of the exponential factor, provided that
$u\bnk$ is a smooth function of $\kk$ \cite{blount-ssp62}.

It is important to realize that the gauge freedom of
\eqs{gauge-psi}{gauge-u} propagates into the WFs.  That is, different
choices of smooth gauge correspond to different sets of WFs
having in general different shapes and spreads. In this sense, the WFs are ``more
non-unique'' than the Bloch functions, which only acquire a phase
change.  We also emphasize that there is no
``preferred gauge'' assigned by the Schr\"odinger equation to the
Bloch orbitals.  Thus, the non-uniqueness of the WFs resulting from 
\equ{wanniertransform} is unavoidable.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Multiband case}
\label{sec:multi-band}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Before discussing how this non-uniqueness might be resolved, we first
relax the condition that band $n$ be a single isolated band,
and consider instead a manifold of $\nbw$ bands
%
%\dvm{\scr\textcolor{red}{I implemented Ivo's suggestion about notation
%  for number of bands and WFs, but with a twist: I propose to use:\\
%  $\backslash$nbw when we don't need to distingluish between number
%  of bands and WFs.\\
%  $\backslash$nbands and $\backslash$nwann when we do want to
%  distinguish them, as for disentanglement.\\
%  Currently these are $J$, $\cal J$, and $J$, respectively, but
%  it is easy to redefine them in defs.tex.\\
%  I tried to insert these macros everywhere there were needed in
%  Ch.~2, but I might have missed one or two.  In any case, we
%  need to propagate them into later chapters too.}}
%
that remain separated with respect to any lower or higher bands
outside the manifold.  Internal degeneracies and crossings among the
$\nbw$ bands may occur in general.  In the simplest case this manifold
corresponds to the occupied bands of an insulator, but more generally
it consists of any set of bands that is separated by a gap from both
lower and higher bands everywhere in the Brillouin zone.  Traces over
this band manifold are invariant with respect to any unitary
transformation among the $\nbw$ occupied Bloch orbitals at a given
wavevector, so it is natural to generalize the notion of a ``gauge
transformation'' to
%
\begin{equation}
\vert\,\tilde{\psi}_{n{\bf k}} \,\rangle=
\sum_{m=1}^\nbw\,U_{mn}^{({\bf k})}\,\vert\,\psi_{m{\bf k}} \,\rangle
\;\;.
\label{eq:blochlike}
\end{equation}
%
Here $U_{mn}^{({\bf k})}$ is a unitary matrix of dimension $\nbw$ that
is periodic in $\bf k$, with
\equ{gauge-psi} corresponding to the special case of a diagonal
$U$ matrix.  It follows that the projection operator onto this
band manifold at wavevector $\kk$,
%
\begin{equation}
{P}_{\,\kk}=\sum_{n=1}^\nbw\ket{\psi\bnk}\bra{\psi\bnk}
           =\sum_{n=1}^\nbw\ket{\tilde{\psi}\bnk}\bra{\tilde{\psi}\bnk}
\label{eq:proj-equiv}
\end{equation}
%
is invariant, even though the $\vert\tilde{\psi}_{n{\bf k}}
\rangle$ resulting from \equ{blochlike} are no longer generally
eigenstates of ${H}$, and $n$ is no longer a band index in
the usual sense.
%
% \footnote{This formulation may also be applied to the case of
% a partially occupied manifold, but in this case the
% occupation factor $f_{n\kk}$ appearing in the Hamiltonian gauge
% has to be replaced by an occupation matrix $\tilde{f}_{mn,\kk}$ in
% the transformed gauge.}


Our goal is again to construct WFs out of these
transformed Bloch functions using \equ{wanniertransform}.
Figs.~\ref{fig:si-gaas}(a) and (b) show, for example, what the result
might eventually look like for the case of the four occupied
valence bands of Si or GaAs, respectively.  From these four bands,
one obtains four equivalent WFs per unit cell, each
localized on one of the four nearest-neighbor Si-Si or Ga-As bonds.
The presence of a bond-centered inversion symmetry for Si, but not
GaAs, is clearly reflected in the shapes of the WFs.

\begin{figure}
\begin{center}
\includegraphics[angle=0,width=5cm]{figs/rmp_j1}
\end{center}
\caption{\label{fig:si-gaas} (Color online) MLWFs constructed from the
four valence bands of Si (left) and GaAs (right; Ga at
upper right, As at lower left), having the
character of $\sigma$-bonded combinations of $sp^3$ hybrids. The red and
blue isosurfaces correspond to two opposite values for the amplitudes of the
real-valued MLWFs.}
\end{figure}

Once again, we emphasize that the gauge freedom expressed in
\equ{blochlike} implies that the WFs are strongly
non-unique.  This is illustrated in Fig.~\ref{fig:gaas}, which
shows an alternative construction of WFs for GaAs.
The WF on the left was constructed from the lowest
valence band $n$=1, while the one on the right is one of three
constructed from bands $n$=2-4.  The former has primarily
As $s$ character and the latter has primarily As $p$ character,
although both (and especially the latter) contain some Ga $s$ and $p$
character as well.  The WFs of Figs.~\ref{fig:si-gaas}(b)
and Fig.~\ref{fig:gaas} are related to each other by a
certain manifold of 4$\times$4 unitary matrices $U_{nm}^{(\kk)}$
relating their Bloch transforms in the manner of \equ{blochlike}.

\begin{figure}[b]
\begin{center}
\includegraphics[angle=0,width=9cm]{figs/rmp_j2}
\end{center}
\caption{\label{fig:gaas} (Color online) MLWFs constructed from
the $s$ band (left) or from the three $p$ bands (right) of GaAs.}
\end{figure}

However, before we can arrive at well-localized WFs like
those shown in Figs.~\ref{fig:si-gaas} and \ref{fig:gaas}, we again
have to address questions of smoothness of the gauge choice expressed
in \equ{blochlike}.  This issue is even more profound in the
present multiband case, since this smoothness criterion
is generally {\it incompatible} with the usual construction of
Bloch functions.  That is, if we simply insert the usual Bloch
functions $\ket{\psi\bnk}$, defined to be eigenstates of ${H}$,
into the right-hand side of \equ{wanniertransform}, it is
generally {\it not} possible to produce well-localized WFs.
The problem arises when there are degeneracies among
the bands in question at certain locations in the Brillouin zone.
Consider, for example, what happens if we try to
construct a single WF from the highest band $n$=4
in GaAs.  This would be doomed to failure, since this band becomes
degenerate with bands $n$=2-3 at the zone center $\Gamma$ as shown
in Fig.~\ref{fig:gaas}, and as a result, band 4 is non-analytic in
$\kk$ in the vicinity of $\Gamma$.  The Fourier transform of
\equ{wanniertransform} would then result in a poorly localized
object having power-law tails in real space.

In such cases, therefore, the extra unitary mixing expressed in
\equ{blochlike} is mandatory, even if it may be optional in the
case of a set of discrete bands that do not touch anywhere in the BZ.
So, generally speaking, our procedure must be that we
%
start from a set of Hamiltonian eigenstates
$\ket{\psi\bnk}$ that are not per se smooth in
$\kk$, and introduce unitary rotations $U_{mn}^{(\kk)}$ that
``cancel out'' the discontinuities in such a way that smoothness
is restored, i.e., that the resulting
$\ket{\tilde{\psi}\bnk}$ of \equ{blochlike} obey the smoothness
condition that $\nabla_\kk\ket{\tilde{\psi}\bnk}$ remains
regular at all $\kk$.  Then, when these $\ket{\tilde{\psi}\bnk}$ are
inserted into \equ{wanniertransform} in place of the $\ket{\psi\bnk}$,
well-localized WFs should result.  Explicitly,
this results in WFs constructed according to
%
\begin{equation}
\vert\,{{\bf R}n}\,\rangle\;= \;\ibz d{\bf k}
\,e^{-i{\bf k\cdot {\bf R}}}
\,\sum_{m=1}^\nbw\,U_{mn}^{({\bf k})}\,\vert\,\psi_{m{\bf k}} \,\rangle
\;\;.
\label{eq:wannier}
\end{equation}

The question remains how to choose the unitary rotations
$U_{mn}^{(\kk)}$ so as to accomplish this task.
We will see that one way to do this is to use a projection
technique, as outlined in the next section.  Ideally, however,
we would like the construction to result in
WFs that are ``maximally localized'' according
to some criterion.  Methods for accomplishing this are discussed
in Sec.~\ref{sec:maxloc}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Normalization conventions}
\label{sec:norm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In the above equations, we have adopted the convention that
Bloch functions are normalized to one unit cell,
$\int_V d\rr\,|\psi_{n\kk}(\rr)|^2=1$, even though they extend
over the entire crystal.  With this notation,
$\inprod{\psi_{n\kk}}{\psi_{n\kk}}$ is not unity; instead,
it diverges according to the rule
%
\beq
\inprod{\psi_{n\kk}}{\psi_{m\kk'}}=\frac{(2\pi)^3}{V}\,\delta_{nm}
\,\delta^3(\kk-\kk') \;\;.
\eeq
%
With these conventions it is easy to check, for example, that
Eqs.~(\ref{eq:wanniertransform}) and (\ref{eq:wannier-inverse})
are consistent.

It is often more convenient to work on a discrete
uniform $\kk$ mesh instead of continuous $\kk$
space.\footnote{The discretization of $\kk$-space amounts to imposing
  periodic boundary conditions on the Bloch wavefunctions over a
  supercell in real space, and it should be kept in mind that
  technically speaking the WFs given by
  \eqs{discrete-abnormal}{discrete} are not localized, as they also
  display the supercell periodicity. Under these circumstances the
  notion of ``Wannier localization'' refers to localization {\it
    within} one supercell. This is sharply defined provided that a
  sufficiently large supercell has been chosen to insure negligible
  overlap between a WF and its periodic images. \label{foot:discretization}}  Letting $N$ be the
number of unit cells in the 
periodic supercell, or equivalently, the number of
mesh points in the BZ, it is possible to keep the conventions close to
the continuous case by defining the Fourier transform pair as
%
\bea
\label{eq:discrete-abnormal}
\ket{\psi_{n\kk}}&=&\sum_\RR\,e^{i\kk\cdot\RR}\,\ket{\RR n} \; \nonumber\\
&\Updownarrow& \\
\ket{\RR n}&=&\frac{1}{N}\sum_\kk \,e^{-i\kk\cdot\RR}\, \ket{\psi_{n\kk}}
\nonumber
\eea
%
with $\inprod{\psi_{n\kk}}{\psi_{m\kk'}}=N\delta_{nm}\,\delta_{\kk\kk'}$.
More commonly, however, one writes
%
\bea
\label{eq:discrete}
\ket{\psi_{n\kk}}&=&
   \frac{1}{\sqrt{N}}\sum_\RR\,e^{i\kk\cdot\RR}\,\ket{\RR n} \; \nonumber\\
&\Updownarrow& \\
\ket{\RR n}&=&
   \frac{1}{\sqrt{N}}\sum_\kk \,e^{-i\kk\cdot\RR}\, \ket{\psi_{n\kk}} \nonumber
\eea
%
with $\inprod{\psi_{n\kk}}{\psi_{m\kk'}}=\delta_{nm}\,\delta_{\kk\kk'}$, and this is
the convention we'll use in this work (from its first instance \equ{r} onwards).
\equa{bandproj}, which expresses the equivalence of the projectors
onto Bloch and Wannier manifolds, is generalized to the multiband case
in the present notation as
%
\begin{equation}
{P}=\sum\bnk \ket{\psi\bnk}\bra{\psi\bnk}
       =\sum_{n\RR} \ket{\RR n}\bra{\RR n}\;.
\label{eq:mbandproj}
\end{equation}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Wannier functions via projection}
\label{sec:projection}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
A simple yet often effective approach for constructing a smooth gauge
in $\kk$, and a corresponding set of well-localized WFs, is by 
projection - an approach that finds its roots in the analysis of 
\textcite{cloizeaux-pr64a}. Here, as discussed e.g. in Sec. IV.G.1 of
\textcite{marzari-prb97}, one starts from a set of $\nbw$ trial
orbitals $g_n(\bf r)$ corresponding to some rough guess for the WFs in
the home unit cell.  These $g_n(\bf r)$ are projected onto the Bloch
manifold at wavevector $\bf k$ to obtain
%
\beq
\vert\phi_{n\bf k}\rangle=\sum_{m=1}^\nbw\,\vert\psi_{m\bf k} \rangle
\langle\psi_{m\bf k}\vert g_n\rangle\;,
\label{eq:proj}
\eeq
%
which are typically smooth in $\kk$-space, albeit not orthonormal.  We
note that in actual practice such projection is achieved by first
computing a matrix of inner products $(A_{\bf
  k})_{mn}=\langle\psi_{m\bf k}\vert g_n\rangle$ and then using these
in \equ{proj}.  The overlap matrix $(S_{\bf k})_{mn}=\langle\phi_{m\bf
  k}\vert\phi_{n\bf k}\rangle =(A_\kk^\dagger
A_\kk^{\phantom{\dagger}})_{mn}$ is then computed and used to
construct the L\"owdin-orthonormalized Bloch-like states
%
\begin{equation}
\vert\tilde{\psi}_{n\bf k}\rangle = \sum_{m=1}^\nbw\, 
\vert\phi_{m\bf k}\rangle(S_{\bf k}^{-1/2})_{mn} \;.
\label{eq:symm}
\end{equation}
%
These $\vert\tilde{\psi}_{n\bf k}\rangle$ have now a smooth gauge in $\kk$, 
are related to the original
$\vert\psi_{n\bf k} \rangle$ by a unitary transformation,\footnote{One
  can prove that this transformation is unitary by performing the
  singular value decomposition $A=ZDW^\dagger$,
  with $Z$ and $W$ unitary and $D$ real and diagonal. It follows that
  $A(A^\dagger A)^{-1/2}$ is equal to $ZW^\dagger$, and thus
  unitary.}  and when substituted into \equ{wanniertransform} in place
of the $\vert\psi_{n\kk}\rangle$ result in a set of well-localized WFs.
We note in passing that the $\vert\tilde{\psi}_{n\bf k}\rangle$ 
 are uniquely defined by the trial orbitals $g_n(\bf r)$ and the
chosen (isolated) manifold, since any arbitrary unitary rotation among the
$\vert\psi_{n\bf k}\rangle$ orbitals cancels out exactly and does not
affect the outcome of Eq. \ref{eq:proj}, eliminating thus any gauge
freedom.

We emphasize that the trial functions do not need to resemble
the target WFs closely; it is often sufficient
to choose simple analytic functions (e.g., spherical harmonics times
Gaussians) provided
they are roughly located on sites where WFs are expected to be centered and
have appropriate angular character.  The method is successful
as long as the inner-product matrix $A_\kk$ does not
become singular (or nearly so) for any $\kk$, which can be insured by
checking that the ratio of maximum and minimum values of $\det S_\kk$
in the Brillouin zone does not become too large.  For example,
spherical ($s$-like) Gaussians located at the bond centers will
suffice for the construction of well-localized WFs, akin to those
shown in Fig.~\ref{fig:si-gaas}, spanning the four occupied valence bands
of semiconductors such as Si and GaAs.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Maximally localized Wannier functions}
\label{sec:maxloc}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The projection method described in the
previous subsection has been used by many authors
\cite{stephan-prb00,ku-prl02,lu-prb04,qian-prb08},
as has a related method involving downfolding of the band structure
onto a minimal basis \cite{andersen-prb00,zurek-c05}; some of these approaches
will also be discussed in Sec.~\ref{sec:proj-down}.  Other authors have made
use of symmetry considerations, analyticity requirements, and
variational procedures
\cite{sporkmann-prb94,sporkmann-jpcm97,smirnov-prb01}.
A very general and now widely used approach, however, has been that 
developed by \textcite{marzari-prb97} in which localization
is enforced by introducing a well-defined {\it localization
criterion}, and then refining the $U_{mn}^{({\bf k})}$ in order
to satisfy that criterion.  We first give an overview of this approach,
and then provide details in the following subsections.

First, the localizaton functional
\begin{eqnarray}
\Omega&=&\sum_n\left[    
\langle\,{{\bf 0}n}\,\vert\,r^2\,\vert\,{{\bf 0}n}\,\rangle\,-\,
\langle\,{{\bf 0}n}\,\vert\,{\bf r}\,\vert\,{{\bf 0}n}\,\rangle^2\,
\right]\nonumber\\
&=&\sum_n\left[    
\langle r^2\rangle_n-\bar{\bf r}_n^2
\right]
\label{eq:omega}
\end{eqnarray}
%
is defined, measuring of the sum of the quadratic spreads of the $\nwann$
WFs in the home unit cell around their centers.  This turns out to be
the solid-state equivalent of the Foster-Boys criterion of quantum
chemistry \cite{boys-rmp60,foster-rmp60a,foster-rmp60b,boys-66}.  The
next step is to express $\Omega$ in terms of the Bloch functions. This
requires some care, since expectation values of the position operator
are not well defined in the Bloch representation.  The needed
formalism will be discussed briefly in Sec.~\ref{sec:rsr} and more
extensively in Sec.~\ref{sec:bp}, with much of the conceptual work
stemming from the earlier development the modern theory of
polarization
\cite{blount-ssp62,resta-f92,king-smith-prb93,vanderbilt-prb93,resta-rmp94}.

Once a $\kk$-space expression for $\Omega$ has been derived, maximally
localized WFs are obtained by minimizing it with respect to the
$U_{mn}^{({\bf k})}$ appearing in \equ{wannier}.  This is done as a
post-processing step after a conventional electronic-structure
calculation has been self-consistently converged and a set of
ground-state Bloch orbitals $\vert\,\psi_{m{\bf k}}\,\rangle$ has been
chosen once and for all. The $U_{mn}^{({\bf k})}$ are then iteratively
refined in a direct minimization procedure of the localization
functional that is outlined in Sec.~\ref{sec:local} below.  This
procedure also provides the expectation values $\langle r^2\rangle_n$
and $\bar{\bf r}_n$; the latter, in particular, are the primary
quantities needed to compute many of the properties, such as the
electronic polarization, discussed in Sec.~\ref{sec:pol}.  If desired,
the resulting $U_{mn}^{({\bf k})}$ can also be used to construct
explicitly the maximally localized WFs via \equ{wannier}.  This step
is typically only necessary, however, if one wants to visualize the
resulting WFs or to use them as basis functions for some subsequent
analysis.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Real-space representation}
\label{sec:rsr}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

An interesting consequence stemming from the choice of
(\ref{eq:omega}) as the localization functional is that it allows
a natural decomposition of the functional into gauge-invariant
and gauge-dependent parts.  That is, we can write
%
\begin{equation}
\label{eq:spread-decomp}
\Omega\,=\,\Omega_{\rm\,I}\,+\,\widetilde{\Omega}
%\,=\, \Omega_{\rm\,I}\,+\,\Omega_{\rm D}\,+\,\Omega_{\rm OD}
\end{equation}
%
where
%
\begin{equation}
\Omega_{\rm\,I} = \sum_n \left[ \langle 
\,{{\bf 0}n}\,\vert\,r^2\,\vert\,{{\bf 0}n}\,
\rangle\,-\,\sum_{{\bf R}m}
\,\bigl\vert \langle{{\bf R}m}\vert{\bf r}\vert
{\bf 0}n\rangle\bigr\vert^2 \right]
\label{eq:OmegaI}
\end{equation}
%
and
%
\begin{equation}
\widetilde{\Omega}=\sum_n \sum_{{\bf R}m\ne{\bf
0}n} \bigl\vert \langle{{\bf R}m}\vert{\bf r}\vert
{\bf 0}n\rangle\bigr\vert^2 \,.
\label{eq:Omegatil}
\end{equation}
%
It can be shown that not only $\widetilde{\Omega}$ but also
$\Omega_{\rm\,I}$ is {\it positive definite}, and moreover that
$\Omega_{\rm\,I}$ is {\it gauge-invariant}, i.e., invariant under any
arbitrary unitary transformation (\ref{eq:wannier}) of the Bloch
orbitals \cite{marzari-prb97}.  This follows straightforwardly from
recasting \equ{OmegaI} in terms of the band-group projection operator
$P$, as defined in \equ{mbandproj}, and its complement $Q=1-P$:
%
\begin{eqnarray}
\Omega_{\rm\,I} &=& \sum_{n\alpha}
\langle{{\bf 0}n}\vert r_\alpha Q r_\alpha \vert{\bf 0}n\rangle
\nonumber \\
&=& \sum_{\alpha} {\rm Tr}_{\rm c}\, [P r_\alpha Q r_\alpha ]
\;.
\label{eq:posdef}
\end{eqnarray}
%
The subscript `c' indicates trace per unit cell.  Clearly
$\Omega_{\rm\,I}$ is gauge invariant, since it is expressed in terms of
projection operators that are unaffected by any gauge transformation.
It can also be seen to be positive definite by using the idempotency
of $P$ and $Q$ to write
$\Omega_{\rm\,I}
= \sum_\alpha {\rm Tr}_{\rm c}\, [(P r_\alpha Q)(P r_\alpha Q)^\dagger ]
= \sum_\alpha ||Pr_\alpha Q||_{\rm c}^2$.

The minimization procedure of $\Omega$ thus actually corresponds to
the minimization of the 
non-invariant
%non gauge-invariant 
part $\widetilde\Omega$
only.  At the minimum, the off-diagonal elements $\left\vert
\langle{{\bf R}m}\vert{\bf r}\vert {\bf 0}n\rangle\right\vert^2 $ are
as small as possible, realizing the best compromise in the
simultaneous diagonalization, within the subpace of
the Bloch bands considered, of the three position operators $x$, $y$
and $z$, which do not in general commute when projected 
onto this space.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Reciprocal-space representation}
\label{sec:recipspace}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

As shown by \textcite{blount-ssp62}, matrix elements of the position
operator between WFs take the form
%
\begin{equation}
\langle{\bf R}n\vert{\bf r}\vert{\bf 0}m\rangle = i\,{\frac{V}{(2\pi)^3}}
\int d{\bf k} \, e^{i{\bf k}\cdot{\bf R}}
\langle u_{n{\bf k}}\vert\nabla_{\bf k}\vert u_{m{\bf k}}\rangle 
\label{eq:rmatel}
\end{equation}
and
\begin{equation}
\langle{\bf R}n\vert r^2 \vert{\bf 0}m\rangle = -
{\frac{V}{(2\pi)^3}} \int d{\bf k}\,
e^{i{\bf k}\cdot{\bf R}}
\langle u_{n{\bf k}}\vert\nabla_{\bf k}^2\vert u_{m{\bf k}}\rangle
\;\;.
\label{eq:rrmatel}
\end{equation}
These expressions provide the needed connection with our underlying Bloch
formalism, since they allow to express the localization functional 
$\Omega$
in terms of the matrix elements of $\nabla_{\bf k}$ and $\nabla_{\bf k}^2$.
In addition, they allow to calculate the effects on the localization of
any unitary transformation of the $\vert u_{n{\bf k}}\rangle$  
without having to recalculate expensive (especially when
plane-wave basis sets are used) scalar products.
We thus determine the Bloch orbitals $\vert u_{m{\bf k}} \rangle$ 
on a regular mesh
of {\bf k}-points, and use finite differences to evaluate the above
derivatives. 
In particular, we make the assumption 
that the BZ has been discretized into a uniform 
Monkhorst-Pack mesh, and the Bloch orbitals have been determined on
that mesh.\footnote{Even the case of $\Gamma$-sampling --
where the Brillouin zone is sampled with a single {\bf k}-point -- is 
encompassed by the above formulation. 
In this case the neighboring {\bf k}-points
are given by reciprocal lattice vectors ${\bf G}$
and the Bloch orbitals differ only by phase factors
$\exp(i{\bf G\cdot r})$ from their counterparts at $\Gamma$. 
The algebra does become simpler, though, as will be discussed
in Sec.~\ref{sec:gpoint}.}

For any $f(\bf k)$ that is a smooth function of $\bf k$, it can be shown
that its gradient can
be expressed by finite differences as
%
\begin{equation}
\nabla f({\bf k})=\sum_{\bf b}\, w_b \, {\bf b}
 \,[f({\bf k+b})-f({\bf k})]
+{\cal O}(b^2)
 \;\;
\label{eq:grad}
\end{equation}
calculated on stars (``shells") of near-neighbor {\bf k}-points;
here $\bf b$ is a vector connecting a $\kk$-point to one of its 
neighbors, $w_b$ is an appropriate geometric factor that
depends on the number of points in the star and its
geometry (see Appendix B in \textcite{marzari-prb97} 
and \textcite{mostofi-cpc08} for a detailed
description). 
%
In a similar way,
%
\begin{equation}
\vert\nabla f({\bf k})\vert^2=\sum_{\bf b}\, w_b
 \,[f({\bf k+b})-f({\bf k})]^2 
+{\cal O}(b^3)
\;\;.
\label{eq:nabla}
\end{equation}

It now becomes straightforward to calculate the matrix elements in
\eqs{rmatel}{rrmatel}.  All the information needed for the the
reciprocal-space derivatives is encoded in the overlaps
between Bloch orbitals at neighboring {\bf k}-points
%
\begin{equation}
M_{mn}^{(\bf k,b)} = \langle u_{m\bf k}\vert u_{n,\bf k+b}\rangle\;\;.
\label{eq:Mkb}
\end{equation}
%
These overlaps play a central role in the formalism, since all 
contributions to the localization functional can be expressed 
in terms of them. Thus, once the $M_{mn}^{(\bf k,b)}$ have been calculated,
no further interaction with the electronic-structure code that calculated
the ground state wavefunctions is necessary - making the entire
Wannierization procedure a code-independent post-processing step
\footnote{In particular, see \textcite{ferretti-jpcm07} for the
extension to ultrasoft pseudopotentials and the projector-augmented wave method,
and \textcite{posternak-prb02,freimuth-prb08,kunes-cpc10} for the full-potential linearized augmented
planewave method.}.
There is no unique form for the localization functional in terms of the overlap
elements, as it is possible to
write down many alternative finite-difference expressions for
$\bar{\bf r}_n$ and $\langle r^2 \rangle_n$ which agree numerically to
leading order in the mesh spacing $b$ (first and second order for
$\bar{\bf r}_n$ and $\langle r^2 \rangle_n$ respectively).  We give
here the expressions of \textcite{marzari-prb97}, which have the
desirable property of transforming correctly under gauge
transformations that shift $\ket{\0 n}$ by a lattice vector.  They are
%
\begin{equation}
\bar{\bf r}_n = - {\frac{1}{N}} \sum_{\bf k,b}\, w_b \, {\bf b} \,
{\rm Im}\,\ln M_{nn}^{(\bf k,b)} 
\label{eq:r}
\end{equation}
(where we use, as outlined in Sec.~\ref{sec:norm}, the convention of \equ{discrete}), and
%
%
\begin{equation}
\langle r^2 \rangle_n = {\frac{1}{N}} \sum_{\bf k,b}\, w_b \,
\left\{
  \left[
     1-\vert M_{nn}^{(\bf k,b)}\vert^2
  \right] +
  \left[
     {\rm Im}\,\ln M_{nn}^{(\bf k,b)}
  \right]^2 
\right\} 
\;\;.
\label{eq:rr}
\end{equation}
%
The corresponding expressions for the gauge-invariant and
gauge-dependent parts of the spread functional are
%
\begin{equation}
\label{eq:omega_i}
\Omega_{\rm I} = {\frac{1}{N}} \sum_{\bf k,b}\, w_b \,
\Big( \nbw-\sum_{mn} \vert M_{mn}^{(\bf k,b)} \vert^2 \Big)
\end{equation}
%
and 
%
\begin{eqnarray}
\label{eq:omtilde}
\widetilde\Omega &=& 
\frac{1}{N} \sum_{\bf k,b} w_b
\sum_{m\ne n} \vert M_{mn}^{\bf (k,b)} \vert^2 \\ \nonumber
&+& {\frac{1}{N}} \sum_{\bf k,b} w_b
\sum_n \left( - {\rm Im}\,\ln M_{nn}^{\bf (k,b)} 
 - {\bf b}\cdot\bar{\bf r}_n \right)^2 \;\;.
\end{eqnarray}
%
As mentioned, it is possible to write down alternative discretized
expressions which agree numerically with
Eqs.~(\ref{eq:r})--(\ref{eq:omtilde}) up to the orders indicated in
the mesh spacing $b$; at the same time, one needs to be careful in
realizing that certain quantities, such as the spreads, will display
slow convergence with respect to the BZ sampling (see \ref{sec:gpoint}
for a discussion), or that some exact results (e.g. that the sum of
the centers of the Wannier functions is invariant with respect to
unitary transformations) might acquire some numerical noise.  In particular,
\textcite{stengel-prb06} showed how to modify the above expressions in
a way that renders the spread functional strictly invariant under BZ
folding.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Localization procedure}
\label{sec:local}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In order to minimize the localization functional, we 
consider the first-order change of the spread functional
$\Omega$ arising from an infinitesimal gauge transformation
$U_{mn}^{\bf(k)}=\delta_{mn}+dW_{mn}^{\bf(k)}\,,$
where $dW$ is an infinitesimal anti-Hermitian matrix,
$dW^\dagger=-dW$, so that
$\vert u_{n\bf k}\rangle \;\rightarrow\; \vert u_{n\bf k}\rangle
+ \sum_m dW_{mn}^{\bf(k)} \, \vert u_{m\bf k}\rangle\;.$
%
We use the convention
%
\begin{equation}
\left({\frac{d\Omega}{dW}}\right)_{nm}= {\frac{d\Omega}{dW_{mn}}}
\label{eq:dfdw}
\end{equation}
%
(note the reversal of indices)
and introduce $\cal A$ and $\cal S$ as the superoperators
${\cal A}[B]=(B-B^\dagger)/2$ and ${\cal S}[B]=(B+B^\dagger)/2i$.
Defining
%
%
\begin{equation}
q_n^{\bf(k,b)} = {\rm Im}\,\ln M_{nn}^{\bf(k,b)}
+ {\bf b\cdot\bar{r}}_n \,,
\label{eq:qdef}
\end{equation}
%
\begin{equation}
R_{mn}^{\bf(k,b)} = M_{mn}^{\bf(k,b)} M_{nn}^{\bf(k,b)*} \,,
\label{eq:Rdef}
\end{equation}
%
\begin{equation}
T_{mn}^{\bf(k,b)} = 
\frac{M_{mn}^{\bf(k,b)}}{M_{nn}^{\bf(k,b)}} \; q_n^{\bf(k,b)} \,,
\label{eq:Tdef}
\end{equation}
%
and referring to \textcite{marzari-prb97} for the details, we arrive
at the explicit expression for the gradient $G^{\bf(k)}=d\Omega/dW^{\bf(k)}$
of the spread functional $\Omega$ as
%
\begin{equation}
G^{\bf(k)} = 4 \sum_{\bf b}\, w_b
\left( {\cal A}[R^{\bf(k,b)}] - {\cal S}[T^{\bf(k,b)}]
\right) \;\;.
\label{eq:grado}
\end{equation}
%
This gradient is used to drive the evolution of the
$U_{mn}^{({\bf k})}$ (and, implicitly, of the
$\vert\,{{\bf R}n}\,\rangle$ of \equ{wannier}) towards the minimum
of $\Omega$. A simple steepest-descent implementation, for example,
takes small finite steps in the direction
opposite to the gradient $G$ until a minimum is reached.

For details of the minimization strategies and the enforcement of
unitarity during the search, the reader is referred to
\textcite{mostofi-cpc08}.
%
We should like to point out here, however, that most of the
operations can be performed using inexpensive matrix algebra on
small matrices.  The most computationally demanding parts of the
procedure are typically the calculation of the self-consistent
Bloch orbitals $\vert u_{n\bf k}^{(0)}\rangle$ in the first place,
and then the computation of a set of overlap matrices
%
\begin{equation}
\label{eq:overlapmv}
M_{mn}^{(0)\bf(k,b)} = \langle u_{m\bf k}^{(0)}
\vert u_{n,\bf k+b}^{(0)}\rangle
\end{equation}
%
that are constructed once and for all from the $\vert u_{n\bf
  k}^{(0)}\rangle$.  After every update of the unitary matrices
$U^{\bf(k)}$, the overlap matrices are updated with inexpensive matrix
algebra
\begin{equation}
M^{\bf(k,b)} = U^{\bf(k)\dagger} \, M^{(0)\bf(k,b)} \, U^{\bf(k+b)}
\label{eq:newM}
\end{equation}
%
without any need to access the Bloch wavefunctions themselves.
This not only makes the algorithm computationally fast and efficient,
but also makes it independent of the basis used to represent the
Bloch functions.  That is, any electronic-structure code package
capable of providing the set of overlap matrices $M^{\bf(k,b)}$
can easily be interfaced to a common Wannier maximal-localization code.

It should be stressed that the localization functional sometimes
displays multiple local minima in addition to the desired global
minimum.  These false minima typically correspond to the formation
of topological defects (e.g., ``vortices'')  in the gauge field in
discrete $\kk$ space, and the resulting WFs are
poorly localized.  The best way to avoid these local
minima is to initialize the procedure with a gauge choice that is
already fairly smooth.  For this purpose, the projection method
already described in Sec.~\ref{sec:projection} is ideal. Local minima
can also arise when the branch cuts for the complex logarithms in
\equ{r} and \equ{rr} are inconsistent, i.e. when at any given {\bf
  k}-point the contributions from different {\bf b}-vectors differ by
random amounts of 2$\pi$ in the logarithm. A locally appropriate
choice of branch cuts can always be performed, and so this problem is
less severe than the topological problem alluded to before.
Therefore, the minimization procedure is initialized by first using the
projection method to generate a set of reference Bloch orbitals, as
discussed more fully in \textcite{marzari-prb97} and
\textcite{mostofi-cpc08}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{The limit of isolated systems or large supercells}
\label{sec:isolated}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The formulation introduced above can be significantly simplified in
two important and related cases, which merit a separate
discussion. The first is the case of open boundary conditions: this is
the most appropriate choice for treating finite, isolated systems
(e.g., molecules and clusters) using localized basis sets, and is a
common approach in quantum chemistry.  The localization procedure can
then be entirely recast in real space, and corresponds exactly to
determining Foster-Boys localized orbitals.  The second is the case of
systems that can be described using very large periodic
supercells. This is the most appropriate strategy for non-periodic
bulk systems, such as amorphous solids or liquids (see
Fig.~\ref{fig:a-si} for a paradigmatic example), but obviously
includes also periodic systems with large unit cells.  In this
approach, the Brillouin zone is considered to be sufficiently small
that integrations over ${\bf k}$-vectors can be approximated with a
single ${\bf k}$-point (usually at the $\Gamma$ point, i.e., the
origin in reciprocal space).  The localization procedure in this
second case is based on the procedure for periodic boundary conditions
described above, but with some notable simplifications.  Isolated
systems can also be artificially repeated and treated using the
supercell approach, although care may be needed in dealing with the
long-range electrostatics if the isolated entities are charged or have
significant dipole or multipolar character \cite{makov-prb95,
  dabo-prb08}.
\begin{figure}
\begin{center}
\includegraphics[width=6.0cm]{figs/a-si}
\end{center}
\caption{\label{fig:a-si} (Color online) MLWFs
in amorphous Si, either around distorted but fourfold coordinated atoms,
or in the presence of a fivefold defect.
Adapted from \textcite{fornari-cms01}.}
\end{figure}

\subsubsection{Real-space formulation for isolated systems}

For an isolated system, described with open boundary conditions, all
orbitals are localized to begin with, and further localization 
is achieved via unitary transformations within this set.
We adopt a simplified notation
$\vert{\bf R}n\rangle\rightarrow\vert w_i\rangle$ to refer
to the localized orbitals of the isolated system
that will become maximally localized. 
We decompose again the localization functional
$\Omega=\sum_i [\langle r^2\rangle_i-\bar{\bf r}_i^2]$ into a
term $\Omega_{\rm I}=\sum_\alpha{\rm tr}\,[Pr_\alpha
Qr_\alpha]$ (where $P=\sum_i\vert w_i\rangle\langle w_i\vert$,
$Q={\bf 1}-P$, and `$\rm tr$' refers to a sum over all the states $w_i$) that is invariant under unitary rotations, and 
a remainder 
$\widetilde\Omega = \sum_\alpha \sum_{i\ne j} \, \vert\langle w_i
\vert r_\alpha \vert w_j \rangle \vert^2$ that needs to be minimized.
Defining the matrices $X_{ij}=\langle w_i \vert x \vert w_j \rangle$,
$X_{{\rm D},ij}=X_{ij}\,\delta_{ij}$, $X^\prime=X-X_{\rm D}$ (and
similarly for $Y$ and $Z$), $\widetilde\Omega$ can be rewritten as
%
\begin{equation}
\label{eq:omegatil}
\widetilde\Omega = 
{\rm tr}\,[X^{\prime\,2}+Y^{\prime\,2}+Z^{\prime\,2}] \;\;.
\end{equation}
%
If $X$, $Y$, and $Z$ could be simultaneously diagonalized, then
$\widetilde\Omega$ would be minimized to zero (leaving only
the invariant part). 
This is straighforward in one dimension, but is not generally possible
in higher dimensions.
The general solution to the three-dimensional problem consists instead
in the optimal, but approximate, simultaneous co-diagonalization of
the three Hermitian matrices $X$, $Y$, and $Z$ by a single unitary
transformation that minimizes the numerical value of the localization
functional.  Although a formal solution to this problem is missing,
implementing a numerical procedure (e.g., by steepest-descent or
conjugate-gradients minimization) is fairly straightforward. It should
be noted that the problem of simultaneous co-diagonalizaton arises
also in the context of multivariate analysis \cite{flury-sjssc86} and
signal processing \cite{cardoso-sjmaa96}, and has been recently
revisited in relation with the present localization approach
\cite{gygi-cpc03} (see also Sec.~IIIA in \textcite{berghold-prb00}).

To proceed further, we write
%
\begin{equation}
d\Omega = 2\,{\rm tr}\,[X^{\prime}dX+Y^{\prime}dY+Z^{\prime}dZ] \;\;,
\label{eq:dotr}
\end{equation}
exploiting the fact that ${\rm tr}\,[X^{\prime}X_{\rm D}]=0$, and similarly for
$Y$ and $Z$.
%
We then consider an infinitesimal unitary transformation $\vert
w_i\rangle \;\rightarrow\; \vert w_i\rangle +\sum_j dW_{ji} \vert
w_j\rangle$ (where $dW$ is antihermitian),
from which $dX=[X,dW]$, and similary for $Y$, $Z$.  Inserting in \equ{dotr}
and using ${\rm tr}\,[A[B,C]]=\,{\rm tr}\,[C[A,B]]$ and
$[X^{\prime},X]=[X^{\prime},X_{\rm D}]$, we obtain
$d\Omega={\rm tr}\,[dW\,G]$ where the gradient of the spread functional
$G={d\Omega/dW}$ is given by
%
\begin{equation}
G= 2\,\Bigl\{ [X^{\prime},X_{\rm D}] + [Y^{\prime},Y_{\rm D}]
+ [Z^{\prime},Z_{\rm D}]\,dW \Bigr\} \;\;.
\label{eq:Greal}
\end{equation}
%
The minimization can then be carried out using a
procedure similar to that outlined above for
periodic boundary conditions.
Last, we note that minimizing 
$\widetilde\Omega$ is equivalent to maximizing 
${\rm tr}\,[X^2_{\rm D}+Y^2_{\rm D}+Z^2_{\rm D}]$, since
${\rm tr}\,[X^{\prime}X_{\rm D}]={\rm tr}\,[Y^{\prime}Y_{\rm D}]={\rm tr}\,[Y^{\prime}Y_{\rm D}]=0$.

\subsubsection{$\Gamma$-point formulation for large supercells}
\label{sec:gpoint}

A similar formulation applies in reciprocal space when 
dealing with isolated or very large systems
in periodic boundary conditions,
i.e., whenever it becomes
appropriate to sample the wavefunctions only at the $\Gamma$-point of 
the Brillouin zone.
For simplicity, we start with the case of a simple cubic lattice
of spacing $L$, and 
define the matrices $\cal{X}$, $\cal{Y}$, and $\cal{Z}$ as
\begin{equation}
{\cal{X}}_{mn}=\langle  w_m | e^{-2\pi ix/L} | w_n \rangle
\end{equation}
and its periodic
permutations (the extension to supercells of arbitrary symmetry has been derived by
\textcite{silvestrelli-prb99}).\footnote{We point out that the definition of the 
$\cal{X}, \cal{Y}, \cal{Z}$ matrices
for extended systems, now common in the literature, is different from the one used in the previous 
subsection for the $X, Y, Z$ matrices for isolated system.
We preserved these two notations for the sake of consistency with published work, 
at the cost of making less evident the close similarities that exist
between the two minimization algorithms.}
The coordinate $x_n$ of the $n$-th WF
center (WFC) can then be obtained from
\begin{equation}
x_n = -{\frac{L}{{2\pi}}}\,{\rm Im}\; {\rm ln} 
\langle  w_n | e^{-i{\frac{2\pi}{L}} x} | w_n \rangle 
= -{\frac{L}{{2\pi}}}\,{\rm Im}\; {\rm ln}\; {\cal{X}}_{nn} \;,
\label{eq:rcenter}
\end{equation}
and similarly for $y_n$ and $z_n$.
\equa{rcenter} has been shown by \textcite{resta-prl98}
to be the correct definition of the expectation value of 
the position operator for a system with periodic boundary conditions,
and had been introduced several years ago to deal with the 
problem of determining the average position of a single electronic 
orbital
in a periodic supercell \cite{selloni-prl87} (the
above definition becomes self-evident in the limit where $w_n$ tends
to a Dirac delta function).

Following the derivation of the previous subsection, or of
\textcite{silvestrelli-ssc98}, it can be shown that the
maximum-localization criterion is equivalent to maximizing the functional
\begin{equation}
\Xi = \sum_{n=1}^\nwann \left(|{\cal{X}}_{nn}|^2+|{\cal{Y}}_{nn}|^2+|{\cal{Z}}_{nn}|^2\right)\;.
\label{eq:ximax}
\end{equation}
The first term of the gradient $d\Xi/dA_{mn}$ is given by
$[{\cal{X}}_{nm}({\cal{X}}_{nn}^*-{\cal{X}}_{mm}^*)-{\cal{X}}_{mn}^*({\cal{X}}_{mm}-{\cal{X}}_{nn})]$, and 
similarly
for the second and third terms.  Again, once the gradient 
is determined, minimization can be performed using a steepest-descent
or conjugate-gradients algorithm; as always, the computational cost
of the localization procedure is small, given that it involves only
small matrices of dimension $\nbw \times \nbw$, 
once the scalar products needed to construct the initial
${\cal{X}}^{(0)}$, ${\cal{Y}}^{(0)}$ and ${\cal{Z}}^{(0)}$ have been calculated, which takes
an effort of order $\nbw^2*N_{\rm basis}$.
We note that in the 
limit of a single $\kk$ point the distinction between Bloch orbitals
and WFs becomes irrelevant, since no Fourier transform
from ${\bf k}$ to ${\bf R}$ is involved in the transformation
\equ{wannier}; rather, we want to find the optimal unitary
matrix that rotates the ground-state self-consistent orbitals
into their maximally localized representation, making this problem exactly equivalent to the one
of isolated systems.
%

Before concluding, we note that care should be taken when comparing
the spreads of MLWFs calculated in supercells of different sizes.
The Wannier centers and the general shape of the MLWFs often
converge rapidly as the cell size is increased, but
even for the ideal case of an isolated molecule,
the total spread $\Omega$ often displays much slower convergence.
This behavior derives from the finite-difference representation of the
invariant part $\Omega_{\,\rm I}$ of the localization functional
(essentially, a second derivative); while $\Omega_{\,\rm I}$
does not enter into the localization procedure, it does contribute
to the spread, and in fact usually represents the largest term.  This
slow convergence was noted by \textcite{marzari-prb97}
when commenting on the convergence properties of $\Omega$
with respect to the spacing of the Monkhorst-Pack mesh, and
has been studied in detail by others \cite{umari-prb03,stengel-prb06}.
For isolated systems in a supercell, this problem can be avoided
simply by using a very large $L$ in \equ{rcenter}, since in practice
the integrals only need to be calculated in the small region where
the orbitals are non-zero \cite{wu-phd04}. For extended bulk systems,
this convergence problem can be ameliorated significantly by
calculating the position operator using real-space integrals
\cite{lee-prl05,stengel-prb06,lee-phd06}.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Reality, uniqueness, and exponential localization}
\label{sec:wfprop}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In calculations on ordinary non-magnetic insulators without spin-orbit
interaction, in practice one always seems to find the MLWFs to
be real (more precisely, each can be chosen real by trivially removing
an overall phase).  In the $\Gamma$-sampling case this happens
whenever the starting orbitals are chosen to be real (as is always the
case in the presence of time-reversal symmetry), since the
localization gradient becomes real, and no symmetry breaking towards
complex wavefunctions can take place.  More often, however, the MLWFs
are built from Bloch functions on a mesh of ${\bf k}$-points where the
Bloch orbitals are in general complex, and such result is not so
obvious - \textcite{brouder-prl07} show that it holds true when there
  is a unique minimum in the localization functional.

  Another related question has to do with the existence of competing
  local minima of the localization functional in the space of the
  unitary matrices, as already mentioned in Sec.~\ref{sec:local}.
  This difficulty never seems to arise when the $\Gamma$-sampling
  formulation of Sec.~\ref{sec:gpoint} is used, provided only six
  neighbours are used for the reciprocal-space gradients; in this case
  the functional appears to be well-behaved, with a single basin of
  attraction leading to the global minimum.
%
  When using a $\kk$-point mesh, however, it is quite common to find
  several false local minima competing with the desired global
  minimum.  We attribute this to the complex and random phase
  relationships that can occur between orbitals in the
  three-dimensional topology of the Brillouin zone.  Typically we find
  the MLWFs to be intrinsically complex at the false local minima,
  while they are real at the global minimum, providing a hint for
  identifying the desired minimum.  However, it is much better to
  avoid the multiple-minimum problem if possible, and indeed this is
  one of the primary motivations for using the projection scheme of
  Sec.~\ref{sec:projection} to initialize the minimization.  When this
  is done, the minimization always seems to proceed to the desired
  global minimum, even with a rather poor choice of initial trial
  functions.

  Last, smoothness in reciprocal space translates in faster decay in
  direct space, and asymptotic expansions of Fourier transforms have
  been provided by \textcite{duffin-dmj53,duffin-dmj60}, with well-defined
  relations between the rank of the points of non-analyticity and
  their consequences on asymptotic decay.  
    In the systems considered it can be expected that WFs
  could be chosen to be exponentially localized.\footnote{One should
    note that MLWFs put all their emphasis on the quadratic spreads
    around their centers, and so it is conceivable that they could
    have a slower asymptotic decay than other methods that would focus
    on minimizing e.g. the decay in the tails.}

  Formal proofs that it is possible to find a sufficiently smooth
  gauge to insure the exponential localization are now available.  The
  exponential localization was proven for 1D systems already by
  \textcite{kohn-pr59}, and this result has been extended to the
  determination of the algebraic prefactors modulating the exponential
  decay by \textcite{he-prl01}.  A proof of the exponential
  localization of the projection operator for a band manifold in 3D
  was provided by \textcite{cloizeaux-pr64a,cloizeaux-pr64b}, and this
  was finally extended to the exponential localization of the WFs
  themselves in 3D by \textcite{brouder-prl07}, following a proof of
  the triviality of Bloch bundles \cite{panati-ahp07}. Last, exponential localization
of Wannier functions can even be achieved in the
case of metals, provided that empty states are also
considered, to disentangle a manifold of states that 
is only partially filled, but smooth. This is the case, e.g., of the $p_z$ MLWFs
describing the $\pi/\pi^\star$ manifold of graphene (see Fig.~\ref{fig:graph_pz}).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Hybrid Wannier functions}
\label{sec:hybrid}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

It is sometimes useful to carry out the Wannier transform in one
spatial dimension only, leaving wavefunctions that are still
delocalized and Bloch-periodic in the remaining directions
\cite{sgiarovello-prb01}.  Such orbitals are usually denoted as
``hermaphrodite'' or ``hybrid'' WFs.  Explicitly,
Eq.~(\ref{eq:wannier}) is replaced by the hybrid WF definition
%
\begin{equation}
\vert\,{l,n\kk_\parallel}\,\rangle\;=\;\frac{c}{2\pi}
\;\int_0^{2\pi/c} \ket{\,\psi\bnk\,}
\,e^{-ilk_\perp c}\,dk_\perp \;,
\label{eq:hybrid}
\end{equation}
%
where $\kk_\parallel$ is the wavevector in the plane (delocalized
directions) and $k_\perp$, $l$, and $c$ are the wavevector, cell
index, and cell dimension in the direction of localization.  The 1D
Wannier construction can be done independently for each
$\kk_\parallel$ using direct (i.e., non-iterative) methods as
described in Sec.~IV\,C\,1 of \textcite{marzari-prb97}.

Such a construction has proved useful for a variety of purposes, from
verifying numerically exponential localization in 1 dimension, to
treating electric polarization or applied electric fields along a
specific spatial direction
\cite{giustino-prl03,giustino-prb05,wu-prl06,stengel-prb06,murray-prb09}
or for analyzing aspects of topological insulators
\cite{coh-prl09,soluyanov-prb11a,soluyanov-prb11b}.  Examples will be
discussed in Secs.~\ref{sec:layer} and \ref{sec:surface-bands}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Entangled bands}
\label{sec:entangled}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The methods described in the previous sections were designed with {\it
  isolated} groups of bands in mind, separated from all other bands by
finite gaps throughout the entire Brillouin zone.  However, in many
applications the bands of interest are not isolated. This can happen,
for example, when studying electron transport, which is governed by
the partially filled bands close to the Fermi level, or when dealing
with empty bands, such as the four low-lying antibonding bands of
tetrahedral semiconductors, which are attached to higher conduction
bands. Another case of interest is when a partially filled manifold is
to be downfolded
%
into a basis of WFs - e.g.  to construct model Hamiltonians. In all
these examples the desired bands lie within a limited energy range but
overlap and hybridize with other bands which extend further out in
energy.  We will refer to them as {\it entangled bands}.

The difficulty in treating entangled bands stems from the fact that it
is unclear exactly which states to choose to form $\nwann$ WFs,
particularly in those regions of $\kk$-space where the bands of interest
are hybridized with other unwanted bands.  Before a Wannier
localization procedure can be applied, some prescription is needed for
constructing $\nwann$ states per $\kk$-point from a linear combination
of the states in this larger manifold.  Once a suitable
$\nwann$-dimensional Bloch manifold has been identified at each $\kk$,
the same procedure described earlier for an isolated group of bands
can be used to generate localized WFs spanning that manifold.

The problem of computing well localized WFs starting from entangled
bands is thus broken down into two distinct steps, subspace selection
and gauge selection.  As we will see, the same guiding principle can
be used for both steps, namely, to achieve ``smoothness'' in
$\kk$-space.  In the subspace selection step a $\nwann$-dimensional
Bloch manifold which varies smoothly as function of $\kk$ is
constructed. In the gauge-selection step that subspace is represented
using a set of $\nwann$ Bloch functions which are themselves smooth
functions of $\kk$, such that the corresponding WFs are well
localized.

\subsubsection{Subspace selection via projection}
\label{sec:subspace_proj}

The projection technique discussed in Section \ref{sec:projection} can
be easily adapted to produce $\nwann$ smoothly-varying Bloch-like
states starting from a {\it larger} set of Bloch
bands~\cite{souza-prb01}.  The latter can be chosen, for example, as
the bands lying within a given energy window, or within a specified
range of band indices.  Their number $\nbands_\kk\geq \nwann $ is not
required to be constant throughout the BZ.

We start from a set of $\nwann$ localized trial orbitals $g_n(\bf r)$
and project each of them onto the space spanned by the chosen
eigenstates at each $\kk$,
%
\begin{equation}
  \vert\phi_{n\bf k}\rangle=\sum_{m=1}^{\nbands_\kk}\, \vert\psi_{m\bf k}\rangle\,\langle\psi_{m\bf k}\vert g_n\rangle\;.
\label{eq:disen_proj}
\end{equation}
%
This is identical to \equ{proj}, except for the fact that the overlap
matrix $(A_{\bf k})_{mn}=\langle\psi_{m\bf k}\vert g_n\rangle$ has
become rectangular with dimensions $\nbands_\kk \times\nwann$.  We
then orthonormalize the resulting $\nwann$ orbitals using \equ{symm},
to produce a set of $\nwann$ smoothly-varying Bloch-like states across
the BZ, 
%
\begin{equation}
\vert\tilde{\psi}_{n\bf k}\rangle = \sum_{m=1}^\nwann\, 
\vert\phi_{m\bf k}\rangle(S_{\bf k}^{-1/2})_{mn}\;.
\label{eq:disen_orth}
\end{equation}
%
As in \equ{symm}, $(S_{\bf k})_{mn}=\langle\phi_{m\bf
  k}\vert\phi_{n\bf k}\rangle =(A_\kk^\dagger
A_\kk^{\phantom{\dagger}})_{mn}$, but with rectangular $A_\kk$
matrices.

The above procedure achieves simultaneously the two goals of subspace
selection and gauge selection, although neither of them is performed
optimally.  The gauge selection can be further refined by minimizing
$\widetilde{\Omega}$ within the projected subspace. It is also
possible to refine iteratively the subspace selection itself, as will
be described in the next section. However, for many applications this
one-shot procedure is perfectly adequate, and in some cases it may
even be preferable to more sophisticated iterative approaches. For
example, it often results in ``symmetry-adapted'' WFs which inherit
the symmetry properties of the trial orbitals~\cite{ku-prl02}.

As an example, we plot in Fig.~\ref{fig:si_proj} the eight
disentangled bands obtained by projecting the band structure of
silicon, taken within an energy window that coincides with the entire
energy axis shown, onto eight atomic-like $sp^3$ hybrids. The
disentangled bands, generated using Wannier interpolation
(Sec.~\ref{sec:band-interpolation}), are shown as blue triangles,
along with the original first-principles bands (solid lines).  While
the overall agreement is quite good, significant deviations can be
seen wherever higher unoccupied and unwanted states possessing some
significant $sp^3$ character are admixed into the projected
manifold. This behavior can be avoided by forcing certain Bloch states
to be preserved identically in the projected manifold - we refer to
those as belonging to a frozen ``inner'' window, since this is often
the simplest procedure for selecting them.  The placement and range of
this inner window will depend on the problem at hand. For example, in
order to describe the low-energy physics for, e.g., transport
calculations, the inner window would typically include all states in a
desired range around the Fermi level.

We show as red circles in Fig.~\ref{fig:si_proj} the results obtained
by forcing the entire valence manifold to be preserved, leading to a
set of eight projected bands that reproduce exactly the four valence
bands, and follow quite closely the four low-lying conduction bands.
For the modifications to the projection algorithm required to enforce
this frozen inner window, we refer to Sec.~III.G of
\textcite{souza-prb01}.

\begin{figure}
\begin{center}
\includegraphics[width=8.7cm]{figs/silicon_band}
\end{center}
\caption{\label{fig:si_proj} (Color online) Solid black lines: band
  structure of bulk crystalline Si. Blue triangles: band
  structure for the subspace selected by projection onto atomic $sp^3$
  orbitals. Red circles: band structure for the subspace selected by
  projection onto atomic $sp^3$ orbitals and forcing the valence
  manifold to be reproduced exactly using the frozen window
  indicated.}
\end{figure}

Projection techniques can work very well, and an application of this
approach to graphene is shown in Fig.~\ref{fig:graph_pz}, where the
$\pi/\pi^\star$ manifold is disentangled with great accuracy by a
simple projection onto atomic $p_z$ orbitals, or the entire occupied
manifold together with $\pi/\pi^\star$ manifold is obtained by
projection onto atomic $p_z$ and $sp^2$ orbitals (one every other
atom, for the case of the $sp^2$ orbitals - albeit bond-centered $s$
orbitals would work equally well).

\begin{figure}
\begin{center}
\includegraphics[width=8cm]{figs/graphite_band}
\end{center}
\caption{\label{fig:graph_pz} (Color online) Solid black lines: band
  structure of graphene. Blue triangles: band structure for the
  subspace selected by projection onto atomic $p_z$ orbitals. Red
  circles: band structure for the subspace selected by projection onto
  atomic $p_z$ orbitals on each site and sp$^2$ orbitals on alternate
  sites, and using the frozen window indicated.  The lower panels show
  the MLWFs obtained from the standard localization procedure applied
  to these two projected manifolds.}
\end{figure}

Projection methods have been extensively used to study
strongly-correlated systems \cite{ku-prl02,anisimov-prb05a,ku-prl11},
in particular to identify a ``correlated subspace'' for LDA+U or DMFT
calculations, as will be discussed in more detail in
Sec.~\ref{sec:basis}. It has also been argued~\cite{ku-prl10} that
projected WFs provide a more appropriate basis for the study of
defects, as the the persuit of better localization in a MLWF scheme
risks defining the gauge differently for the defect WF as compared to
the bulk. Instead, a straightforward projection approach ensures the
similarity between the WF in the defect (supercell) and in the
pristine (primitive cell) calculations, and this has been exploited to
develop a scheme to unfold the band-structure of disordered supercells
into the Brillouin zone of the underlying primitive cell, allowing a
direct comparison with angular resolved photoemission spectroscopy
experiments~\cite{ku-prl10}.

\subsubsection{Subspace selection via optimal smoothness}
\label{sec:disentanglement}

The projection onto trial orbitals provides a simple and effective way
of extracting a smooth Bloch subspace starting from a larger set of
entangled bands. The reason for its success is easily understood: the
localization of the trial orbitals in real space leads to smoothness
in $\kk$-space.  In order to further refine the subspace selection
procedure, it is useful to introduce a precise measure of the
smoothness in $\kk$-space of a manifold of Bloch states. The search for
an optimally-smooth subspace can then be formulated as a minimization
problem, similar to the search for an optimally-smooth gauge.

As it turns out, smoothness in $\kk$ of a Bloch space is precisely
what is measured by the funtional $\Omega_{\rm I}$ introduced in
Sec.~\ref{sec:rsr}. We know from \equ{spread-decomp} that the
quadratic spread $\Omega$ of the WFs spanning a Bloch space of
dimension $\nwann$ comprises two positive-definite contributions, one
gauge-invariant ($\Omega_{\rm I}$), the other gauge-dependent
($\widetilde{\Omega}$).  Given such a Bloch space (e.g., an isolated
group of bands, or a group of bands previously disentangled via
projection), we have seen that the optimally smooth gauge can be found
by minimizing $\widetilde{\Omega}$ with respect to the unitary mixing
of states within that space.

From this perspective, the gauge-invariance of $\Omega_{\rm I}$ means
that this quantity is insensitive to the smoothness of the individual
Bloch states $\ket{\tilde{u}_{n\kk}}$ chosen to represent the Hilbert
space.  But considering that $\Omega_{\rm I}$ is a part of the spread
functional, it must describe smoothness in some other sense.  What
$\Omega_{\rm I}$ manages to capture is the {\it intrinsic} smoothness
of the underlying Hilbert space. This can be seen starting from the
discretized $\kk$-space expression for $\Omega_{\rm I}$, \equ{omega_i},
and noting that it can be written as
%
\beq
\label{eq:omega_i_spillage}
\Omega_{\rm I}=\frac{1}{N}\sum_{{\bf k},{\bf b}} \,
w_b \, T_{{\bf k},{\bf b}}
\eeq
%
with
%
\beq
\label{eq:spillage}
T_{{\bf k},{\bf b}}= {\rm Tr} [{P}_{\bf k}\,{Q}_{{\bf k}+{\bf b}}],
\eeq
%
where ${P}_{\bf k}=\sum_{n=1}^\nwann \, \ket{\tilde{u}_{n \bf
    k}}\bra{\tilde{u}_{n \bf k}}$ is the gauge-invariant projector
onto the Bloch subspace at $\kk$, ${Q}_{\bf k}={\bf 1}-{P}_{\bf k}$,
and ``Tr'' denotes the electronic trace over the full Hilbert
space. It is now evident that $T_{{\bf k},{\bf b}}$ measures the
degree of mismatch (or ``spillage'') between the neighboring Bloch
subspaces at $\kk$ and $\kk+{\bf b}$, vanishing when they are
identical, and that $\Omega_{\rm I}$ provides a BZ average of the
local subspace mismatch.

The optimized subspace selection procedure can now be formulated as
follows~\cite{souza-prb01}.  A set of $\nbands_\kk \geq \nwann$ Bloch
states is identified at each point on a uniform BZ grid, using, for
example, an energy window.  An iterative procedure is then used to
extract self-consistently across the BZ the $\nwann$-dimensional
subspaces having collectively the smallest possible value of
$\Omega_{\rm I}$ (typically the minimization starts from an initial
guess for the target subspaces given, e.g., by projection onto trial
orbitals).  Viewed as function of $\kk$, the Bloch subspace obtained
at the end of the minimization is ``optimally smooth'' in that it
changes as little as possible with $\kk$. The minimization algorithm
can be easily modified in order to preserve identically the Bloch
eigenstates within an inner energy window.

As in the case of the one-shot projection, the outcome of this
iterative procedure is a set of $\nwann$ Bloch-like states at each
$\kk$ which are linear combinations of the initial $\nbands_\kk$
eigenstates.  One important difference is that the resulting states
are not guaranteed to be individually smooth, and the minimization of
$\Omega_{\rm I}$ must therefore be followed by a gauge-selection step,
which is in every way identical to the one described earlier for
isolated groups of bands. Alternatively, it is possible to combine the
two steps, and minimize $\Omega=\Omega_{\rm I}+\widetilde{\Omega}$
simultaneously with respect to the choice of Hilbert subspace and the
choice of gauge~\cite{thygesen-prl05,thygesen-prb05}; this should lead
to the most-localized set of $\nwann$ WFs that can be constructed from
the initial $\nbands_\kk$ Bloch states.  In all three cases, the
entire process amounts to a linear transformation taking from
$\nbands_\kk$ initial eigenstates to $\nwann$ smooth Bloch-like
states,
%
\beq
\label{eq:wannierization-matrix}
\ket{\tilde{\psi}_{n\kk}}=\sum_{m=1}^{\nbands_\kk}\,\ket{\psi_{m\kk}}
V_{\kk,mn}.
\eeq
%
In the case of the projection procedure, the explicit expression for
the $\nbands_{\kk}\times\nwann$ matrix $V_\kk$ can be surmised from
\eqs{disen_proj}{disen_orth}.

Let us compare the one-shot projection and iterative procedures for
subspace selection, using crystalline copper as an example.  Suppose
we want to disentangle the five narrow $d$ bands from the wide $s$
band that crosses and hybridizes with them, to construct a set of
well-localized $d$-like WFs. The bands that result from projecting
onto five $d$-type atomic orbitals are shown as blue triangles in
Fig.~\ref{fig:cu_5}.  They follow very closely the first-principles
bands away from the $s$-$d$ hybridization regions, where the
interpolated bands remain narrow.

\begin{figure}
\begin{center}
\includegraphics[angle=0,width=7cm]{figs/copper_band}
\end{center}
\caption{\label{fig:cu_5} (Color online) Solid black lines: band
  structure of bulk crystalline Cu. Blue triangles: band structure
  for the subspace selected by projection onto atomic $3d$
  orbitals. Red circles: band structure for the subspace derived from
  the previous one, once the criterion of optimal smoothness has been
  applied.}
\end{figure}

The red circles show the results obtained using the iterative scheme
to extract an optimally-smooth five-dimensional manifold.  The maximal
``global smoothness of connection'' is achieved by keeping the five
$d$-like states and excluding the $s$-like state. This happens because
the smoothness criterion embodied by \eqs{omega_i_spillage}{spillage}
implies that the orbital character is preserved as much as possible
while traversing the BZ.  Inspection of the resulting MLWFs confirms
their atomic $d$-like character.  They are also significantly
more localized than the ones obtained by projection and the
corresponding disentangled bands are somewhat better behaved,
displaying less spurious oscillations in the hybridization regions.

In addition, there are cases where the flexbility of the minimization
algorithm leads to surprising optimal states whose symmetries would
not have been self-evident  in advance.
One case is shown in Fig.~\ref{fig:cu_7}.
Here we wish to construct a minimal Wannier basis for copper,
describing both the narrow $d$-like bands and the wide
free-electron-like band with which they hybridize. By choosing
different dimensions for the disentangled subspace, it was found
that the composite set of bands is faithfully represented by seven
MLWFs, of which five are the standard $d$-like orbitals, and the
remaining two are $s$-like orbitals centered at the tetrahedral
interstitial sites of the fcc structure. The latter arise from the
constructive interference between $sp^3$ orbitals that would be part
of the ideal $sp^3d^5$ basis set; in this case, bands up to 20~eV
above the Fermi energy can be meaningfully described with a minimal
basis set of seven orbitals, that would have been difficult to
identify using only educated guesses for the projections.

\begin{figure}
\begin{center}
\includegraphics[angle=0,width=7cm]{figs/rmp_j5}
\end{center}
\caption{\label{fig:cu_7} (Color online) Solid black lines: band
  structure of bulk crystalline Cu. Coloured lines: band structure
  for the subspace selected by optimal smoothness, and a target
  dimensionality of 7, giving rise to 5 atom-centered $d$-like MLWFs
  and two $s$-like MLWFs in the tetrahedral interstitials, shown
  below.  The colour coding represents the projection of the
  disentangled bands onto these MLWFs, smoothly varying from red
  (representing $s$-like interstitial MLWFs) to blue (atom-centered
  $d$-like MLWFs).}
\end{figure}

  The concept of a natural dimension for the disentangled manifold has
  been explored further by
  \textcite{thygesen-prl05,thygesen-prb05}. As illustrated in
  Fig.~\ref{fig:thyg}, they showed that by
  minimizing $\Omega=\Omega_{\rm I}+\widetilde{\Omega}$ for different
  choices of $\nwann$, one can find an optimal $\nwann$
  such that the resulting ``partially-occupied''
  Wannier functions are most localized (provided enough bands are used
  to capture the bonding and anti-bonding combinations of those
  atomic-like WFs).

\begin{figure}
\begin{center}
\includegraphics[angle=270,width=6cm]{figs/thyg}
\end{center}
\caption{\label{fig:thyg} 
  Plots of average $\Xi$ per Wannier function vs.~size of
  the Wannier space for several molecules.  $\Xi$
  is defined in \equ{ximax}; larger value indicates greater
  localization. $N_{\rm occ}$ is the number of occupied
  states, $J$ is the size of the Wannier space, and $N$ is
  the total number of states included in the DFT calculation.
  Adapted from \textcite{thygesen-prl05}.}
\end{figure}


%On the other hand, both for band-interpolation and to construct model
%Hamiltonians, it is often required to have WFs representing the
%disentangled subspace that are as localized as possible. Such
%criterion corresponds to requiring that the disentangled subspace be
%as smooth as possible; also, heuristically, smoothness implies that
%the character of bands is preserved as much as possible while the
%entire Brillouin zone is spanned. Such concepts will be formally
%defined in the following, and used to identify a criterion and an
%algorithm that allows to disentangle target manifolds.
%
%\nmm{we'll need to put a note later on mentioning
%    how this "fails" when we start to be high in energy, where there
%    are a lot of spaghetti bands, and projection might be more robust}

%For definiteness, let us suppose we want to disentangle the five
%narrow $d$ bands of copper from the wide $s$ band which crosses them
%and construct a set of well-localized WFs associated with the $d$
%bands.  As mentioned in the previous subsection, the difficulty arises
%because there are regions of $\kk$-space where all six bands are close
%together, so that as a result of hybridization the distinction between
%$d$-band and $s$-band levels is not meaningful.  In order to separate
%the two manifolds, we rely on a numerical definition of smoothness and
%differentiability, and we extract a five-dimensional manifold by
%requiring that its change of character be as small as possible as we
%span the entire Brillouin zone. In a finite-sampling perspective, this
%is equivalent to requiring that there be maximal overlap between the
%disentangled states at one point in the Monkhorst-Pack mesh, and its
%closest neighbors.

%Such approach is best explained describing the procedure that we
%follow.  We again select an energy window or a group of bands that
%encompasses all the manifold of interest, and at each $\kk$-point we
%want to select $\nwann$ states that are a linear combination of the
%equal or larger number $\nbands_{\bf k}$ of bands under consideration.
%We denote by ${\cal F}({\bf k})$ the $\nbands_{\bf k}$-dimensional
%Hilbert space spanned by the eigenstates $|u_{n \bf k}\rangle$ within
%the window. Then, at any ${\bf k}$, if $\nbands_{\bf k} > \nwann$ we
%extract the $\nwann$-dimensional subspace ${\cal S}({\bf k}) \subseteq
%{\cal F}({\bf k})$ that, among all possible $\nwann$-dimensional
%subspaces of ${\cal F}({\bf k})$, leads to the smallest $\Omega_{\rm
%  I}$ [Eq.~(\ref{eq:omega_i})] (note that if $\nbands_{\bf k}=\nwann$
%there is nothing to do).  The choice of minimizing $\Omega_{\rm I}$ is
%easily explained, since $\Omega_{\rm I}$ measures the change of
%character of the states across the Brillouin zone. Indeed,
%Eq.~(\ref{eq:omega_i}) shows that $\Omega_{\rm I}$ is small whenever
%${| \langle \tilde{u}_{n \bf k} | \tilde{u}_{m,{\bf k}+{\bf b}}
%  \rangle |}^2$, the square of the magnitude of the overlap between
%states at nearby $\kk$-points, is large. Thus by minimizing $\Omega_{\rm
%  I}$ we are choosing self-consistently at every ${\bf k}$ the
%subspace ${\cal S}({\bf k})$ that changes as little as possible with
%${\bf k}$, i.e., has minimum ``spillage'' or mismatch with neighboring
%subspaces. In the present example this maximal ``global smoothness of
%connection'' will be achieved by keeping the five $d$-like states and
%excluding the $s$-like state.

%What is meant by spillage becomes clear once we rewrite
%Eq.~(\ref{eq:omega_i}) for $\Omega_{\rm I}(\{\tilde{u}_{n \bf k}\})$ as


\subsubsection{Iterative minimization of $\Omega_{\rm I}$}

The minimization of $\Omega_{\rm I}$ inside an energy window is
conveniently done using an algebraic algorithm~\cite{souza-prb01}. The
stationarity condition $\delta \Omega_{\rm I}(\{\tilde{u}_{n \bf
  k}\})=0$, subject to orthonormality constraints, is equivalent to
solving the set of eigenvalue equations
%
\begin{equation}
\label{eq:stationary_cond}
\Big[
  \sum_{\bf b} \, w_b \, {P}_{{\bf k}+{\bf b}}
\Big]
\left| \left. \tilde{u}_{n{\bf k}} \right\rangle  \right.=
\lambda_{n \bf k} \left| \left. \tilde{u}_{n{\bf k}} \right\rangle  \right..
\end{equation}
%
Clearly these equations, one for each $\kk$-point, are coupled, so that
the problem has to be solved self-consistently throughout the
Brillouin zone.  This can be done using an iterative procedure: on the
$i$-th iteration go through all the $\kk$-points in the grid, and for
each of them find $\nwann$ orthonormal states
$\left|\left.\tilde{u}_{n \bf k}^{(i)}\right\rangle \right.$, defining
a subspace
%${\cal S}^{(i)}({\bf k}) \subseteq {\cal F}({\bf k})$ such that the 
whose spillage over the neighboring subspaces 
%${\cal S}^{(i-1)}({\bf k}+{\bf b})$
from the previous iteration is as small as possible.  At each step the
set of equations
%
\begin{equation}
\label{eq:stationary_cond_iter}
\Big[
  \sum_{\bf b} \, w_b \, {P}_{{\bf k}+{\bf b}}^{(i-1)}
\Big]
\left| \left. \tilde{u}_{n{\bf k}}^{(i)} \right\rangle  \right.=
\lambda_{n \bf k}^{(i)} \left| \left. 
\tilde{u}_{n{\bf k}}^{(i)} \right\rangle  \right.
\end{equation}
%
is solved, and the $\nwann$ eigenvectors 
with largest eigenvalues are selected. That choice ensures that at
self-consistency the stationary point corresponds to the absolute
minimum of $\Omega_{\rm I}$.

In practice \equ{stationary_cond_iter} is solved in the basis of
the original $\nbands_{\bf k}$ Bloch eigenstates $|u_{n \bf k}\rangle$
inside the energy window.  Each iteration then amounts to
diagonalizing the following $\nbands_{\bf k} \times \nbands_{\bf k}$
Hermitian matrix at every ${\bf k}$:
%
\begin{equation}
\label{eq:z_matrix}
Z_{mn}^{(i)}({\bf k})=
\Big<  u_{m \bf k}
\Big| \sum_{\bf b} w_b 
  {\left[ {P}^{(i-1)}_{{\bf k}+{\bf b}} \right]}_{\rm in}
\Big| u_{n \bf k} \Big> .
\end{equation}
%
%
Since these are small matrices,
each step of the iterative procedure is computationally inexpensive.

As a final comment, we mention that 
in the case of degeneracies or quasi degeneracies in the
spreads of orbitals centered on the same site, the localization
algorithm will perform rather arbitrary mixing of these (as can be the
case, e.g., for the $d$ or $f$ electrons of a transition-metal ion, or
for its $t_2g$ or $e_g$ groups). A solution to this problem is to
diagonalize an operator with the desired symmetry in the basis of the
Wannier functions that have been obtained (see
\textcite{posternak-prb02} for the example of the $d$ orbitals in
MnO).

%We contrast in Fig.~\ref{fig:cu_5} the exact band structure of copper
%with that of 5-dimensional manifolds obtained by projection or by
%optimal smoothness - it can be seen that \hl{even in this
%  high-symmetry case}
%
%\ism{Again high symmetry is mentioned, but I don't quite follow.}
%
%the additional criterion of minimum spillage provides a more faithful
%representation of the bands. In addition, there are cases where the
%flexbility of the algorithms, or the lack of self-evident symmetries,
%provide firm guidance in selecting the optimal states - one case is
%shown in Fig.~\ref{fig:cu_7}, where the band structure of copper is
%faithfully represented with seven MLWFs, of which five are the
%standard $d$-like orbitals, but the remaining two are $s$-like
%\textcolor{blue}{orbitals centered at the tetrahedral interstitial
%  sites of the fcc structure}. These latter arise from the
%constructive interference of the $sp^3$ orbitals that would be part of
%the ideal basis set for copper; in this case, bands up to 20~eV above
%the Fermi energy can be meaningfully described with a minimal basis
%set of seven orbitals, that would have been difficult to identify
%using only educated guesses for the projections.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Many-body generalizations}
\label{sec:manybody}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The concept of WFs is closely tied to the framework
of single-particle Hamiltonians.  Only in this case can we
define $\nbw$ occupied single-particle Bloch functions at each
wavevector $\kk$ and treat all $\nbw$ of them on an equal footing,
allowing for invariance with respect to unitary mixing among them.
Once the two-particle electron-electron interaction is formally
included in the Hamiltonian, the many-body wavefunction cannot be
reduced to any simple form allowing for the construction of
WFs in the usual sense.

One approach is to consider the reduced one-particle density matrix
%
\beq
n(\rr,\rr') =\int \Psi^*(\rr,\rr_2...)\,
 \Psi(\rr',\rr_2,...)\, d\rr_2\,d\rr_3 ...
\eeq
%
for a many-body insulator.  Since $n(\rr,\rr')$ is periodic in the
sense of $n(\rr+\RR,\rr'+\RR)=n(\rr,\rr')$, its eigenvectors -- the
so-called ``natural orbitals'' -- have the form of Bloch functions
carrying a label $n,\kk$.  If the insulator is essentially a
correlated version of a band insulator having $\nbw$ bands, then at
each $\kk$ there will typically be $\nbw$ occupation eigenvalues
$\nu_{n\kk}$ that are close to unity, as well as some small ones that
correspond to the quantum fluctuations into conduction-band states.
If one focuses just on the subspace of one-particle states spanned by
the $\nbw$ valence-like natural orbitals, one can use them to
construct one-particle WFs following the methods described earlier, as
suggested by \textcite{koch-ssc01}.  However, while such an approach
may provide useful qualitative information, it cannot provide the
basis for any exact theory.  For example, the charge density, or
expectation value of any other one-particle operator, obtained by
tracing over these WFs will not match their exact many-body
counterparts.

A somewhat related approach, adopted by \textcite{hamann-prb09},
is to construct WFs out
of the quasiparticle states that appear in the GW approximation
\cite{aryasetiawan-ropp98}.  Such an approach will be described
more fully in Sec.~\ref{sec:gw-bands}.
Here again, this approach may give useful physical and chemical
intuition, but the one-electron quasiparticle wavefunctions do not
have the physical interpretation of occupied states, and charge
densities and other ground-state properties cannot be computed
quantitatively from them.

Finally, a more fundamentally exact framework for a many-body
generalization of the WF concept, introduced
in \textcite{souza-prb00b}, is to consider
a many-body system with twisted boundary conditions applied to
the many-body wavefunction in a supercell.  For example, consider
%
$M$ electrons in a supercell consisting of $M_1\times M_2\times M_3$
primitive cells, and impose the periodic boundary condition
%
\beq
\Psi_\qq(..., \rr_j+\RR, ...)=e^{i\qq\cdot\RR}\,\Psi_\qq(..., \rr_j, ...)
\label{eq:superper}
\eeq
%
for $j=1,...,M$,
where $\RR$ is a lattice vector {\it of the superlattice}.
One can then construct a single ``many-body WF''
in a manner similar to \equ{wanniertransform}, but with
$\kk\rightarrow\qq$ and $\ket{\psi_{n\kk}}\rightarrow\ket{\Psi_q}$
on the right side.  The resulting many-body WF is
a complex function of $3M$ electron coordinates, and as such it
is unwieldy to use in practice.  However, it is closely related
to Kohn's theory of the insulating state \cite{kohn-pr64},
and in principle it can be used to formulate many-body versions of
the theory of electric polarization and related quantities, as shall
be mentioned in Sec.~\ref{sec:pol-manybody}.

