% \section{WANNIER FUNCTIONS AS BASIS FUNCTIONS}

% \label{sec:basis}


In Ch.~\ref{sec:winterp}, we described the use of Wannier functions as
a compact tight-binding basis that represents a given set of energy
bands exactly, and that can be used to calculate a variety of
properties efficiently, and with very high accuracy. This is possible
because (a) WFs and bands are related by unitary transformations, and
(b) WFs are sufficiently localized that any resulting tight-binding
representation may be truncated with little loss of accuracy. In this
Chapter, we review two further general approaches to the use of WFs as
optimal and compact basis functions for electronic-structure
calculations that exploit their localization and transferability.

The first category includes methods in which WFs are used to go up in the
length scale of the simulations, using the results of electronic-structure 
calculations on small systems in order to construct accurate
models of larger, often meso-scale, systems. Examples include using
WFs to construct tight-binding Hamiltonians for large, structurally
complex nanostructures (in particular for studying quantum transport
properties), to parametrize semi-empirical force-fields, and to
improve the system-size scaling of quantum Monte Carlo (QMC)
and GW calculations, and the evaluation of exact-exchange integrals.

The second category includes methods in which WFs are used to
identify and focus on a desired, physically
relevant subspace of the electronic degrees of freedom that is singled
out (``downfolded'') for further detailed analysis or special treatment with a more
accurate level of electronic structure theory, an approach that is
particularly suited to the study of strongly-correlated electron
systems, such as materials containing transition metal, lanthanoid, or
even actinoid ions.


\subsection{WFs as a basis for large-scale calculations}
\label{sec:wfs-large-scale}

Some of the first works on linear-scaling electronic structure algorithms
\cite{yang-prl91,baroni-el92,galli-prl92,hierse-prb94} highlighted the connection 
between locality in electronic structure, which underpins linear-scaling algorithms
(Sec.~\ref{sec:linear-scaling}), and the
transferability of information across length-scales. In particular,
\textcite{hierse-prb94} 
considered explicitly two large systems A and B, different globally but which
have a certain similar local feature, such as a particular chemical
functionalization and its associated local environment, which we will
call C. They argued that the local electronic structure information
associated with C should be similar whether calculated from system A
or system B and, therefore,
that it should be possible to transfer this information from a
calculation on A in order to construct a very good
approximation to the electronic structure in the locality of feature C
in system B. In this way, large computational savings could be made on
the self-consistency cycle, enabling larger-scale calculations. 

The units of electronic structure information that
\textcite{hierse-prb94} used were localized non-orthogonal orbitals,
optimized in order to variationally minimize the total energy 
of the system. These orbitals are also referred to in the literature as ``support
functions''~\cite{hernandez-prb95} or ``non-orthogonal generalized
Wannier functions''~\cite{skylaris-prb02}. 

\subsubsection{Quantum transport}

A local representation of electronic structure is particularly
suited to the study of quantum transport, as we illustrate here in the case of
quasi-one-dimensional systems.
%
\begin{figure}
\begin{center}
\includegraphics[width=3.2in]{figs/lcr}
\caption{Top: Illustration of the lead-conductor-lead geometry
  used in quantum transport calculations. The conductor ($C$), left lead
  ($L$) and right lead ($R$) are described by Hamiltonians $H_C$, $H_L$
  and $H_R$, respectively. The coupling between adjacent regions is
  described by matrices $h_{LC}$ and $h_{CR}$. Bottom: The leads are
  split into principal layers corresponding to the sub-matrices in
  \equ{hamloc}. From \textcite{shelley-cpc11}.}
\label{fig:LCR}
\end{center}
\end{figure}
%
We consider a system composed of a conductor $C$ connected to two
semi-infinite leads, $R$ and $L$, as shown in Fig.~\ref{fig:LCR}. 
The conductance ${\cal G}$ through
a region of interacting electrons is related to the scattering
properties of the region itself via the Landauer
formula~\cite{landauer-pm70} 
%
\beq
{\cal G}(E) = \frac{2 e^2}{h} {\cal T}(E),
\label{eq:qc}
\eeq
%
where the transmission function ${\cal T}(E)$ is the probability that
an electron with energy $E$ injected at one end of the conductor will 
transmit to the other end. 
This transmission function can be expressed in terms of the Green's 
functions of the conductors and the coupling of the conductor to the
leads \cite{datta-95,fisher-prb91},
%
\beq
{\cal T}(E) = {\rm Tr}(\Gamma_{L} G_{C}^{r} \Gamma_{R} G_{C}^{a}),
\eeq
%
where $G_{C}^{\{r,a\}}$ are the retarded (r) and advanced (a) Green's
functions of the conductor, and $\Gamma_{\{L,R\}}$ are functions that
describe the coupling of the conductor to the left ($L$) and right ($R$)
leads. Since $G^{a}=(G^{r})^{\dagger}$, we consider $G^{r}$ only and drop
the superscript.

Expressing the Hamiltonian $H$ of the system in terms of a localized,
real-space basis set enables it to be partitioned without ambiguity
into sub-matrices that correspond to the individual subsystems. 
A concept that is particularly useful is that of a \emph{principal
  layer}~\cite{lee-prb81} (PL), which is a section of lead that is
sufficiently long such that $\matel{\chi_i^n}{H}{\chi_j^m}\simeq
0$ if $|m-n| \ge 2$, where $H$ is the Hamiltonian operator of
the entire system and $\ket{\chi_i^n}$ is the $i^{\rm th}$ basis
function in the $n^{\rm th}$ PL. Truncating the matrix elements of the
Hamiltonian in this way 
incurrs a small error which is systematically controlled by increasing
the size of the PL. The Hamiltonian matrix in this basis then takes
the block diagonal form (see also Fig.~\ref{fig:LCR})
%
\beq
H=\left( 
\begin{array}{ccccccc}
\ddots & \vdots & \vdots & \vdots & \vdots & \vdots & \iddots \\
\cdots & H_L^{\bar{0}\bar{0}} & H_L^{\bar{1}\bar{0}} & 0 & 0 & 0 &
\cdots \\
\cdots & H_L^{\bar{1}\bar{0}\dag}	& H_L^{\bar{0}\bar{0}}& h_{LC}
& 0 & 0 & \cdots \\
\cdots & 0 & h_{LC}^\dag & H_C & h_{CR} & 0 & \cdots \\
\cdots & 0 & 0 & h_{CR}^{\dag}	& H_R^{00} & H_R^{01} & \cdots \\
\cdots & 0 & 0 & 0 & H_R^{01\dag} & H_R^{00} & \cdots \\
\iddots & \vdots & \vdots & \vdots & \vdots &\vdots & \ddots \\
\end{array}\right),
\label{eq:hamloc}
\eeq
%
where $H_C$ represents the Hamiltonian matrix of the
conductor region, $H_L^{\bar{0}\bar{0}}$ and $H_R^{00}$ 
are those of each PL of the left and right leads, respectively,
$H_L^{\bar{1}\bar{0}}$ and $H_R^{01}$ are couplings between adjacent
PLs of lead, and $h_{LC}$
and $h_{CR}$ give the coupling between the conductor and the leads.

In order to compute the Green's function of the conductor one starts
from the equation satisfied by the Green's function $G$ of the whole
system, 
%
\beq
(\epsilon - H)G = \mathbb{I}
\label{eq:green}
\eeq
%
where $\mathbb{I}$ is the identity matrix, and $\epsilon = E + i\eta$,  where
$\eta$ is an arbitrarily small, real constant.

From \eqs{green}{hamloc}, it can be shown that
the Green's function of the conductor is then given by~\cite{datta-95}
%
\beq
G_{C} = (\epsilon -H_{C} -\Sigma_{L} -\Sigma_{R})^{-1} ,
\label{eq:gconduct}
\eeq
%
where we define $\Sigma_{L} = h_{LC}^\dagger g_{L} h_{LC}$ and $\Sigma_{R}
= h_{RC} g_{R} h_{RC}^\dagger$, the self-energy terms due to the
semi-infinite leads, and $g_{\{L,R\}}=(\epsilon-H_{\{L,R\}})^{-1}$,
the surface Green's functions of the leads. 

The self-energy terms can be viewed as effective Hamiltonians that
arise from the coupling of the conductor with the leads. Once the
Green's functions of the leads are known, the coupling functions
$\Gamma_{\{L,R\}}$ can be easily obtained as~\cite{datta-95} 
%
\beq
\Gamma_{\{L,R\}} = i [\Sigma_{\{L,R\}}^{r} - \Sigma_{\{L,R\}}^{a}] ,
\eeq
%
where $\Sigma_{\{L,R\}}^{a}=(\Sigma_{\{L,R\}}^{r})^{\dagger}$. 

As for the surface Green's functions $g_{\{L,R\}}$ of the
semi-infinite leads, it is well-known that any solid (or surface) can
be viewed as an infinite (semi-infinite in the case of surfaces) stack
of principal layers with nearest-neighbor
interactions~\cite{lee-prb81}. This corresponds to transforming the
original system into a linear chain of PLs. Within this
approach, the matrix elements of \equ{green} between layer orbitals
will yield a set of coupled equations for the Green's functions which
can be solved using an efficient iterative scheme due to
\textcite{lopezsancho-jpf84,lopezsancho-jpf85}. Knowledge of the
finite Hamiltonian sub-matrices in \equ{hamloc}, therefore, is
sufficient to calculate the conductance of the open
lead-conductor-lead system given by \equ{qc}.

%Knowledge of the bulk Green's function $G$ gives also direct
%information on the electronic spectrum via the spectral density of
%bulk electronic states $N(E) = -(1/\pi) \imag ({\rm Tr}G(E))$.

There are a number of possibilities for the choice of localized basis
$\ket{\chi}$. Early work used model tight-binding
Hamiltonians~\cite{chico-prb96, saito-prb96, anantram-prb98,
nardelli-prb99}, but the increasing 
sophistication of computational methods meant that more realistic first-principles
approaches could be adopted~\cite{nardelli-prb01,fattebert-book03}. 
Maximally-localized Wannier functions
were first used in this 
context by \textcite{calzolari-prb04}, who studied Al and C chains and
a (5,0) carbon nanotube with a single Si substitutional defect, and
by \textcite{lee-prl05}, who studied covalent functionalizations of metallic
nanotubes - capabilities now encoded in the open-source packages Wannier90 \cite{mostofi-cpc08} and
WanT \cite{want}. This
was quickly followed by a number of applications to ever more realstic
systems, studying transport through
molecular junctions~\cite{thygesen-cp05, strange-jcp08},
decorated carbon nanotubes and nanoribbons~\cite{lee-prl06, cantele-nl09, 
rasuli-prb10, li-acn11}, organic 
monolayers~\cite{bonferroni-n08}, 
and silicon nanowires~\cite{shelley-epl11}, as well as
more methodological work on improving the description of electron
correlations within this formalism~\cite{calzolari-n07,bonferroni-n08,
ferretti-prl05,ferretti-prb05}.

The formulation described above relies on a localized description of 
the electronic-structure problem, and it should be noted that several 
approaches to calculating electronic transport properties have been developed 
using localized basis sets rather than MLWFs, ranging from Gaussians~\cite{hod-jcp06} to 
numerical atomic orbitals~\cite{brandbyge-prb02, markussen-prb06, dasilva-prl08}.

The ideas of transferability of
local electronic structure information, discussed earlier in this section,
have been explored and developed in the
context of electronic-structure and transport calculations: the
suggestion that MLWFs can act as the building blocks for the 
Hamiltonians of large-scale nanostructures~\cite{lee-prl05}
has been applied to complex systems containing tens of thousands of
atoms~\cite{lee-prl06,cantele-nl09,shelley-cpc11,shelley-epl11,li-acn11}. 
In these approaches, the Hamiltonian matrices of
large systems for which a full,
conventional DFT calculation would be intractable, are built using
first-principles calculations performed on smaller fragments of the larger
system. The Hamiltonian matrix elements in the basis of MLWFs that are obtained from
the calculations on these fragments can be used to constract the entire 
Hamiltonian matrix of the desired system, 
and where the size of the fragments is a
systematically controllable determinant of the accuracy of the
method~\cite{shelley-cpc11}. 

%As noted in Sec.~\ref{sec:surface-bands},
%this ab initio tight-binding approach has been used recently in the 
%context of surface properties of topological insulators~\cite{zhang-njp10}.

An issue that arises when combining matrix elements from more than one 
DFT calculation into a single tight-binding Hamiltonian is that a common
reference potential must be defined.
For example, consider combining the results of a calculation on a 
perfect bulk material and one on the same material with an isolated 
structural defect. In the latter case, the diagonal (on-site) matrix
elements $\matel{w_n}{H}{w_n}$
%
in the system with the defect should converge to the same value as 
those in the pristine system as one goes further away from the location
of the defect.
%
With periodic boundary conditions, however, this is not guaranteed:
the difference between the matrix elements in the respective calculations
will, in general, tend to a non-zero constant, this being the difference
in reference potential between the two calculations.
%
In other words, the difference in value between the on-site matrix elements from 
the bulk calculation and 
those from the calculation with the defect, but far away from the location 
of the defect, give a direct measure of the potential mismatch between the 
two calculations. This offset can then be used to align all of the diagonal
matrix elements from one calculation with those of the other, prior to the 
them being combined in a tight-binding Hamiltonian for a larger nanostructure\footnote{
More generally, one could envisage an approach where electrostatic
self-consistency is attained solving the Poisson equation
for the entire nanostructure\cite{leonard-prl99}.}
%

The potential alignment approach described above has been used not only in 
the context of transport calculations, but also for the
determination of band offsets for calculating formation energies of point
defects~\cite{corsetti-prb11}, for developing tight-binding models of surfaces 
of topological insulators~\cite{zhang-njp10} (see 
Sec.~\ref{sec:surface-bands}), and for computing Schottky 
barriers~\cite{singhmiller-phd09}. 
%

\subsubsection{Semi-empirical potentials}

First-principles molecular dynamics simulations of large-scale
(thousands of atoms) systems for long (nanoseconds) time-scales are
computationally costly, if not intractable. Molecular dynamics
simulations with empirical interatomic potentials are a feasible
alternative and there is an ongoing effort in developing potentials
that are more accurate, more transferable and, therefore, more
predictive.
%
One approach in this direction is to fit the parameters that appear
within empirical potentials so that they reproduce target properties,
such as forces and stresses, obtained from accurate DFT
calculations on a large number of atomic 
configurations~\cite{ercolessi-el94}. 
%
In the particular class of ionic condensed matter systems, e.g.,
first and second row metal oxides, it is well-known that the
electronic properties of an ion can be significantly affected by its
coordination environment and, therefore, that it is also important to 
include an accurate description of polarization effects in any
interatomic potential. 
%
While simple potentials may attempt to account for these many-body
effects in an average manner, the result is that the potential loses
transferability and, hence, predictive power. 
%
As a result, there has been some effort in developing potentials such
that they are also fitted to information from DFT calculations 
regarding the electron distribution around each ion, in particular,
dipoles and quadrupoles. 

In this vein, \textcite{aguado-fd03} introduced the use of MLWFs in
order to calculate dipole and quadrupole moments in MgO and used these
to construct an interatomic potential based on the aspherical ion
method (AIM) potential~\cite{rowley-jcp98}. In an ionic system such as
this, MLWFs are found to be well-localized close to the ions and each can
therefore be associated unambiguously with a particular ion. 
The dipole moment ${\bf \mu}^{I}$ of each ion $I$ is
then calculated as 
%
\beq
\mu^{I}_{\alpha} = -2 \sum_{n\in I} \bar{r}_{n\alpha} + Z_I R^{I}_{\alpha},
\eeq
%
where $\alpha$ is a Cartesian component, $Z_I$ is the charge of ion
$I$ at position $\RR_I$, $\bar{r}_{n\alpha}$ is the centre of
the $n^{\mathrm{th}}$ MLWF and is given by
\equ{r} or \equ{rcenter}, and the factor of two accounts for spin degeneracy.

For the quadrupole moments $\theta^{I}_{\alpha\beta}$, the fact that
the MLWFs are localized within the simulation cell is exploited in
order to explicitly compute the real-space integral
%
\beq
\theta^{I}_{\alpha\beta} = -2\sum_{n\in I} \int_{V^{I}_{c}} 
\! d\rr \; \left|w_{n}(\rr)\right|^{2} \left[ 3
  r^{I}_{\alpha}r^{I}_{\beta} - (r^{I})^{2}\delta_{\alpha\beta}\right]/2,
\eeq
%
where $\rr^{I}=\rr - \RR^{I}$, and the integral is performed over a
sphere $V^{I}_{c}$, centered on $R^{I}$, such that the integral of the
electron density associated with the MLWF within this sphere is greater
than some threshold. The potential obtained
exhibits good transferability and the method has been used
to parametrize similar potentials for other alkaline earth
oxides~\cite{aguado-jcp03}, Al$_2$O$_3$~\cite{jahn-prb06}, and the
CaO-MgO-Al$_2$O$_3$-SiO$_2$ (CMAS) mineral
system~\cite{jahn-poteapi07}. 

The use of MLWFs for attempting to obtain better interatomic potentials
has not been limited to the 
solid state. In biomolecular simulations, an important factor
in developing accurate force-fields is having an accurate description
of the electrostatic potential. Starting from DFT calculations on
isolated molecules, \textcite{sagui-jcp04} 
partition the electronic charge density into contributions from
individual MLWFs and calculate the multipoles of each using an
order-by-order expansion of gauge-invariant
cumulants~\cite{resta-prl98, souza-prb00b} (the reader is referred to
\textcite{sagui-jcp04} for full details). Using fast particle mesh
Ewald and multigrid methods, these multipoles can then be 
used to generate the electrostatic potential. \textcite{sagui-jcp04} show that higher
order multipoles, e.g., up to hexadecapole, may be incorporated without
computational or numerical difficulty and that agreement with the
`exact' potential obtained from DFT is very good.

Finally, we note that recently
\textcite{rotenberg-prl10} have proposed force-fields whose
parametrization is based entirely on a partitioning of the
electronic density in terms of MLWFs. Their method has been applied
successfully to water and molten salts. It
remains to be seen whether the approach is extensible to more complex
or anisotropic systems. 
%
% AM: too much detail.
% In this so-called
% Wannier orbitals force-fields (WOFF) approach, the density due to each
% MLWF is replaced by an analytic frozen density that corresponds to a
% Slater orbital with the same center and spread as the MLWF. A further
% orientational averaging within a cut-off radius around each ion is
% performed, which defines a frozen electronic charge density associated with
% each ion. These frozen densities are then used to determine
% exchange-repulsion, electrostatic and dispersion
% interactions~\cite{silvestrelli-prl08}, while induction effects are
% included by calculating the response of the MLWFs to an external
% electric field. 


\subsubsection{Improving system-size scaling}

The localized nature of MLWFs in real-space makes them a natural and
appealing choice of basis for electronic structure calculations as the
sparsity exhibited by operators that are 
represented in a localized basis may
be exploited in order to achieve computational efficiencies and
improved scaling of numerical algorithms with respect to system
size. Recently, therefore, MLWFs have been used for this very purpose in a
number of contexts and we mention, in brief, some of them here.
% 

In quantum Monte Carlo (QMC) calculations~\cite{foulkes-rmp01}, a
significant computational 
effort is expended in evaluating the Slater determinant for a given
electronic configuration of $N$ electrons. These Slater determinants are usually
constructed from a set of extended single-particle states, obtained
from, e.g., a DFT or Hartree-Fock
calculation, represented in a basis set of, e.g., extended
plane-waves. This gives rise to $\mathcal{O}(N^3)$ scaling of the
computational cost of evaluating such a determinant.
%
\textcite{williamson-prl01} suggested, instead, to use
MLWFs that were smoothly truncated to zero beyond a certain cut-off
radius that is independent of system size. This ensures that each
electron falls only within the localization region of a fixed number
of MLWFs, thus reducing the asymptotic scaling by one factor of
$N$. Furthermore, by representing the MLWFs in a basis of localized spline
functions, rather than plane-waves or even Gaussian functions, the
evaluation of each orbital is rendered independent of system size,
thereby reducing the overall cost of computing the determinant of the
given configuration to $\mathcal{O}(N)$. 
%
More recently, rather than truncated MLWFs, the use of non-orthogonal
orbitals obtained by projection~\cite{reboredo-prb05} or other
localization criteria~\cite{alfe-jpcm04} has also been suggested.

In another development, \textcite{wu-prb09} use MLWFs in order to
compute efficiently Hartree-Fock exact-exchange integrals in extended
systems. Hybrid exchange-and-correlation
functionals~\cite{becke-jcp93} for DFT calculations, in which some
proportion of Hartree-Fock exchange is included in order to alleviate
the well-known problem of self-interaction that exists in local and semi-local
functionals such as the local-density approximation and its
generalized gradient-dependent variants, have been used relatively
little in extended systems.
%
This is predominantly due to the computational
cost associated with evaluating the exchange integrals between
extended eigenstates that are represented in a plane-wave basis
set. \textcite{wu-prb09} show that by performing a unitary
transformation of the eigenstates to a basis of MLWFs, and working in
real-space in order to exploit the fact that spatially distant MLWFs
have exponentially vanishing overlap, the number of such overlaps that
needs to be calculated scales linearly, in the limit of large system-size, 
with the number of orbitals (as opposed to
quadratically), which is a sufficient improvement to enable Car-Parrinello
molecular dynamics simulations with hybrid functionals.

Similar ideas that exploit the locality of
MLWFs have been applied  to many-body perturbation theory approaches
for going beyond DFT and Hartree-Fock calculations, for example, in the
contexts of the GW approximation~\cite{umari-prb09}, the evaluation
of local correlation in extended systems~\cite{buth-prb05,
  pisani-jcp05}, and the Bethe-Salpeter
equation~\cite{sasioglu-prb10}. The improved scaling and efficiency of
these approaches open the possibility of such calculations on larger
systems than previously accessible. 

Finally, we note that MLWFs have been used recently to compute van der
Waals interactions in an approximate but efficient
manner~\cite{silvestrelli-prl08, andrinopoulos-jcp11,
  silvestrelli-arxiv11}. 

% % AM: too much detail
% The GW approximation~\cite{hedin-pr65, rieger-cpc99} to many-body
% perturbation theory provides a means by which to obtain excitation
% spectra. Despite its success, such calculations come with a relatively high
% computational cost. A large portion of this cost is from determining
% density response functions, for which a basis of product functions between
% valence and conduction eigenstates of the Hamiltonian are required. In
% the case of using extended eigenstates, the number of product
% wavefunctions scales quadratically with the number of eigenstates (and
% hence the system size). \textcite{umari-prb09} show that using
% MLWFs instead, for example, introduces sparsity which, in the limit of 
% sufficiently large system sizes, reduces this scaling to linear, which
% opens the possibility of GW calculations on larger systems than currently
% accessible. 


\subsection{WFs as a basis for strongly-correlated systems}
\label{sec:wfs-min-models}

For many strongly-correlated electron problems, the essential physics
of the system can be explained by considering only a subset of the
electronic states. A recent example is understanding the behaviour
of unconventional (high-$T_{\rm c}$) superconductors, in which
a great deal of insight can be gained by considering only the MLWFs of
$p$ and $d$ character on Cu and O, respectively, for
cuprates~\cite{sakakibara-prl10}, and those on As and Fe,
respectively, for  the iron pnictides~\cite{kuroki-prl08,
  cao-prb08, hu-jpcc10, suzuki-jpsj11}. 
%
Other strongly-correlated materials for which MLWFs have been used to
construct minimal models to help understand the physics include
manganites~\cite{kovacik-prb10}, topological
insulators~\cite{zhang-np09, zhang-prb09} (see also
Sec.~\ref{sec:surface-bands}), and polyphenylene vinylene (PPV), in
particular relating to electron-hole excitations~\cite{karabunarliev-jcp03a,
  karabunarliev-jcp03b}.  

Below we outline some of the general principles behind the
construction and solution of such minimal models. 

\subsubsection{First-principles model Hamiltonians}
\label{sec:mod-ham}
Consider a strongly-correlated electron system described by a 
Hamiltonian of the form
%
\beq
\label{eq:scham1}
\hat{H} = \hat{H}_{0} + \hat{H}_{\rm int},
\eeq
%
where $\hat{H}_{0}$ contains the one-particle terms and $\hat{H}_{\rm
  int}$ the interaction (e.g., Coulomb repulsion) terms. 
%
In second-quantized notation and expressed in terms of a complete
tight-binding basis, this may be expressed as
%
\beq
\begin{split}
\hat{H} =& \sum_{ij}\sum_{\RR,\RR'} h_{ij}^{\RR\RR'}
            (\hat{c}_{i\RR}^{\dagger} \hat{c}_{j\RR'} + {\rm H.c.})
            \\& + \frac{1}{2}\sum_{ijkl} \sum_{\RR\RR'\RR''\RR'''}
            U_{ijkl}^{\RR\RR'\RR''\RR'''}
            \hat{c}_{i\RR}^{\dagger} \hat{c}_{j\RR'}^{\dagger}
            \hat{c}_{k\RR''}\hat{c}_{l\RR'''},  
\label{eq:scham2}
\end{split}
\eeq
%
where $\RR$ labels the correlated site, lower-case indices (such as
$i$ and $j$) refer to the orbital and spin degrees of freedom,
$\hat{c}_{i\RR}^{\dagger}$ ($\hat{c}_{i\RR}$) creates (annihilates) an
electron in the orbital $w_{i\RR}(\rr)$, and ${\bf h}$ and ${\bf U}$
contain the matrix elements of the single-particle and (screened)
interaction parts of the Hamiltonian, respectively.

Usually, a complete tight-binding representation of all the states of the
system is not required, and the essential physics can be described by
a smaller set of physically relevant orbitals. For example, those
states close to the Fermi level, or those of a particularly symmetry,
or those localized on specific sites, may be sufficient. In this way,
the size of the basis used to represent the many-body Hamiltonian is
greatly reduced, rendering \equ{scham2} more amenble to solution
(see, e.g., \textcite{solovyev-jpcm08}).

%
In this spirit, MLWFs obtained from DFT calculations have been used as
the orbital basis for these minimal models.\footnote{An alternative
  approach is to obtain localized orbitals via the downfolding method,
  discussed in Ch.~\ref{sec:locorb}.}
%
Advantages of using MLWFs include the fact that they can be
chosen to span precisely the required energy range (using the
disentanglement procedure outlined in Sec.~\ref{sec:disentanglement}),
and that they naturally incorporate hybridization 
and bonding appropriate to their local environment.

The single-particle hopping parameters of the model Hamiltonian are
obtained easily from the matrix elements of the DFT Hamiltonian
represented in the basis of MLWFs, using \equ{ham-R}.
%
The interaction parameters of the model Hamiltonian can be calculated,
for example, from either constrained DFT~\cite{dederichs-prl84,
  nakamura-prb06}, within the random phase
approximation~\cite{springer-prb98, aryasetiawan-prb04,
  solovyev-prb05}, or by direct calculation of the matrix elements of
a suitable screened Coulomb interaction between, for example,
MLWFs~\cite{miyake-prb06, nakamura-prb06}. It is interesting to note
that numerical evidence suggests that on-site Coulomb interactions
(both screened and bare), are maximized when calculated with a basis
of MLWFs~\cite{miyake-prb08} and, therefore, that MLWFs may be an
optimally localized basis for this purpose. This is perhaps not
surprising given the broad similarities between MLWFs and WFs obtained
via the Edmiston-Ruedenberg localization scheme~\cite{edmiston-rmp63},
discussed in Sec.~\ref{sec:alt-loc}, which maximizes the electronic 
Coulomb self-interaction of each orbital.

%
Once the parameters of the model have been determined, the model
Hamiltonian is then said to be ``from first-principles'', in the sense
that the parameters of the model are determined from DFT rather than
by fitting to experiments. 
%
The many-body Hamiltonian in the minimal basis of MLWFs may be then
solved by direct diagonalization, or one of a number of other
approaches that are too numerous to review here but which include, for
example, the dynamical mean-field theory (DMFT) approach. 
DMFT maps
the many-body problem on to an Anderson impurity 
model~\cite{anderson-pr61} in which on-site correlation is treated
non-perturbatively and the correlated sites are coupled to a
self-consistent energy bath that represents the rest of the system.
%
The impurity sites, also known as the ``correlated subspaces'', are
defined by localized projector functions and MLWFs are a common
choice~\cite{lechermann-prb06, trimarchi-jpcm08, weber-prb10}. In
particular, one would typically choose orbitals of $d$ or $f$
character associated with transition metal, lanthanoid or actinoid ions.
%
The Green's function for the impurity site is calculated
self-consistently, for example, by a numerical functional integration
(which constitutes the bulk of the computation). Further
self-consistency with the DFT ground-state may also be attained by
using the solution to the impurity problem to update the 
electronic density that is then used to construct an updated Kohn-Sham
potential, which determines a new set of eigenstates, MLWFs and,
hence, model Hamiltonian parameters that can then be fed back in to
the DMFT cycle. The reader is referred to \textcite{kotliar-romp06}
and \textcite{held-aip07} for further details; other examples that use
localized Wannier functions or generalized tight-binding models to
address correlated electrons problems can be found in  
\textcite{amadon-prb08, anisimov-prb05a, held-jpcm08, korotin-epjb08,
  korshunov-prb05, ku-prl02}.

\subsubsection{The DFT + Hubbard $U$ method}

In the model Hamiltonian approaches above, the results of a DFT
calculation are used to parameterize the model Hamiltonian of a
strongly-correlated electron system. In constrast, in the DFT +
$U$~\cite{anisimov-prb91, anisimov-prb93} method, 
the energy functional in DFT is explicitly augmented with a term that
aims to account for the inadequacies of local and semi-local exchange and
correlation functionals for describing localized electronic states,
such as those associated with $d$ and $f$ electrons,
which are most likely to suffer from spurious self-interaction effects
and the absence of appropriate on-site Coulomb repulsion. 
%

In DFT + $U$, these sites are defined by a set of
``Hubbard projectors'', which are localized, atomic-like
orbitals of $d$ or $f$ character, centered on the ions
\footnote{The strength of the 
  correction $U$, while usually treated as an adjustable parameter
  that is chosen
  in order to reproduce experimental observables such as band-gaps,
  may also be calculated from first-principles using,
  e.g., constrained DFT~\cite{dederichs-prl84} or
  linear-response theory~\cite{cococcioni-prb05,kulik-prl06}, thus
  removing any arbitariness in its value.}.
%
MLWFs are an appealing choice of Hubbard projectors,
since they are able to capture hybridization effects of the correlated
sites with their local
environment, while retaining the desirable property of
localization that is the hallmark of orbitals affected by large
self-interaction errors. Recently, they have been used successfully for DFT + $U$
calculations by a  number of authors~\cite{fabris-prb05, miyake-prb08,
  anisimov-jpcm07}.
Last, it has been noted that the precise choice of projectors,
which is to some extent arbitrary, can affect calculated observables
and, therefore, constitutes an additional adjustable
parameter in the calculation~\cite{oregan-prb10}.

%{\bf Also mention self-interaction correction
%  (SIC)~\cite{perdew-prb81} and MLWFs~\cite{stengel-prb08}?}

% AAM: TOO MUCH DETAIL.
% The correction
% to the DFT energy then takes the form\footnote{\equ{E_U} is the
%   simplified, rotationally-invariant form of the Hubbard correction,
%   suggested by \textcite{dudarev-prb98}, and includes a term that
%   approximately subtracts double-counting of Coulomb interactions
%   already included in the DFT exchange and correlation
%   functional. For more general rotationally-invariant forms of $E_U$,
%   the reader is referred to \textcite{liechtenstein-prb95,
%     anisimov-jpcm97b}.} 
% %
% \beq
% E_{U} = \frac{1}{2}\sum_{I\sigma} U^{I\sigma} 
% {\rm tr} \left[ {\bf n}^{I\sigma} (\mathbb{I} - {\bf{n}}^{I\sigma}) \right],
% \label{eq:E_U}
% \eeq
% %
% where $U^{I\sigma}$ is the strength of the screened interaction on
% site $I$ for electrons of spin $\sigma$, and ${\bf n}^{I\sigma}$ is the
% occupancy matrix of site $I$ for electrons of spin $\sigma$, given by
% %
% \beq
% n_{mm'}^{I\sigma} = \matel{\chi_{m}^{I}}{\hat{\rho}^{\sigma}}{\chi_{m'}^{I}},
% \eeq
% %
% where $\hat{\rho}^{\sigma}$ is the single-particle density matrix for
% electrons of spin $\sigma$. 

% \subsubsection{Van der Waals interactions}

% \dvm{I tentatively moved this here from the summary chapter.
% Can it be made to work here?}

% In a different context, a manifestation of many-body correlations
% is the presence of van der Waals (also known as dispersion)
% interactions.  [Discuss here
% Silvestrelli's work on van der Waals interactions
% \cite{
% silvestrelli-prl08,
% silvestrelli-cpl09,
% silvestrelli-jpca09,
% silvestrelli-jcp09}
% ]
% \amm{Now inline:}
% Silvestrelli's work on vdW doesn't fit particularly well elsewhere 
% and is very recent, so perhaps worth a mention in the summary
% chapter; my student has a paper in review with J Chem Phys at the
% moment on the problems associated with Silvestrelli's method and how to cure
% them [Andrinopoulos, Hine and Mostofi, submitted to J Chem Phys
% (2011)].






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%% Original text from Nicola from here onwards %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% \subsection{First-principles tight-binding}
% \label{s:trans}
% \subsubsection{Electronic structure}

% \subsubsection{Landauer conductance}
% We consider a system composed of a conductor $C$ connected to two
% semi-infinite leads, $R$ and $L$.  
% The conductance through
% a region of interacting electrons is
% related to the scattering properties of the region itself via the
% Landauer formula \cite{landauer-pm70}:
% \begin{equation}
% {\cal C} = {2 e^2 \over h} {\cal T},
% \end{equation}
% where ${\cal T}$ is the transmission function and ${\cal C}$ is the
% conductance. ${\cal T}$ is the probability that an electron
% injected at one end of the conductor will transmit to the other
% end. This transmission function can be expressed in terms of the Green's
% functions of the conductors and the coupling of the conductor to the
% leads \cite{datta-95,fisher-prb91}:
% \begin{equation}
% {\cal T} = {\rm Tr}(\Gamma_L G_C^r \Gamma_R G_C^a),
% \end{equation}
% where $G_C^{\{r,a\}}$ are the retarded and advanced Green's functions of
% the conductor, and $\Gamma_{\{L,R\}}$ are functions that describe the
% coupling of the conductor to the leads. To compute the Green's
% function of the conductor we start from the equation for the Green's
% function of the whole system:
% \begin{equation}
% (\epsilon - H)G = I
% \label{green}
% \end{equation}
% where $\epsilon = E+{\rm i}\eta$ with $\eta$ arbitrarily small
% and $I$ is the identity matrix. In the hypothesis that the
% Hamiltonian of the system can be expressed in a discrete real-space
% matrix representation, the previous equation corresponds to the
% inversion of an infinite matrix for the open system, consisting of the
% conductor and the semi-infinite leads. The above Green's function can
% be partitioned into sub-matrices that correspond to the individual
% subsystems:
% %yes
% \begin{equation}
% \left(
% \begin{array}{lll}
% G_L  & G_{LC} &  G_{LCR}
% \\
% G_{CL} & G_C & G_{CR} 
% \\
% G_{LRC} & G_{RC} & G_R
% \end{array}
% \right)
% =
% \end{equation}
% $$
% \left(
% \begin{array}{ccc}
% (\epsilon -H_L) & h_{LC} & 0 
% \\
% h_{LC}^\dagger & (\epsilon -H_C) & h_{CR}
% \\
% 0 & h_{CR}^\dagger & (\epsilon -H_R)
% \end{array}
% \right)^{-1},
% \label{condmatr}
% $$
% where the matrix $(\epsilon -H_C)$ represents the finite isolated
% conductor, $(\epsilon -H_{\{R,L\}})$ represent the infinite leads, and
% $h_{CR}$ and $h_{LC}$ are the coupling matrices that will be non-zero
% only for adjacent points in the conductor and the leads,
% respectively. From this equation it is straightforward to obtain an
% explicit expression for $G_C$ \cite{datta-95}:
% %yes
% \begin{equation}
% G_C = (\epsilon -H_C -\Sigma_L -\Sigma_R)^{-1}
% \label{gconduct}
% \end{equation}
% where we define $\Sigma_{L} = h_{LC}^\dagger g_L h_{LC}$ and $\Sigma_R
% = h_{RC} g_R h_{RC}^\dagger$ as the self-energy terms due to the
% semi-infinite leads and $g_{\{L,R\}}=(\epsilon -H_{\{L,R\}})^{-1}$ are
% the leads' Green's functions. The self-energy terms can be viewed as
% effective Hamiltonians that arise from the coupling of the conductor
% with the leads. Once the Green's functions are known, the coupling
% functions $\Gamma_{\{L,R\}}$ can be easily obtained as \cite{datta-95}:
% \begin{equation}
% \Gamma_{\{L,R\}} = {\rm i}[\Sigma_{\{L,R\}}^r - \Sigma_{\{L,R\}}^a]
% \end{equation}
% where the advanced self-energy $\Sigma_{\{L,R\}}^a$ is the Hermitian
% conjugate of the retarded self-energy $\Sigma_{\{L,R\}}^r$. The core
% of the problem lies in the calculation of the Green's functions of the
% semi-infinite leads. In what follows we will present an efficient
% approach to compute the self-energy terms in the general case of an
% arbitrary localized-orbital Hamiltonian.

% \subsubsection{\it Transmission through a bulk system}
% \label{ss:bulk}

% It is well known that any solid (or surface) can be viewed as an
% infinite (semi-infinite in the case of surfaces) stack of principal
% layers with nearest-neighbor interactions \cite{lee-prb81}. This
% corresponds to transforming the original system into a linear chain of
% principal layers. Within this approach, the matrix elements of
% Eq. (\ref{green}) between layer orbitals will yield a set of equations
% for the Green's functions:
% %yes
% \begin{equation}
% \begin{array}{ccl}
% (\epsilon -H_{00}) G_{00} & = & I + H_{01}G_{10}\\
% (\epsilon -H_{00}) G_{10} & = & H_{01}^\dagger G_{00} + H_{01}G_{20}\\
% \multicolumn{3}{c}{\dotfill}\\
% (\epsilon -H_{00}) G_{n0} & = & H_{01}^\dagger G_{n-1,0} + H_{01}G_{n+1,0}
% \end{array}
% \label{serie}
% \end{equation}
% where $H_{nm}$ and $G_{nm}$ are the matrix elements of the Hamiltonian
% and the Green's function between the layer orbitals, and we assume that in
% a bulk system $H_{00} =H_{11}=\dots$ and
% $H_{01}=H_{12}=\dots$. Following Lopez-Sancho {\it et
% al.} \cite{lopezsancho-jpf84,lopezsancho-jpf85}, this chain can be transformed in order to
% express the Green's function of an individual layer in terms of the
% Green's function of the preceding (or following) one. This is done
% via the introduction of the transfer matrices $T$ and $\overline{T}$,
% defined such that $G_{10}=TG_{00}$ and $G_{00}=\overline{T}G_{10}$.
% The transfer matrix can be easily computed from the Hamiltonian matrix
% elements via an iterative procedure, as outlined in
% Ref. \cite{lopezsancho-jpf84,lopezsancho-jpf85}. 
% In particular $T$ and $\overline T$ can be written as
% \[
% \begin{array}{c}
% T = t_0 + \tilde{t}_0t_1 + \tilde{t}_0\tilde{t}_1t_2+\ldots+
% \tilde{t}_0\tilde{t}_1\tilde{t}_2\cdots t_n \\
% \overline T = \tilde{t}_0 + t_0\tilde{t}_1 +t_0t_1\tilde{t}_2 +\ldots+
% t_0t_1t_2\cdots\tilde{t}_n
% \end{array}
% \]
% where $t_i$ and $\tilde{t}_i$ are defined via the recursion formulas:
% \[
% \begin{array}{l}
% t_i = (I -t_{i-1}\tilde{t}_{i-1} - \tilde{t}_{i-1}t_{i-1})^{-1} 
% t_{i-1}^2, \\
% \tilde{t}_i = (I -t_{i-1}\tilde{t}_{i-1} - \tilde{t}_{i-1}t_{i-1})^{-1} 
% \tilde{t}_{i-1}^2
% \end{array}
% \]
% \noindent and
% \[
% \begin{array}{l}
% t_0 = (\epsilon - H_{00})^{-1} H_{01}^\dagger,\\
% \tilde{t}_0 = (\epsilon - H_{00})^{-1} H_{01}.
% \end{array}
% \]
% The process is repeated until $t_n,\tilde{t}_n \le \delta$ with
% $\delta$ arbitrarily small.
% With this proviso, we can write the bulk Green's function as:
% \begin{equation}
% G(E) = (\epsilon - H_{00} - H_{01}T -H_{01}^\dagger\overline{T})^{-1}.
% \end{equation}
% If we compare the previous expression with Eq. (\ref{condmatr}), in the
% hypothesis of leads and conductors being of the same material (bulk
% conductivity), we can identify the present bulk system,
% or rather one of its principal layers, with the conductor $C$, so that
% $H_{00}\equiv H_C$, $H_{01}\equiv h_{CR}$ and $H_{01}^\dagger\equiv
% h_{LC}^\dagger$.  In particular, by comparing with Eq. \ref{gconduct},
% we obtain the expression of the self-energies of the conductor-leads
% system:
% \begin{equation}
% \Sigma_L = H_{01}^\dagger \overline T, ~~~~~~~~~
% \Sigma_R = H_{01} T.
% \end{equation}
% The coupling functions are then obtained from the sole knowledge of the
% transfer matrices and the coupling Hamiltonian matrix elements:
% $\Gamma_L=-{\rm Im}(H_{01}^\dagger \overline T)$ and $\Gamma_R=-{\rm
% Im}(H_{01} \overline T)$. 
% The knowledge of the bulk Green's
% function $G$ gives also direct informations on the electronic spectrum
% via the spectral density of bulk electronic states:
% $
% N(E) = -(1/\pi) {\rm Im}({\rm Tr}G(E)).
% $
% \subsection{Model Hamiltonians}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% End of original text from Nicola %%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

