% \section{RELATION TO OTHER LOCALIZED ORBITALS}
% \label{sec:locorb}


\subsection{Alternative localization criteria}
\label{sec:alt-loc}


As we have seen, WFs are inherently non-unique and, in practice, some
strategy is needed to determine the gauge. Two possible approaches
were already discussed in Sec.~\ref{sec:theory}, namely, projection
and maximal localization. The latter approach is conceptually more
satisfactory, as it does not depend on a particular choice of trial
orbitals. However, it still does not uniquely solve the problem of
choosing a gauge, as different localization criteria are possible and
there is, \emph{a priori}, no reason to choose one over another.

While MLWFs for extended
systems have been generated for the most part by minimizing the sum of
quadratic spreads, Eq. (17), a variety of other localization criteria
have been used over the years for molecular systems. We will briefly
survey and compare some of the best known schemes below. What they all
have in common is that the resulting localized molecular orbitals
(LMOs) $\phi_i(\rr)$ can be expressed as linear combinations of a set
of molecular eigenstates $\psi_i(\rr)$ (the ``canonical'' MOs),
typically chosen to be the occupied ones,
%
\beq
\phi_i({\rr})=\sum_{j=1}^{\nwann}U_{ji}\psi_j({\rr}) .
\label{eq:lmos}
\eeq
%
The choice of gauge then arises from minimizing or maximizing some
appropriate functional of the LMOs with respect to the coefficients
$U_{ij}$, under the constraint that the transformation (\ref{eq:lmos})
is unitary, which ensures the orthonormality of the resulting LMOs.


{\it The Foster-Boys criterion} (FB) \cite{boys-rmp60,foster-rmp60a,boys-66}.
This is essentially the same as the
Marzari-Vanderbilt criterion of minimizing the sum of the quadratic
spreads, except that the sum runs over the orbitals in the molecule,
rather than in one crystal cell,
%
\begin{equation}
\Omega_{\rm FB}=\sum_{i=1}^{\nwann} [\langle \phi_i|{\rr}^2|\phi_i\rangle -
\langle \phi_i|{\rr}|\phi_i\rangle ^2].
\label{eq:omega_boys}
\end{equation}
%
Interestingly, this criterion is equivalent to minimizing the
``self-extension'' of the orbitals \cite{boys-66},
%
\beq
\sum_{i=1}^{\nwann}\int d{\rr}_1 d{\rr}_2 \, |\phi_i({\rr}_1)|^2
 ({\rr}_1-{\rr}_2)^2 |\phi_i({\rr}_1)|^2
%\sum_{i}^{\nwann}\int d{\rr}_1 d{\rr}_2 \phi_i^{*}({\rr}_1)
%\phi_i^{*}({\rr}_2) ({\rr}_1-{\rr}_2)^2 \phi_i({\rr}_1)
%\phi_i({\rr}_2),
\label{eq:omega_boysa}
\eeq
and also to maximizing the sum of the squares of the distance between
the charge centers of the orbitals
\begin{equation}
\sum_{i,j=1}^{\nwann} |\langle \phi_i|{\rr}|\phi_i\rangle - \langle
\phi_j|{\rr}|\phi_j\rangle |^2 ,
\label{eq:omega_boys1}
\end{equation}
%
which is closely related to \equ{ximax} of Sec.~\ref{sec:gpoint}.
The relation between Eqs.~(\ref{eq:omega_boys})--(\ref{eq:omega_boys1}) is discussed in \textcite{boys-66}.

{\it The Edmiston-Ruedenberg criterion} (ER). Here localization is achieved by
maximizing the Coulomb self-interaction of the orbitals
\cite{edmiston-rmp63}
%
\begin{equation}
\Omega_{\rm ER}=
\sum_{i=1}^{\nwann}\int d{\rr}_1 d{\rr}_2 \, |\phi_i({\rr}_1)|^2
({\rr}_1-{\rr}_2)^{-1} |\phi_i({\rr}_1)|^2.
\label{eq:omega_er}
\end{equation}
%
From a computational point of view, the ER approach is more demanding 
(it scales as $\nwann^4$ versus $\nwann^2$ for FB),
but recently more efficient implementations have been developed
\cite{subotnik07}.

%
{\it The von Niessen criterion} (VN). The %idea
goal here is to maximize the density overlap of the
orbitals \cite{niessen-jcp72}
\begin{equation}
\Omega_{\rm VN}=
\sum_{i=1}^{\nwann}\int d{\rr}_1 d{\rr}_2 \, |\phi_i({\rr}_1)|^2
 \delta({\rr}_{1}-{\rr}_2) |\phi_i({\rr}_1)|^2.
\label{eq:omega_vn}
\end{equation}

{\it The Pipek-Mezey criterion} \cite{pipek-jcp89} (PM).
This approach differs from the
ones above in that it makes reference to some extrinsic objects. The
idea is to maximize the sum of the squares of the Mulliken atomic
charges \cite{mulliken-jcp55}. These are obtained with respect to a set
of atomic orbitals $\chi_{\mu}$ centered on atomic sites A. We define
\beq
\label{eq:PM}
P_{A}=\sum_{\mu\in A}\sum_{\nu} D_{\mu\nu}S_{\mu\nu}
\eeq
where the sum over $\mu$ involves all of the atomic states on atom
site A, $D_{\mu\nu}$ is the density matrix in the atomic basis, and
$S_{\mu\nu}=\inprod{\chi_{\mu}}{\chi_{\nu}}$ is the overlap operator.
The functional to be maximized is given by
%
\begin{equation}
\Omega_{\rm PM}=\sum_{i=1}^{\nwann} \sum_{A=1}^{N_A}
|\langle\phi_{i}|P_{A}|\phi_{i}\rangle |^{2}.
\label{eq:omega_pm}
\end{equation}
%
This is somewhat similar in spirit to the
projection scheme discussed in Sec.~\ref{sec:projection}, except that
it is not a one-shot procedure.

As indicated in \equ{lmos}, all of these schemes amount to unitary
transformations among the canonical MOs and, as such, they
give rise to representations of the electronic
  structure of the system that are equivalent to that provided by the
  original set of eigenstates.
For the purpose of providing chemical intuition, the usefulness of a
given scheme depends on how well it matches a particular
viewpoint of bonding. There have been few studies of the VN scheme, but
the FB, ER, and PM schemes have been extensively compared. In many
cases all three approaches lead to similar localized orbitals which
agree with the simple Lewis picture of bonding.
A notable
exception is for systems that exhibit both $\sigma$ and $\pi$
bonding. For a double bond, both the FB and ER schemes mix the
$\sigma$ and $\pi$ orbitals, giving rise to two bent ``banana bond''
(or $\tau$) orbitals \cite{pauling-book60} as shown in
Fig~\ref{fig:banana_bond}.  When applied to benzene,
both schemes give alternating $\sigma$ and $\tau$ orbitals. Trivially,
there are two equivalent sets of these orbitals, reminiscent of the
two Kekul\'{e} structures for benzene.  In contrast to the FB and ER
schemes, PM provides a clear separation of $\sigma$ and $\pi$
orbitals. For example, in benzene PM gives a framework of six $\sigma$
orbitals and three $\pi$ orbitals around the ring.

In some situations the FB and ER orbitals have been found to be quite
different. This has been observed when the bonding in a molecule can
be represented as two distinct resonant structures. The ER scheme
generally gives orbitals corresponding to one of the possible
structures, whilst the FB orbitals correspond a hybrid of the
structures.  An extreme example is CO$_2$ \cite{brown-jacs77}. In
agreement with the O=C=O Lewis picture, ER gives two lone pairs on
each oxygen and two $\tau$ orbitals between each carbon and oxygen. In
contrast, the FB scheme gives a single lone pair on each oxygen and
three highly polarized $\tau$ orbitals between the carbon and each
oxygen, as shown in Fig.~\ref{fig:banana_bond}.

\begin{figure}
\begin{center}
\includegraphics[width=7cm]{figs/rmp_ch3_fig1}
\end{center}
\caption{\label{fig:banana_bond} (Color online) (a) One of the two
  ``banana'' or $\tau$ orbitals in ethene (b) Centers of the MLWF in
  ethene (c) $\tau$ MLWF in CO$_2$ (d) Centers of the MLWF in
  CO$_2$. Atoms colors are: hydrogen (white), carbon (yellow), oxygen
  (red). MLWF centers are shown in black. MLWF are computed using the
  scheme of Marzari and Vanderbilt \cite{marzari-prb97} with a large
  vacuum supercell and a $\Gamma$-point sampling of the BZ.}
\end{figure}
MLWFs, which may be thought of
  as the solid-state equivalent of FB orbitals, have been
widely used to examine chemical bonding, as will be discussed in
detail in Sec.~\ref{sec:bonding}. The ER scheme has not been
extensively examined, an isolated exception being the work of
\textcite{miyake-prb08} who proposed a method to maximize the
Coulomb self-interaction of the orbitals in periodic systems (see also
Sec.~\ref{sec:mod-ham}).
They examined a range of bulk transition metals
and SrVO$_3$ and in each case found that the resulting WFs were
essentially identical to MLWFs. 
\subsection{Minimal-basis orbitals}
\label{sec:proj-down}

In the same way as was described for the Marzari-Vanderbilt scheme of
Sec.~\ref{sec:maxloc}, the alternative localization criteria
described above can be applied in a solid-state context to isolated
groups of bands. One is often interested in the more general situation
where the bands of interest are not isolated. The challenge then is to
generate a ``disentangled''
set of $\nwann$ localized WFs starting from some larger set
of bands.  Two procedures for doing so were discussed in
Sec.~\ref{sec:entangled}, namely, projection and iterative
minimization of the spread functional $\Omega$, not only with respect
to the choice of gauge, but also with respect to the choice of Hilbert
space.

 Here we discuss two further procedures which achieve the same goal by
different means. They have in common the fact that the resulting
orbitals constitute a minimal basis of atomic-like orbitals.

%\jym{The relation between the projection and
%partitioning/downfolding approach has been discussed by
%\textcite{solovyev-prb07} where under certain conditions the approaches
%are shown to be equivalent.  }
% We note that in general projection or downfolding schemes may
% produce orbitals which are neither orthonormal or localised according
% to a particular criteria. However, if this is
%required the orbitals could straight-forwardly be symmetrically
%orthonormalised, and/or used as input to one of the localization scheme
%in the previous section.
\subsubsection{Quasiatomic orbitals}

The quasiatomic orbitals (QOs)
scheme \cite{lu-prb04,chan-prb07,qian-prb08}
is a projection-based approach that
has the aim of extracting a minimal tight-binding basis
%
from an initial first-principles calculation.
%
%
In this regard it is similar in spirit to the Wannier interpolation
techniques discussion in Sec.~\ref{sec:winterp}.  Unlike WFs, however,
the quasiatomic orbitals form a {\it non-orthogonal} basis of {\it
  atom-centered} functions, each with specified angular
character. Their radial part, and hence their detailed local shape,
depends on the bonding environment.

As in the Wannier scheme, the first step is to construct a suitable
$\nwann$-dimensional subspace $\{\phi_{n}\}$ starting from a larger
set of $\nbands$ Bloch eigenstates $\{ \psi_{n}\}$. For simplicity
we consider a $\Gamma$-point only sampling of the BZ and hence omit the
$\kk$ index. 
%
The goal is to
identify a disentangled subspace
%that has, as much as possible, 
with
the desired atomic-orbital character, as specified by a
  given set of $J$ atomic orbitals (AOs) $\ket{A_i}$, where $i$ is a
composite site and angular-character index. One possible strategy
would be to employ the one-shot projection approach of
Sec.~\ref{sec:subspace_proj}, using the AOs as trial orbitals.
\textcite{lu-prb04} proposed a more optimized procedure, based on
maximizing the sum-over-square similarity measure\footnote{This is
  similar in spirit to the Pipek-Mazey procedure, \equ{PM}, but
  applied to subspace selection rather than gauge selection.}
%
\beq
\label{eq:similarity}
{\cal L}=\sum_{i=1}^\nwann\left|\left|
\sum_{n=1}^{\nwann}\ket{\phi_{n}}\inprod{\phi_{n}}{A_i}\right|\right|^2
\eeq
%
with respect to the orthonormal set $\{\phi_{n}\}$, expressed as
linear combinations of the original set $\{ \psi_{n}\}$. It is usually
required that a subset of $N\leq \nwann$ of the original
eigenstates are identically preserved (``frozen in'') in the
disentangled subspace, in which case the optimization is performed
with respect to the remaining $\nwann-N$ states $\phi_{n}$. $\nbands$
must be of sufficient size to capture all of the antibonding character
of the AOs.

In later work it was realised
\cite{qian-prb08} that the QOs can be constructed without explicit
calculation of the eigenstates outside the frozen window. The key insight
is to realise that the QOs will only have a contribution from the finite
subset of this basis spanned by the AOs (see \equ{qo} below).  
The component of the AOs projected onto
the $N$ states within the frozen window
$\ket{A_{i}^{\|}}$ is given by
\beq
 \ket{A_{i}^{\|}}=\sum_{n=1}^{N}  \ket{\psi_{n}}\inprod{\psi_{n}}{A_{i}}.
\eeq
The component of $\ket{A_{i}}$
projected onto the states outside the frozen window
%$\ket{A_{i}^{\bot}}$,  
can hence constructed directly using only the AOs and the states
within the frozen window, as
%
\beq
\ket{A_{i}^{\bot}}=\ket{A_{i}}-\ket{A_{i}^{\|}} .
\eeq
%
Using $\ket{A_{i}^{\bot}}$  as a basis, the set of $\{\phi_{n}\}$ which
maximize $\cal{L}$ can be obtained using linear algebra as
reported in \textcite{qian-prb08}. 

Once a subspace with the correct orbital character has been
identified, a basis of
quasiatomic orbitals 
%spanning that subspace 
can be
obtained by retaining the portion of the original AOs that ``lives''
on 
%the 
that
subspace,
%
\beq\label{eq:qo}
  \ket{Q_{i}}=\sum_{n=1}^\nwann\ket{\phi_{n}}\inprod{\phi_{n}}{A_{i}}.
  \eeq 
%
  In general the angular dependence of the resulting QOs will no
  longer that of pure spherical harmonics, but only approximately so.
%~\cite{qian-prb08}. [IS: citation repeated immediately]
In \textcite{qian-prb08}, QOs were obtained for simple metals and
semiconductors. Later applications have used the orbitals for the
study of quantum transport in nano-structures \cite{qian-prb10}.

\subsubsection{NMTO and Downfolding}
%
An alternative scheme for obtaining a minimal basis representation is the
perturbation approach introduced by \textcite{lowdin-jcp51}. Here the
general strategy is to partition a set of orbitals into an ``active'' set that is
intended to describe the states of interest, and a ``passive'' set that
will be integrated out.  Let us write the
Hamiltonian for the system in a block representation,
%
\beq \label{eq:lowdin_part}
H=\left(
  {\begin{array}{cc}
      H_{00} & 0  \\
      0 & H_{11} \\
 \end{array} } \right) +\left( {\begin{array}{cc}
 0 & V_{01}  \\
 V_{10} & 0  \\
\end{array} } \right) ,
\eeq
where $H_{00}$ ($H_{11}$) is the projection onto the active (passive)
subspace, and $V_{01}$ is the coupling between the two subspaces. An
eigenfunction can be written as the sum of its projections onto the
two subspaces $\ket{\psi}=\ket{\psi_0}+\ket{\psi_1}$.  This leads to
%
\beq
(H_{00}-\varepsilon) \ket{\psi_0} + V_{01}\ket{\psi_1}=0,
\eeq
\beq
V_{10}\ket{\psi_0}+(H_{11}-\varepsilon) \ket{\psi_1}=0,
\eeq
where $\varepsilon$ is the eigenvalue
  corresponding to $\ket{\psi}$.
Elimination of $\ket{\psi_1}$ gives an effective Hamiltonian for the
system which acts only on the active subspace:
%
\beq 
H^{\rm eff}_{00}(\varepsilon)=H_{00}-V_{01}(H_{11}-\varepsilon)^{-1}V_{10}.
\eeq
%
This apparent simplification has introduced an energy dependence
into the Hamiltonian. One practical way forward is to approximate this
%
as an energy-independent Hamiltonian $H^{\rm
  eff}_{00}(\varepsilon_0)$, choosing the reference energy
$\varepsilon_0$ to be the average energy of the states of
interest. This approach has been used to contruct tight-binding
Hamiltonians from full electronic structure calculations \cite{solovyev_prb04}.

A form of L\"owdin partitioning has
been widely used in connection with
the linear-muffin-tin-orbital (LMTO) method, particularly in its most
recent formulation, the $N^{\rm th}$-order muffin-tin-orbital (NMTO)
approach \cite{andersen-prb00,zurek-c05}.
In this context it is usually referred to as ``downfolding'',
although we note that some authors use this term to refer to any scheme to produce a
minimal basis-set representation.
%

In studies of complex materials there may be a significant number of
MTOs, typically one for each angular momentum state ($s$, $p$, $d$) on
every atomic site. One may wish to construct a minimimal basis to
describe states within a particular energy region, e.g., the occupied
states, or those crossing the Fermi level.  Let us assume we have of
basis orbitals (MTOs in this case) which we wish to partition into
``active'' and ``passive'' sets. Using the notation of
\equ{lowdin_part}, L\"owdin partitioning gives a set of
energy-dependent orbitals $\phi_0(\varepsilon,\rr)$ for the active
space according to~\cite{zurek-c05}
%
\beq
\phi_0(\varepsilon,\rr)=\phi_0(\rr)-\phi_1(\rr) (H_{11}-\varepsilon)^{-1}V_{10}.
\eeq
%
Taking into account the energy dependence, this reduced set of
orbitals spans the same space as the original full set of orbitals,
and can be seen to be the original orbitals of the
active set dressed by an energy-dependent linear combination of
orbitals from the passive set. In the NMTO scheme, the next step is
to form an energy-independent set of orbitals through an $n$th-order
polynomial fit to the energy dependence.

To give a specific example, accurate calculations on tetrahedral
semiconductors will require the inclusion of $d$ states in the MTO
basis, i.e., nine states ($s$, $p$, $d$) per site. However, it would
be desirable to construct a minimal basis to describe the valence and
lower conduction states with only four states ($s$, $p$) on each site
\cite{lambrecht_prb86}. We therefore designate $s$ and $p$ as
``active'' channels and $d$ as ``passive''.  Downfolding will result
in an MTO with either $s$ or $p$ character on a given site, with
inclusion of $d$ character on neighbouring sites. In other words, the
tail of the MTO is modified to ``fold in'' the character of the
passive orbitals.

In Fig.~\ref{fig:nmto_graph}, the band structure of graphite
calculated using a full $s$, $p$, and $d$ basis on each carbon atom is shown in
black \cite{zurek-c05}. The red bands are obtained by choosing $s$,
$p_x$, and $p_y$ states on every second
carbon atom as the active channels and downfolding all other
states. The energy mesh spans the energy range of the $sp^2$ bonding
states (shown on the right-hand-side of the band-structure plot in
Fig.~\ref{fig:nmto_graph}). For these bands the agreement with the
full calculation is perfect to within the resolution of the
figure.  Symmetric othonormalization of the three NMTOs gives the
familiar set of three symmetry related $\sigma$-bonding orbitals (compare
with the MLWF of graphene in Fig.~\ref{fig:graph_pz}).
%
\begin{figure}
\begin{center}
\includegraphics[width=8.5cm]{figs/rmp_nmto_graph}
\end{center}
\caption{\label{fig:nmto_graph} (Color online)
The band structure of graphite calculated with a full NMTO spd basis is
given in black. The red bands have been calculated with an s, px, and
py orbital on
every second carbon atom (shown above the band structure). Also shown is one of the sp2-bond
orbitals which arise by symmetrical orthonormalization of the s, px,
and py orbitals. The energy meshes used for each calculation
are given to the right of the band structure. From \textcite{zurek-c05}.}
\end{figure}
%
\textcite{lechermann-prb06} have compared MLWFs and NMTO orbitals for
a set of t$_{2g}$ states located
around the Fermi level SrVO$_3$. It was found that both schemes gave
essentially identical orbitals.


\subsection{Non-orthogonal orbitals and linear scaling}
\label{sec:linear-scaling}
%
In recent years there has been much progress in the development
of practical linear-scaling methods for electronic structure
calculations, that is,
methods in which the computational cost of the calculation grows only
linearly with the size of the system. The fundamental principle behind
such approaches is the fact that electronic structure is inherently
local~\cite{heine-ssp80}, or `nearsighted'~\cite{kohn-prl96}. This
manifests itself in the exponential localization of the WFs
in insulators 
%
(discussed in Sec.~\ref{sec:wfprop}) and, more generally, in the
localization properties of the single-particle density matrix
%
\beq
\rho(\rr,\rr') = \matel{\rr}{P}{\rr'}
= \ibz d\kk\; \sum_{n} f_{n\kk} \psi_{n\kk}(\rr)\psi^{\ast}_{n\kk}(\rr') , \label{eq:dm-bloch}
\eeq
%
where, following \textcite{janak-prb78}, the projection operator $P$ of
\equ{bandproj} has been generalized to the case of fractional
eigenstate occupancies $f_{n\kk}$. The quantity $\rho(\rr,\rr')$
has been shown to decay exponentially as $\exp(-\gamma|\rr-\rr'|)$
in insulators and semiconductors, where the exponent $\gamma$
can be heuristically related to the direct band gap of the
system~\cite{cloizeaux-pr64a, ismail-prl99, taraskin-prl02}.
%
It has also been shown to take the same form in metals at finite
temperature,\footnote{For metals at zero temperature, the
  discontinuity in occupancies as a function of $\kk$ results in the
  well-known \emph{algebraic} decay of the density matrix.} but with
$\gamma$ determined by
%
the ratio between the thermal energy $k_{\rm B}T$ and the Fermi
velocity~\cite{goedecker-prb98}.

Exponential localization may seem a surprising result given that the
Bloch eigenstates extend across the entire system.  Expressing the
density matrix in terms of WFs using \equ{wannier}, we find
%
\beq \rho(\rr,\rr') = \sum_{ij}\sum_{\RR\RR'}
w_{i\RR}(\rr)K_{ij}(\RR'-\RR)w^{\ast}_{j\RR'}(\rr'),
\label{eq:dm-wannier}
\eeq
%
where we have defined the 
%
\emph{density kernel}\footnote{This term
  was, to the best of our knowledge, first used by
  \textcite{mcweeny-romp60}.}
%
\beq K_{ij}(\RR) = \ibz d\kk\; e^{-i\kk\cdot\RR} 
\sum_{n} [U^{(\kk)\dagger}]_{in} f_{n\kk} [U^{(\kk)}]_{nj} ,
\label{eq:kij1}
\eeq
%
which reduces to $\delta_{ij}\delta_{\RR\0}$ in the case of a set of
fully occupied bands. 
%
We now see that the spatial localization of the density
matrix is closely linked to that of the Wannier functions
themselves. This locality is exploited in
linear-scaling methods by retaining an amount of information in the
density matrix that scales only linearly with system size.

Many different linear-scaling DFT approaches exist;
for comprehensive reviews the reader is referred to
\textcite{galli-cossms96},
%
\textcite{goedecker-romp99}, and \textcite{bowler-arxiv11}. Many of
them are based on the variational minimization of an energy
%
functional expressed either in terms of localized Wannier-like
orbitals or the density operator itself.  The common point between
these variational methods is that the idempotency of the density
operator or the orthogonality of the Wannier orbitals is not
imposed explicitly during the minimization procedure. Instead,
the energy functionals are constructed such that these properties
are satisfied automatically at the minimum, which coincides with
the true ground state of the system.

Many of these methods also make use of non-orthogonal
localized orbitals, referred to as ``support
functions''~\cite{hernandez-prb95} or ``non-orthogonal generalized
Wannier functions'' (NGWFs)~\cite{skylaris-prb02}, in contrast to
canonical WFs, which are orthogonal. 
The density matrix in \equ{dm-wannier} 
can be generalized so as to be represented in terms of a set
of non-orthogonal localized orbitals $\{\phi_{\alpha\RR}(\rr)\}$ and a
corresponding non-unitary (and, in general, non-square) transformation
matrix $M^{(\kk)}$, which take the place of $\{w_{i\RR}(\rr)\}$ and
$U^{(\kk)}$, respectively. 
Two main benefits arise from allowing non-orthogonality. First, it is
no longer necessary to enforce explicit orthogonality constraints on
the orbitals during the energy minimization
procedure~\cite{galli-prl92, mauri-prb93, ordejon-prb93, hierse-prb94,
  hernandez-prb96}.
Second, a non-orthogonal representation can be
more localized than an essentially equivalent orthogonal
one~\cite{anderson-prl68, he-prl01}. 
%%%%%%%%%%%%%%%%%%%%%%
% AAM: EXTRA DETAIL IN CASE NEEDED
%\beq
%\rho(\rr,\rr') = \sum_{\alpha\beta} \sum_{\RR\RR'}
%\phi_{\alpha\RR}(\rr)K^{\alpha\beta}(\RR-\RR')\phi^{\ast}_{\beta\RR'}(\rr'), 
%\label{dm-ngwf1}
%\eeq
%%
%in which 
%%
%\beq
%K^{\alpha\beta}(\RR-\RR') = \ibz \sum_{n} [M^{(\kk)\dagger}]^{\alpha}_{n} f_{n\kk} [U^{(\kk)}]^{\beta}_{n} 
%\: e^{i\kk\cdot(\RR-\RR')},
%\eeq
%AAM: END EXTRA DETAIL
%%%%%%%%%%%%%%%%%%%%%
In practice, linear-scaling methods target large systems, which means
that $\Gamma$-point only sampling of the BZ is usually sufficient. In
this case, the separable form for the density matrix simplifies to
%
\beq
\rho(\rr,\rr') = \sum_{\alpha\beta}
\phi_{\alpha}(\rr)K^{\alpha\beta}\phi^{\ast}_{\beta}(\rr'), 
\label{dm-ngwf2}
\eeq
%
where the density kernel is\footnote{It is worth noting that the
  non-orthogonality of the orbitals results in an important
  distinction between covariant and contravariant quantities, as
  denoted by
% upstairs and downstairs
  raised and lowered
  Greek indices~\cite{artacho-pra91, oregan-prb11}.}
%
\beq
K^{\alpha\beta} = \sum_{n} [M^{\dagger}]^{\alpha}_{n} f_{n} [M]^{\beta}_{n}.
\eeq
%

Minimization of an appropriate energy functional with respect to the
degrees of freedom present in the density matrix leads to
ground-state non-orthogonal orbitals that are very similar in
appearance to (orthogonal) MLWFs.
Fig.~\ref{fig:ngwfs} shows, for example, NGWFs on a Ni atom in bulk
NiO, obtained using the {\sc onetep} linear-scaling DFT
code~\cite{skylaris-jcp05}. A recent comparison of
  static polarizabilities for molecules, calculated using \equ{pol-real} with both
  MLWFs and NGWFs, demonstrates remarkable agreement between the
  two~\cite{oregan-tbc}.
\begin{figure}
\begin{center}
\includegraphics[width=5cm]{figs/ngwfs}
\end{center}
\caption{\label{fig:ngwfs}
(Color online) Isosurfaces of the set of nine non-orthogonal
generalized Wannier functions (NGWFs) on a nickel atom in NiO (shown 
centered on different, symmetrically equivalent, Ni atoms in the
lattice). The isosurface is set to half of the maximum for the $s$ and
$p$-like NGWFs and $10^{-3}$ times the maximum for the $d$-like NGWFs. 
Adapted from \textcite{oregan-prb11}.}
\end{figure}

\subsection{Other local representations}
Over the years, a number of other computational schemes have been
devised to provide a local analysis of the electronic structure in
molecules and solids. Here we briefly mention those most commonly used
in solid-state studies. The first choice is whether to work with the
electronic wavefunction or with the charge density.

One of the earliest and still most widely used wavefunction based
schemes is the ``Mulliken population analysis''
\cite{mulliken-jcp55}. This starts from a representation of the
density operator in an LCAO basis. If an extended basis, such as
planewaves, has been used, this can be obtained after first performing
a projection onto a suitable set of atomic orbitals
\cite{sanchez-portal-ssc95}.
%
Using the quantity $P_{A}$ introduced in \equ{PM} the Mulliken charge
on an atomic site A is given by
\beq
\label{eq:mull_charge}
Q_{A}=\sum_{i=1}^{\nwann} \langle\phi_{i}|P_{A}|\phi_{i}\rangle.
\eeq
The Mulliken scheme also provides a projection into local
angular-momentum eigenstates and an overlap (or bond) population between atom pairs.
%
The major disadvantage of the scheme is the fact that the absolute values
obtained have a marked dependence on the
LCAO basis. In fact, the results tend to become less meaningful as the
basis is expanded, as orbitals on one atomic site contribute to the
wavefunction on neighbouring atoms. However, it is generally accepted
that so long as calculations using the same set of local orbitals are
compared, trends in the values can provide some chemical
intuition \cite{segall-prb96}. An early application was to the study of bonding
at grain boundaries of TiO$_2$ \cite{dawson-prb96}.

An alternative approach is to work directly with the charge density. The
scheme of \textcite{hirshfeld-tca77} attempts to partition the charge density by
first defining a so-called prodensity for the system, typically a
superposition of free atom charge densities $\rho^i({\rr})$.
The ground-state charge density
is then partitioned between atoms according to the proportions of the
procharge at each point in space. This can easily be integrated to
give, for example, a total charge
%
\begin{equation}
Q_{\rm H}^i=\int {\bf dr} \rho({\rr}) \frac{\rho^i({\rr})}{\Sigma_i \rho^i({\rr})}
\end{equation}
%
for each atomic site.  Hirschfield charges have
recently been used to parameterize dispersion corrections to local
density functionals \cite{tkatchenko-prl09}.

Partitioning schemes generally make reference to some arbitrary
auxiliary system; in the case of Hirschfield charges, this is the
free-atom charge density, which must be obtained within some
approximation. In contrast, the ``atoms in molecules'' (AIM) approach developed
by \textcite{bader-cr91} provides a uniquely defined way
 to partition the charge density.  It uses the vector field corresponding to the
 gradient of the charge density.  In many cases the only maxima in the
 charge density occur at the atomic sites. As all field lines must
 terminate on one of these atomic ``attractors'', each point in space can
 be assigned to a particular atom according to the basin of attraction that is
 reached by following the
 density gradient. Atomic regions are now separated by zero-flux
 surfaces ${\rm S}({\rr}_s)$
 defined by the set of points $({\rr}_s)$ at which
%
\begin{equation}
\nabla\rho({\rr}_s)\cdot{\bf n}({\rr}_s)=0 ,
\end{equation}
%
where ${\bf n}({\rr}_s)$ is the unit normal to ${\rm S}({\rr}_s)$. Having
made such a division it is straightforward to obtain values for the
atomic charges (and also dipoles, quadrupoles, and higher moments).
The AIM scheme has been widely used
to analyse bonding in both molecular and solid-state systems, as well as to give
a localized description of response properties such as infra-red absorption
intensities \cite{matta-book-07}.

A rather different scheme is the ``electron localization function'' (ELF)
introduced by \textcite{becke-jcp90} as a simple measure of electron
localization in physical systems. Their original definition is based
on the same-spin pair probability density $P({\rr},{\rr}')$, i.e., the
probability to find two like-spin electrons at positions ${\rr}$ and
${\bf r}'$. \textcite{savin-ac92} introduced a form for the ELF
$\epsilon({\rr})$ which can be applied to an independent-particle
model:
%
\begin{equation}
\epsilon({\rr})=\frac{1}{1+(D/D_h)^2}\,,
\end{equation}
\begin{equation}
D=\frac{1}{2}\sum^{\nbw}_{i=1}|\nabla\psi_i|^2-\frac{1}{8}\frac{|\nabla\rho|^2}{\rho},
\end{equation}
\begin{equation}
D_h=\frac{3}{10}(3\pi^2)^{2/3}\rho^{5/3}, \rho=\sum^{\nbw}_{i=1}|\psi_i|^2,
\end{equation}
%
where
the sum is over all occupied orbitals.  $D$ represents the
difference between the non-interacting kinetic energy and the kinetic
energy of an ideal Bose gas. $D_h$ is the kinetic energy of a
homogeneous electron gas with a density equal to the local density. As
defined, $\epsilon({\rr})$ is a scalar function which ranges from 0 to
1. Regions of large kinetic energy (i.e., electron delocalization) have
ELF values close to zero while larger values correspond to paired
electrons in a shared covalent bond or in a lone pair. In a uniform
electron gas of any density, $\epsilon({\rr})$ will take the value of
$1/2$. Early application of the ELF in condensed phases provided
insight into the nature of the bonding at surfaces of
Al \cite{santis-ss00} and Al$_2$O$_3$ \cite{jarvis-jpc01},
and a large number of other applications have appeared since.
