% \section{WANNIER INTERPOLATION}
% \label{sec:winterp}

%      Band energies, group velocities, effective mass
%      Other quantities [Berry curvature, etc.]
%      El-phonon interactions
%      Wannier interpolation for GW

%      Authors: Ivo 1st, Jonathan 2nd.


%\subsection{Overview}

\begin{figure*}
\begin{center}
\includegraphics[width=13.5cm]{figs/winterp}
\caption{(Color online) Schematic overview of the Wannier
  interpolation procedure.  The left panel shows the 
  BZ mesh $\qq$ used in the first-principles calculation, where the quantity of interest $f(\qq)$ is explicitly
  calculated. The Wannier-transformed quantity $F(\RR)$ is strongly
  localized near the origin of the equivalent supercell, shown in the
  middle panel with WFs at the lattice points. The right panel shows a
  dense mesh of interpolation points $\kk$ in the BZ, where the
  quantity $f(\kk)$ is evaluated at low cost starting from $F(\RR)$.}
\label{fig:interp_schem}
\end{center}
\end{figure*}

Localized Wannier functions are often introduced in textbooks as a
formally exact localized basis spanning a band, or a group of bands,
and their existence provides a rigorous justification for the
tight-binding (TB) interpolation
method~\cite{ashcroft-book76,harrison-book80}.

In this section we explore the ways in which WFs can be used as an
exact or very accurate TB basis, allowing to perform, very efficiently and accurately,
a number of operations on top of a conventional first-principles
calculation. The applications of this ``Wannier interpolation''
technique range from simple band-structure plots to the evaluation of
various physical quantities as BZ integrals.  The method is
particularly useful in situations where a very fine sampling of the BZ
is required to converge the quantity of interest. This is often the
case for metals, as the presence of a Fermi surface introduces sharp
discontinuities in $\kk$-space.

The Wannier interpolation procedure is depicted schematically in
Fig.~\ref{fig:interp_schem}.  The actual first-principles calculation
is carried out on a relatively coarse uniform reciprocal-space mesh
$\qq$ (left panel), where the quantity of interest $f(\qq)$ is
calculated from the Bloch eigenstates. The states in the selected
bands are then transformed into WFs, and $f(\qq)$ is transformed
accordingly into $F(\RR)$ in the Wannier representation (middle
panel). By virtue of the spatial localization of the WFs, $F(\RR)$
decays rapidly with $|\RR|$. Starting from this short-range real-space
representation, the quantity $f$ can now be accurately interpolated
onto an arbitrary point $\kk$ in reciprocal space by carrying out an
inverse transformation (right panel).  This procedure will succeed in
capturing variations in $f(\kk)$ over reciprocal lengths smaller than
the first-principles mesh spacing $\Delta q$, provided that the linear
dimensions $L=2\pi/\Delta q$ of the equivalent supercell are large
compared to the decay length of the WFs.


%%%OLD%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Figure~\ref{fig:interp_schem} depicts schematically the calculation of
%BZ integrals using Wannier interpolation. The actual first-principles
%calculation is carried out on the relatively coarse uniform $\kk$-point
%mesh on the left-hand-side, and the eigenstates in the bands of
%interest are transformed into WFs. The only requirement on the mesh
%spacing $\Delta k$ is that the linear dimension $L=2\pi/\Delta k$ of
%the equivalent supercell should be large compared to the decay length
%of the WFs.  For many applications this first-principles mesh will not
%be sufficiently dense to perform the BZ integration accurately, as the
%integrand may vary significantly over reciprocal lengths smaller than
%$\Delta k$.  Such fine-scale features can be recovered in two steps.

%First, the Wannier matrix elements of the operator(s) associated with
%the quantity of interest are evaluated. Because of the spatial
%localization of the WFs, these matrix elements decay rapidly with the
%separation between the WFs, so that the extra information that
%would be gained by
%further refining the {\it ab initio} $\kk$-mesh is negligible.

%From this short-range real-space representation, the reciprocal-space
%quantity of interest can now be accurately interpolated onto any
%desired $\kk$-point.  The computational cost of the interpolation step
%is very low, allowing to perform the numerical integration on
%ultradense grids. Unlike the {\it ab initio} grid, the interpolation
%grid does not have to be uniform, and if the integrand displays strong
%and rapid variations in certain regions of $\kk$-space, it can be
%refined adaptively in order to accelerate convergence.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\subsection{Band-structure interpolation}
\label{sec:band-interpolation}

The simplest application of Wannier interpolation is to generate
band-structure plots. We shall describe the procedure in some detail, as
the same concepts and notations will reappear in the more advanced
applications to follow. 

From the WFs spanning a group of $\nwann$ bands, a set of Bloch-like states
can be constructed using \equ{wannier-inverse}, which we repeat here
with a slightly different notation, 
%
\beq
\label{eq:bloch-w-gauge}
\ket{\psi_{n\kk}^\ww}=\sum_\RR\,e^{i\kk\cdot\RR}\ket{\RR n}
\hspace{1cm} (n=1,\ldots, \nwann).
\eeq
% 
This has the same form as the Bloch-sum formula in tight-binding
theory, with the WFs playing the role of the atomic orbitals.  The
superscript W serves as a reminder that the states
$\ket{\psi_{n\kk}^\ww}$ are generally not eigenstates of the
Hamiltonian.\footnote{In Ch.~\ref{sec:theory} the rotated Bloch states
  $\ket{\psi_{n\kk}^\ww}$ were denoted by $\ket{\tilde{\psi}_{n\kk}}$,
  see \equ{blochlike}.}  We shall say that they belong to the {\it
  Wannier gauge}.

At a given $\kk$, the Hamiltonian matrix elements in the space of the
$\nwann$ bands is represented in the Wannier gauge by the matrix
%
\beq
\label{eq:h_kw}
H_{\kk,nm}^\ww
=\bra{\psi_{\kk n}^\ww}H\ket{\psi_{\kk m}^\ww}
=\sum_\RR\,e^{i\kk\cdot\RR}\bra{\0 n}H\ket{\RR m}.
%
\eeq
%
In general this is a non-diagonal matrix in the band-like indices,
and the interpolated eigenenergies are obtained by diagonalization,
%
\beq
\label{eq:h_k-diag}
H^\hh_{\kk,nm}
=\left[ U^\dagger_\kk H^\ww_\kk U_\kk\right]_{nm}=
\delta_{nm}\enkh.
\eeq
%

In the following, it will be useful to view the unitary matrices
$U_\kk$ as transforming between the Wannier gauge on the one hand,
and the Hamiltonian (H) gauge (in which the projected Hamiltonian
is diagonal) on the other.\footnote{The
  unitary matrices $U_\kk$ are related to, but not the same as, the
  matrices $U^{(\kk)}$ introduced in \equ{blochlike}. The latter are
  obtained as described in Secs.~\ref{sec:projection} and
  \ref{sec:maxloc}. In the present terminology, they transform from the
  Hamiltonian to the Wannier gauge on the 
  mesh used in the first-principles calculation. Instead, $U_\kk$ transforms from the Wannier to the Hamiltonian
  gauge on the interpolation mesh. That is, the matrix $U_\kk$
  is essentially an interpolation of the
  matrix $\left[ U^{(\kk)}\right]^\dagger$.}
%
From this point forward we adopt a condensed notation in
which band indices are no longer written explicitly, so that,
for example,
$H_{\kk,nm}^\hh=\bra{\psi^\hh_{\kk n}}H\ket{\psi^\hh_{\kk m}}$
is now written as $H_\kk^\hh=\bra{\psi^\hh_\kk}H\ket{\psi^\hh_\kk}$,
and matrix multiplications are implicit.  Then \equ{h_k-diag}
implies that the transformation law for the Bloch states is
%
\beq
\label{eq:psik-h}
\ket{\psi_\kk^\hh}=\ket{\psi_\kk^\ww}U_{\kk}.
\eeq
%

If we insert into Eqs.~(\ref{eq:h_kw}) and (\ref{eq:h_k-diag}) a
wavevector belonging to the first-principles grid, we simply recover
the first-principles eigenvalues $\enk$, while for arbitrary $\kk$ the
resulting $\enkh$ interpolate smoothly between the values on the
grid. (This is strictly true only for an isolated group of bands. When
using disentanglement, the interpolated bands can deviate from the
first-principles ones outside the inner energy window, as discussed in
Sec.~\ref{sec:entangled} in connection with Fig.~\ref{fig:si_proj}.)

Once the matrices $\bra{\0}H\ket{\RR}$ have been tabulated, the
band structure can be calculated very inexpensively by Fourier transforming
[\equ{h_kw}] and diagonalizing [\equ{h_k-diag}] matrices of rank
$\nwann$. Note that $\nwann$, the number of WFs per cell, is typically
much smaller than the number of basis functions (e.g., plane waves)
used in the first-principles calculation.

In practice the required matrix elements are obtained by inverting
\equ{h_kw} over the first-principles grid,\footnote{In this section we
  adopt the normalization convention of \equ{discrete-abnormal}.
}
%
\beq
\label{eq:ham-R}
\begin{split}
\bra{\0}H\ket{\RR}&=\frac{1}{N}\sum_\qq\,e^{-i\qq\cdot\RR}
\bra{\psi_{\qq}^\ww}H\ket{\psi_{\qq}^\ww}\\
&=\frac{1}{N}\sum_\qq\,e^{-i\qq\cdot\RR}
        V_\qq^\dagger E_\qq V_\qq.
\end{split}
\eeq
%
Here $N$ is the number of grid points, $E_\qq$ is the diagonal matrix
of first-principles eigenenergies, and $V_\qq$ is the matrix defined in
\equ{wannierization-matrix}, which converts the $\nbands_\qq$ {\it ab
  initio} eigenstates into the $\nwann\leq \nbands_\qq$ Wannier-gauge
Bloch states,
%
\beq
\label{eq:psik-w}
\ket{\psi_\qq^\ww}=
\ket{\psi_\qq}V_{\qq}.
\eeq


%Recall that the outcome of
%the Wannierization
%
%
%procedure is a set of rectangular matrices ${\cal U}_\qq$ which
%convert the $\nbands_\qq$ first-principles eigenstates at each $\qq$
%into $\nwann\leq \nbands_\qq$ Wannier-gauge Bloch states. In matrix
%notation we have
%%
%\beq
%\label{eq:psik-w}
%\ket{\psi_\qq^\ww}=
%\ket{\psi_\qq}{\cal U}_{\qq}
%\eeq
% 
%so that finally
%%
%\beq
%\label{eq:ham-R}
%\bra{\0}H\ket{\RR}=\frac{1}{N}\sum_\qq\,e^{-i\qq\cdot\RR}
%        {\cal U}_\qq^\dagger E_\qq{\cal U}_\qq,
%\eeq
%
%where $E_\qq$ denotes the diagonal matrix containing the {\it ab
%  initio} eigenenergies $\epsilon_{n\qq}$.

The strategy outlined above~\cite{souza-prb01} is essentially the
Slater-Koster interpolation method.  However, while the
Hamiltonian matrix elements in the localized basis are treated as
adjustable parameters in empirical TB, they are calculated
from first-principles here.  A similar interpolation strategy is widely
used to obtain phonon dispersions starting from the interatomic force
constants calculated with density-functional perturbation
theory~\cite{baroni-rmp01}.  We shall return to this analogy between
phonons and tight-binding electrons \cite{martin-book04} when
describing the interpolation of the electron-phonon matrix elements in
Sec.~\ref{sec:el-ph}.

Wannier band-structure interpolation is extremely accurate. By virtue
of the exponential localization of the WFs within the periodic
supercell (see Footnote \ref{foot:discretization}), the magnitude of the
matrix elements $\bra{\0}H\ket{\RR}$ decreases rapidly with
$|\RR|$, and this exponential localization is preserved even in the
case of metals, provided a smooth subspace has been disentangled.
As the number of lattice vectors included in the summation in
\equ{h_kw} equals the number of first-principles mesh points, beyond a
certain mesh density the error incurred decreases
exponentially~\cite{yates-prb07}.
%the density of the first-principles grid~\cite{yates-prb07},
%
%\ism{Need to say that the magnitude of $\bra{\0}H\ket{\RR}$
%decreases rapidly with $|\RR|$. Nicola asked to mention that the
%criterion for ``optimal interpolators'' is that this decrease is as fast as
%possible (this is not the same as the maxloc requirement we use to
%construct WFs.)}
%
%so that a relatively coarse grid can be used in practice.  
In the following we will illustrate the method with a few
representative examples selected from the literature.

\subsubsection{Spin-orbit-coupled bands of bcc Fe}
\label{sec:fe-bands}

As a first application, we consider the relativistic band structure of
bcc Fe.  Because of the spin-orbit interaction, the spin density is
not perfectly collinear, and the Bloch eigenstates are spinors.  As
mentioned in Sec.~\ref{sec:mag}, spinor WFs can be constructed via a
trivial extension of the procedure described in Sec.~\ref{sec:theory}
for the non-magnetic (spinless) case.  It is also possible to further
modify the wannierization procedure so as to produce two separate
subsets of spinor WFs: one with a predominantly spin-up character, and
the other with a predominantly spin-down character~\cite{wang-prb06}.

\begin{figure}
\begin{center}
\includegraphics[width=3.2in]{figs/Fe-bands}
\caption{Band structure of bcc Fe with spin-orbit coupling
  included. Solid lines: original band structure from a conventional
  first-principles calculation.  Dotted lines: Wannier-interpolated
  band structure. The zero of energy is the Fermi level. From
  \textcite{wang-prb06}.}
\label{fig:Fe-bands}
\end{center}
\end{figure}

Using this modified procedure, a set of nine disentangled WFs per spin
channel was obtained for bcc Fe by \textcite{wang-prb06}, consisting
of three $t_{2g}$ $d$-like atom-centered WFs and six $sp^3d^2$-like
hybrids pointing along the cubic directions.  An inner energy window
was chosen as indicated in Fig.~\ref{fig:Fe-bands}, so that these 18
WFs describe exactly all the the occupied valence states, as well as
the empty states up to $E_{\rm win}$, which was set at approximately
18~eV above the Fermi level.

\begin{figure}[b]
\begin{center}
\includegraphics[width=2.6in]{figs/Fe-avoided-crossing}
\caption{(Color online)  Wannier-interpolated relativistic band
  structure of ferromagnetic bcc Fe along $\Gamma$--H. The bands are
  color-coded according to the expectation value of $S_z$: red for
  spin up and blue for spin down. The energies in eV are referred to
  the Fermi level. The vertical dashed lines indicate $\kk$-points on
  the mesh used in the first-principles calculation for constructing the WFs. From
  \textcite{yates-prb07}.}
\label{fig:Fe-avoided-crossing}
\end{center}
\end{figure}

The interpolated bands obtained using an $8\times 8\times 8$
$\qq$-grid in the full BZ are shown as dashed lines in
Fig.~\ref{fig:Fe-bands}.  The comparison with the first-principles
bands (solid lines), reveals essentially perfect agreement within the
inner window.  This is even more evident in
Fig.~\ref{fig:Fe-avoided-crossing}, where we zoom in on the
interpolated band structure near the Fermi level along $\Gamma$--H, and
color-code it according to the spin-projection along the quantization
axis.  The vertical dotted lines indicate points on the $\qq$-mesh.
For comparison, we show as open circles the eigenvalues calculated
directly from first-principles around a weak spin-orbit-induced
avoided crossing between two bands of opposite spin. It is apparent
that the interpolation procedure succeeds in resolving details of the
true band structure on a scale much smaller than the spacing between
$\qq$-points.

%{\scr \textcolor{red}{ The ability to follow the exact {\it ab
%      initio} bands on very small length scales in $\kk$-space is made
%    possible by the fact that the information about the connectivity
%    of the bands is encoded in the overlaps between Bloch eigenstates
%    at neighboring $\qq$-points [Eq.~(\ref{eq:Mkb})], the basic input
%    to the Wannierization procedure. Instead, traditional
%    interpolation schemes, which rely solely on the knowledge of the
%    energy eigenvalues on a grid, can easily mistake band crossings
%    for avoided crossings, or vice versa. This feature of Wannier
%    interpolation allows for a very stable and accurate determination
%    of band derivatives, as will be discussed in
%    Sec.~\ref{sec:band-der}.  } } 
%%
%\ism{\textcolor{blue}{I am tempted to remove the paragraph in red.}}


\subsubsection{Band structure of a metallic carbon nanotube}
\label{sec:nanotube-bands}

As a second example, we consider Wannier interpolation in large
systems (such as nanostructures), that are often sampled only at the
zone center. We consider here a (5,5) carbon nanotube, studied in a
100-atom supercell (i.e.  five times the primitive unit cell) and with
$\Gamma$-point sampling.  The system is metallic, and the
disentanglement procedure is used to generate well-localized WFs,
resulting in either bond-centered combinations of $sp^2$ atomic
orbitals, or atom-centered $p_z$ orbitals.  The energy bands at any
other $\kk$-points are calculated by diagonalizing \equ{h_kw}, noting
that the size of the supercell has been chosen so that the Hamiltonian
matrix elements on the right-hand-side of this equation are non
negligible only for WFs up to neighboring supercells $\RR^{(\pm)}$ on
either side of $\RR=0$.  Fig.~\ref{fig:nanotube-bands} shows as solid
lines the interpolated bands, unfolded onto the 20-atom primitive
cell.  Even if with this sampling the system has a pseudogap of ~2eV,
the metallic character of the bands is perfectly reproduced, and these
are in excellent agreement with the bands calculated directly on the
primitive cell by direct diagonalization of the Kohn-Sham Hamiltonian
in the full plane-wave basis set (open circles). The vertical dashed
lines indicate the equivalent first-principles mesh obtained by
unfolding the $\Gamma$-point \footnote{When $\Gamma$ sampling is 
used, special care should be used in calculating matrix elements
between WFs, since the center of a periodic image of e.g. the ket
could be closer to the bra that the actual state considered. Similar
considerations apply for transport calculations, and might require 
calculation of the matrix elements in real-space\cite{lee-phd06}.}.
% 
\begin{figure}
\begin{center}
\includegraphics[width=3.2in]{figs/nanotube-bands}
\caption{Left panel: band structure of a (5,5) carbon nanotube,
  calculated by Wannier interpolation (solid lines), or from a full
  diagonalization in a planewave basis set (circles).  The five
  vertical dashed lines indicate the five $\kk$-points corresponding
  to the $\Gamma$ point in a 100-atom supercell.  The middle and right
  panels show the Wannier-based calculation of the quantum conductance
  and the density of states (see Sec. VII), with a perfect match of
  steps and peaks with respect to the exact band structure. From
  \textcite{lee-prl05}.}
\label{fig:nanotube-bands}
\end{center}
\end{figure}

\subsubsection{GW quasiparticle bands 
%of SrTiO$_3$
}
\label{sec:gw-bands}

In the two examples above the WFs were generated from Kohn-Sham Bloch
functions, and the eigenvalues used in \equ{ham-R} were the
corresponding Kohn-Sham eigenvalues. Many of the deficiencies of the
Kohn-Sham energy bands, such as the underestimation of the energy gaps
of insulators and semiconductors, can be corrected using many-body
perturbation theory in the form of the GW approximation
(for a review, see \textcite{aryasetiawan-ropp98}).

One practical difficulty in generating GW band structure plots is that
the evaluation of the quasiparticle (QP) corrections to the
eigenenergies along different symmetry lines in the BZ is
computationally very demanding. At variance with the DFT formalism,
where the eigenenergies at an arbitrary $\kk$ can be found starting
from the self-consistent charge density, the evaluation of the QP
corrections at a given $\kk$ requires a knowledge of the Kohn-Sham
eigenenergies and wavefunctions on a homogeneous grid of points
containing the wavevector of interest.  What is often done instead is
to perform the GW calculation at selected $\kk$-points only, and then
deduce a ``scissors correction,'' i.e., a constant shift to be applied
to the conduction-band Kohn-Sham eigenvalues elsewhere in the
Brillouin zone.

As already mentioned briefly in Sec.~\ref{sec:manybody},
\textcite{hamann-prb09} proposed using Wannier interpolation to
determine the GW QP bands very efficiently and accurately at arbitrary
points in the BZ. The wannierization and interpolation procedures are
identical to the DFT case. The only difference is that the starting
eigenenergies and overlaps matrices over the uniform first-principles
mesh are now calculated at the GW level.  (In the simplest $G_0W_0$
approximation, where only the eigenenergies, not the eigenfunctions,
are corrected, the wannierization is done at the DFT level, and the
resulting transformation matrices are then applied to the corrected QP
eigenenergies.)

% \dvm{I converted Don's figure into a new version srti-gw-new.eps that
% I think should be free of the earlier problems.}

\begin{figure}
\begin{center}
\includegraphics[width=3.2in]{figs/srtio3-gw-new}
\caption{(Color online) Wannier-interpolated upper valence and lower conduction bands of 
SrTiO$_3$ from LDA (solid red) and GW (dashed blue) calculations. The open
circles at the symmetry points denote the exact GW results taken directly
from the first-principles calculation. From \textcite{hamann-prb09}.}
\label{fig:srtio3-gw}
\end{center}
\end{figure}

Figure~\ref{fig:srtio3-gw} shows a comparison between the interpolated
GW (dashed lines) and DFT-LDA (solid lines) bands of
SrTiO$_3$~\cite{hamann-prb09}. Note that the dashed lines pass through
the open circles at the symmetry points, which denote exact (non-interpolated)
GW results.

Among the recent applications of the GW+Wannier method, we mention the
study of the energy bands of zircon and hafnon~\cite{shaltaf-prb09},
and a detailed comparative study between the DFT-LDA,
scissors-shifted, and QP $G_0W_0$ bands of Si and Ge
nanowires~\cite{peelaers-prb11}. In the latter study the authors found
that the simple scissors correction to the DFT-LDA bands is accurate
near the $\Gamma$ point only, and only for bands close to the highest
valence and lowest conduction band. \textcite{kioupakis-prb10} used
the method to elucidate the mechanisms responsible for free-carrier
absorption in GaN and in the In$_{0.25}$Ga$_{0.75}$N alloy, and
\textcite{yazyev-arxiv11} investigated the
quasiparticle effects on the band structure of the topological
insulators Bi$_2$Se$_3$ and Bi$_2$Te$_3$.


\subsubsection{Surface bands of topological insulators}
\label{sec:surface-bands}

Topological insulators (TIs) were briefly discussed in
Sec.~\ref{sec:top} (see \textcite{hasan-rmp10,hasan-arcmp11} for
useful reviews).  Here we focus on the non-magnetic variety, the ${\cal
  Z}_2$-odd TIs. The recent flurry of activity on this class of
materials has been sustained in part by the experimental confirmation
of the ${\cal Z}_2$-odd character of certain quantum-well structures
and of a rapidly increasing number of bulk crystals.

In the case of 3D TIs, the clearest experimental signature of the
${\cal Z}_2$-odd character is at present provided by ARPES
measurements of the surface electron bands. If time-reversal symmetry
is preserved at the surface, ${\cal Z}_2$-odd materials possess
topologically-protected surface states which straddle the bulk gap,
crossing the Fermi level an odd number of times. These surface states
are doubly-degenerate at the time-reversal-invariant momenta of the
surface BZ, and in the vicinity thereof they disperse linearly,
forming Dirac cones.

First-principles calculations of the surface states for known and
candidate TI materials are obviously of great interest for comparing
with ARPES measurements.  While it is possible to carry out a direct
first-principles calculation for a thick slab in order to study the
topologically protected surface states, as done by
\textcite{yazyev-prl10}, such an approach is computationally
expensive. \textcite{zhang-njp10} used a simplified Wannier-based
approach which suceeds in capturing the essential features of the
topological surface states at a greatly reduced computational cost.
%
% \dvm{\tcb{This margin comment moved from Ch.~3.  Jonathan had written:\\
% JY: Move this text elsewhere. Ch2 or Ch 6?.
%   A further advantage of using projected wannier functions with atomic
%   symmetry is that spin-orbit coupling can be easily added to the basis
%   using values of the spin-orbital splitting taken from atomic
%   data. This approach has been used by \textcite{zhang-njp10} to create
%   an effective Hamiltonian for use in the study of the strong topological
%   insulators Bi$_2$Se$_3$, Bi$_2$Te$_3$, and Sb$_2$Te$_3$
%   (see Sec.~\ref{sec:surface-bands}).\\
% DV: Revised text around here.  Is it OK?}}
%
Their procedure is as follows.  First, an inexpensive calculation is
done for the bulk crystal without spin-orbit interaction.
Next, disentangled WFs spanning the upper valence and low-lying
conduction bands are generated, and the corresponding TB
Hamiltonian matrix is constructed.
The TB Hamiltonian is then augmented with spin-orbit couplings
$\lambda\bf L\cdot S$, where $\lambda$ is taken
from atomic data; this is possible because the WFs have been
constructed so as to have specified $p$-like characters.
The augmented TB parameters are then used to construct
sufficiently thick free-standing ``tight-binding slabs'' by a simple
truncation of the effective TB model, and the dispersion relation is
efficiently calculated by interpolation as a function of the
wavevector $\kk_\parallel$ in the surface BZ.

It should be noted that
this approach contains no surface-specific information, being based
exclusively on the bulk WFs. Even if its accuracy is questionable,
however, this method is useful for illustrating the ``topologically
protected'' surface states that arise as a manifestation of the bulk
electronic structure~\cite{hasan-rmp10}.

Instead of applying the naive truncation,
it is possible to refine the procedure so as to incorporate the
changes to the TB parameters near the surface. To do so, the bulk
calculation must now be complemented by a calculation on a thin slab,
again followed by wannierization. Upon aligning the on-site energies
in the interior of this slab to the bulk values, the changes to the TB
parameters near the surface can be inferred.  However,
\textcite{zhang-prl11b} found that the topological surface states are
essentially the same with and without this surface correction.

The truncated-slab approach was applied by \textcite{zhang-njp10} to
the stoichiometric three-dimensional TIs Sb$_2$Te$_3$, Bi$_2$Se$_3$,
and Bi$_2$Te$_3$.  The calculations on Bi$_2$Se$_3$ revealed the
existence of a single Dirac cone at the $\overline{\Gamma}$ point as
shown in Fig.~\ref{fig:BiSe-slab}, in agreement with ARPES
measurements~\cite{xia-np09}.

\begin{figure}
\begin{center}
\includegraphics[width=6.5cm]{figs/BiSe-slab}
\caption{Wannier-interpolated energy bands of a free-standing (111)
  slab containing 25 quintuple layers of
  Bi$_2$Se$_3$~\cite{zhang-njp10}, plotted along the
  $\overline{\Gamma}$--$\overline{\rm K}$ line in the surface BZ.  A
  pair of topologically-protected surface bands can be seen emerging
  from the dense set of projected valence and conduction bulk-like
  bands and crossing at the time-reversal-invariant point
  $\overline{\Gamma}$. Adapted from \textcite{zhang-njp10}.}
\label{fig:BiSe-slab}
\end{center}
\end{figure}

%\ism{\textcite{zhang-njp10} did not use maxloc WFs.  They only carried
%  out the projection step, to ensure atomic-like WFs (instead, in
%  \textcite{zhang-np09} they appear to be using (spinor?) maxloc WFs,
%  but don't provide many details). I think the reason was that they
%  added SOC on top of the Wannier representation, and so they needed
%  atomic orbitals to evaluate the SOC matrix elements
%  analytically. Maybe we can mention this work in either Sec.~II
%  Sec.~III.B.2, as an example of a situation where simple projection
%  is preferable to full-blown maxloc, as it yields quasiatomic
%  orbitals. \textcolor{red}{Jonathan, what do you think?}}


An alternative strategy for calculating the surface bands was used
earlier by the same authors~\cite{zhang-np09}. Instead of explicitly
diagonalizing the Wannier-based Hamiltonian $H(\kk_\parallel)$ of a
thick slab, the Green's function for the semi-infinite crystal as a
function of atomic plane is obtained via iterative
methods~\cite{lopezsancho-jpf84,lopezsancho-jpf85}, using the
approach of \textcite{lee-prl05}.  Here the
localized Wannier representation is used to break down the
semi-infinite crystal into a stack of ``principal layers'' consisting
of a number of atomic planes, such that only nearest-neighbor
interactions between principal layers exist (see Ch.~\ref{sec:basis}
for more details).

Within each principal layer one forms, starting from the
fully-localized WFs, a set of hybrid WFs
%,'' or ``layer Bloch states,''
which are extended (Bloch-like) along the surface but remain localized
(Wannier-like) in the surface-normal direction (see
Secs.~\ref{sec:hybrid} and \ref{sec:layer}).  This is achieved by
carrying out a partial Bloch sum over the in-plane lattice vectors,
%
\beq 
\label{eq:hybrid-blochsum}
\ket{l,n \kk_\parallel}=\sum_{\RR_\parallel}\,
e^{i\kk_\parallel\cdot\RR_\parallel}\ket{\RR n}, 
\eeq
%
where $l$ labels the principal layer, $\kk_\parallel$ is the in-plane
wavevector, and $\RR_\parallel$ is the in-plane component of
$\RR$. The matrix elements of the Green's function in this basis are
%
\beq
G_{ll'}^{nn'}(\kk_\parallel,\epsilon)=\bra{\kk_\parallel ln}
\frac{1}{\epsilon-H}\ket{\kk_\parallel l'n'}.
\eeq
%
The nearest-neighbor coupling between principal layers means that for
each $\kk_\parallel$ the Hamiltonian has a block tri-diagonal form
(the dependence of the Hamiltonian matrix on $\kk_\parallel$ is given
by the usual Fourier sum expression).  This feature can be exploited
to calculate the diagonal elements of the Green's function matrix very
efficiently using iterative
schemes~\cite{lopezsancho-jpf84,lopezsancho-jpf85}.\footnote{ A
  pedagogical discussion, where a continued-fractions expansion is
  used to evaluate the Green's function of a semi-infinite linear
  chain with nearest-neighbor interactions, is given by
  \textcite{grosso-book00}.}
%
From these, the density of states (DOS) projected onto a given atomic
plane ${\cal P}$ can be obtained\cite{grosso-book00} as

\beq 
N_l^{\cal P}(\kk_\parallel,\epsilon)=-\frac{1}{\pi}\,{\rm
  Im}\sum_{n\in {\cal P}}\, G_{ll}^{nn}(\kk_\parallel,\epsilon+i\eta), 
\eeq
%
where the sum over $n$ is restricted to the orbitals ascribed to the
chosen plane and $\eta$ is a positive infinitesimal.

The projection of the DOS onto the outermost atomic plane is shown in
Fig.~\ref{fig:SbTe-SDOS} as a function of energy $\epsilon$ and
momentum $\kk_\parallel$ for the (111) surface of Sb$_2$Te$_3$.  The
same method has been used to find the dispersion of the surface bands
in the TI alloy Bi$_{1-x}$Sb$_x$~\cite{zhang-prb09} and in ternary
compounds with a honeycomb lattice~\cite{zhang-prl11b}.

%
\begin{figure}
\begin{center}
\includegraphics[width=8.0cm]{figs/SbTe-SDOS-cropped}
\caption{(Color online) Surface density-of-states (SDOS) of a
  semi-infinite crystal of Sb$_2$Te$_3$ terminated with a [111]
  surface~\cite{zhang-np09}. Warmer colors represent a higher
  SDOS. The surface states can be clearly seen around
  $\overline{\Gamma}$ as red lines dispersing in the bulk gap.  From
  \textcite{zhang-np09}.}
\label{fig:SbTe-SDOS}
\end{center}
\end{figure}


\subsection{Band derivatives}
\label{sec:band-der}

The first and second derivatives of the energy eigenvalues with
respect to $\kk$ (band velocities and inverse effective masses) appear
in a variety of contexts, such as the calculation of transport
coefficients~\cite{ashcroft-book76,grosso-book00}.  There is therefore
considerable interest in developing simple and accurate procedures for
extracting these parameters from a first-principles band structure
calculation.

A direct numerical differentiation of the eigenenergies calculated on
a grid is cumbersome and becomes unreliable near band crossings.  It
is also very expensive if a Brillouin zone integration is to be
carried out, as in transport calculations. A number of efficient
interpolation schemes, such as the method implemented in the {\tt
  BoltzTraP} package~\cite{madsen-ccp06}, have been developed for this
purpose, but they are still prone to numerical instabilities near band
degeneracies~\cite{uehara-prb00}.  Such instabilities can be avoided
by using a tight-binding parametrization to fit the first-principles
band structure~\cite{schulz-prb92,mazin-prb00}.  As shown by
\textcite{graf-prb92} and \textcite{boykin-prb95}, both the first and
the second derivatives are easily computed within tight-binding
methods, even
in the presence of band degeneracies, and the same can
be done in a first-principles context using WFs. 

Let us illustrate the procedure by calculating the band gradient away
from points of degeneracy; the treatment of degeneracies and second
derivatives is given in \textcite{yates-prb07}. The first step is to
take analytically the derivative $\partial_\alpha=\partial/\partial
k_\alpha$ of \equ{h_kw},
%
\beq
\label{eq:H_alpha}
%\left[ 
H^\ww_{\kk,\alpha}
\equiv \partial_\alpha H^\ww_\kk
%\right]_{nm}
=\sum_\RR\,e^{i\kk\cdot\RR}iR_\alpha\bra{\0}H\ket{\RR}.
\eeq
%
The actual band gradients are given by the diagonal elements of the
rotated matrix,
%
\beq
\label{eq:band-grad}
\partial_\alpha\enkh=
\left[
  U_\kk^\dagger H_{\kk,\alpha}^\ww U_\kk
\right]_{nn}
\eeq
% 
where $U_\kk$ is the same unitary matrix as in \equ{h_k-diag}.

It is instructive to view the columns of $U_\kk$ as orthonormal state
vectors in the $\nwann$-dimensional ``tight-binding space'' defined by
the WFs.  According to Eq.~(\ref{eq:h_k-diag}) the $n$-th column
vector, which we shall denote by $\vt\phi_{n\kk}\ra$, satisfies the
eigenvalue equation
$H_\kk^\ww\vt\phi_{n\kk}\ra=\enkh\vt\phi_{n\kk}\ra$.  Armed with this
insight, we now recognize in Eq.~(\ref{eq:band-grad}) the
Hellmann-Feynman result $\partial_\alpha\enkh= \la \phi_{n\kk}\vt
\partial_\alpha H_\kk^\ww
\vt \phi_{n\kk}\ra$.


\subsubsection{Application to transport coefficients
}
\label{sec:low-field-hall}

Within the semiclassical theory of
transport, the electrical and
thermal conductivities of metals and doped semiconductors can be
calculated from a knowledge of the band derivatives and relaxation
times $\tau_{n\kk}$ on the Fermi surface.  An example is the low-field
Hall conductivity $\sigma_{xy}$ of non-magnetic cubic metals, which in
the constant relaxation-time approximation is independent of $\tau$
and takes the form of a Fermi-surface integral contaning the first and
second band derivatives~\cite{hurd-book72}.

Previous first-principles calculations of $\sigma_{xy}$ using various
interpolation schemes encountered difficulties for materials such as
Pd, where band crossings occur at the Fermi level~\cite{uehara-prb00}.
A Wannier-based calculation free from such numerical instabilities was
carried out by \textcite{yates-prb07}, who obtained
carefully-converged values for $\sigma_{xy}$ in Pd and other cubic
metals.

A more general formalism to calculate the electrical conductivity
tensor in the presence of a uniform magnetic field involves
integrating the equations of motion of a wavepacket under the field to
find its trajectory on the Fermi surface~\cite{ashcroft-book76}.  A
numerical implementation of this approach starting from the
Wannier-interpolated first-principles bands was carried out by
\textcite{liu-prb09}. This formalism is not restricted to cubic
crystals, and the authors used it to calculate the Hall conductivity
of hcp Mg~\cite{liu-prb09} and the magnetoconductivity of
MgB$_2$~\cite{yang-prl08}.

Wannier interpolation has also been used to determine the Seebeck
coefficient in hole-doped LaRhO$_3$ and CuRhO$_2$~\cite{usui-jpcm09},
in electron-doped SrTiO$_3$~\cite{usui-prb10}, in 
SiGe nanowires~\cite{shelley-epl11}, and in ternary 
skutterudites~\cite{volja-arxiv11}.

\subsection{Berry curvature and anomalous Hall conductivity}
\label{sec:curv-ahc}

The velocity matrix elements between Bloch eigenstates take the
form~\cite{blount-ssp62}
%
\beq
\label{eq:velocity}
\bra{\psi_{n\kk}}\hbar v_\alpha\ket{\psi_{m\kk}}=
\delta_{nm}
\partial_\alpha\enk
-i(\epsilon_{m\kk}-\enk)\left[A_{\kk,\alpha}\right]_{nm},
\eeq
%
where 
%
\beq
\label{eq:berry-connection-matrix}
\left[A_{\kk,\alpha}\right]_{nm}=
i\bra{u_{n\kk}}\partial_\alpha u_{m\kk}\rangle
\eeq
%
is the matrix generalization of the Berry
connection of \equ{berry-connection}.

In the examples discussed in the previous section the static transport
coefficients could be calculated from the first term in
\equ{velocity}, the intraband velocity.
The second term describes vertical interband transitions, which
dominate the optical spectrum of crystals over a wide frequency
range. Interestingly, under certain conditions, {\it virtual}
interband transitions also contribute to the dc Hall
conductivity. This so-called {\it anomalous Hall effect} occurs in
ferromagnets from the combination of exchange splitting and spin-orbit
interaction. For a recent review, see \textcite{nagaosa-rmp10}.

In the same way that WFs proved helpful for evaluating $\partial_\kk
\enk$, they can be useful for calculating quantities containing
$\kk$-derivatives of the cell-periodic Bloch states, such as the Berry
connection of \equ{berry-connection-matrix}. A number of properties
are naturally expressed in this form. In addition to the interband
optical conductivity and the anomalous Hall conductivity (AHC), other
examples include the electric polarization (Sec.~\ref{sec:pol-theo})
as well as the orbital magnetization and magnetoelectric coupling
(Sec.~\ref{sec:magnetism}).

Let us focus on the Berry curvature ${\bm{\mathcal{F}}}_{n\kk}$
[\equ{curv-pot}], a quantity with profound effects on the dynamics of
electrons in crystals~\cite{xiao-rmp10}.  ${\bm{\mathcal{F}}}_{n\kk}$
can be nonzero if either spatial inversion or time-reversal symmetries
are broken in the crystal, and when present acts as a kind of
``magnetic field'' in $\kk$-space, with the Berry connection ${\bf
  A}_{n\kk}$ playing the role of the vector potential. This effective
field gives rise to a Hall effect in ferromagnets even in the absence
of an actual applied ${\bf B}$-field (hence the name
``anomalous''). The %intrinsic
AHC is given by\footnote{\equa{ahc} gives the so-called ``intrinsic''
  contribution to the AHC, while the measured effect also contains
  ``extrinsic'' contributions associated with scattering from
  impurities~\cite{nagaosa-rmp10}.}
%
\beq
\label{eq:ahc}
\sigma_{\alpha\beta}^{\rm AH}=-\frac{e^2}{\hbar}\int_{\rm BZ}
\frac{d\kk}{(2\pi)^3}\Omega^{\rm tot}_{\kk,\alpha\beta},
\eeq
%
where
$\Omega^{\rm tot}_{\kk,\alpha\beta}=\sum_n\,f_{n\kk}\Omega_{n\kk,\alpha\beta}$
($f_{n\kk}$ is the Fermi-Dirac distribution function), and we have
rewritten the pseudovector ${\bm{\mathcal{F}}}_{n\kk}$ as an
antisymmetric tensor. 

The interband character of the intrinsic AHC can be seen by using
$\kk\cdot{\bf p}$ perturbation theory to write the $\kk$-derivatives in
\equ{curv-pot}, leading to the ``sum-over-states'' formula
%
\beq
\label{eq:berry-curv-sum}
\Omega^{\rm tot}_{\alpha\beta}=i\sum_{n,m}\,(f_n-f_m)
\frac{\bra{\psi_n}\hbar v_\alpha\ket{\psi_m}
      \bra{\psi_m}\hbar v_\beta\ket{\psi_n}
     }
     {(\epsilon_m-\epsilon_n)^2}.
\eeq
%
The AHC of bcc Fe and SrRuO$_3$ was evaluated from the previous two
equations by \textcite{yao-prl04} and \textcite{fang-science03}
respectively. These pioneering first-principles calculations revealed
that in ferromagnetic metals the Berry curvature displays strong and
rapid variations in $\kk$-space (Fig.~\ref{fig:Fe-curv-spikes}). As a
result, an ultradense BZ mesh containing millions of $\kk$-points is
often needed in order to converge the calculation.  This is the kind
of situation where the use of Wannier interpolation can be most
beneficial.

\begin{figure}
\begin{center}
\includegraphics[width=3.2in]{figs/Fe-curv-spikes}
\caption{Upper panel: band structure near the Fermi level of bcc Fe
  with the spontaneous magnetization along $\hat{\bf z}$. Lower panel:
  the Berry curvature summed over the occupied bands
  [\equ{berry-curv-sum}], plotted along the same symmetry lines. The
  sharp spikes occur when two spin-orbit-coupled bands are separated
  by a small energy across the Fermi level, producing a resonance
  enhancement. Adapted from \textcite{yao-prl04}.}
\label{fig:Fe-curv-spikes}
\end{center}
\end{figure}

The strategy for interpolating the Berry curvature is similar to that
used in Sec.~\ref{sec:band-der} for the band gradient.  One first
evaluates certain objects in the Wannier gauge using Bloch sums, and
then transform to the Hamiltonian gauge. Because the gauge
transformation mixes the bands, it is convenient to introduce a
generalization of \equ{curv-pot} having two band indices instead of
one. To this end we start from \equ{berry-connection-matrix} and
define the matrices
%
\beq
\label{eq:berry-curv-mat}
\Omega_{\alpha\beta}=
\partial_\alpha A_\beta-\partial_\beta A_\alpha=
i\bra{\partial_\alpha u}\partial_\beta u\rangle-
i\bra{\partial_\beta u}\partial_\alpha u\rangle,
\eeq
% 
where every object in this expression should consistently carry either
an H or W label.  Provided that the chosen WFs correctly span all
occupied states, the integrand of \equ{ahc} can now be expressed as
$\Omega^{\rm tot}_{\alpha\beta}=
\sum_{n=1}^\nwann\,f_n\Omega^\hh_{\alpha\beta,nn}$.

A useful expression for $\Omega^\hh_{\alpha\beta}$ can be obtained
with the help of the gauge-transformation law for the Bloch states,
$\ket{u_\kk^\hh}=\ket{u_\kk^\ww}U_{\kk}$
[\equ{psik-h}]. Differentiating both sides with respect to $k_\alpha$
and then inserting into Eq.~(\ref{eq:berry-curv-mat}) yields, after a
few manipulations,
%
\beq
\label{eq:omega-h}
\Omega^\hh_{\alpha\beta}=\overline{\Omega}_{\alpha\beta}
-\left[ D_\alpha,\overline{A}_\beta \right]
+\left[ D_\beta,\overline{A}_\alpha \right]
-i\left[ D_\alpha,D_\beta \right],
\eeq
% 
where $D_\alpha=U^\dagger\partial_\alpha U$, and
$\overline{A}_\alpha$, $\overline{\Omega}_{\alpha\beta}$ are related
to the connection and curvature matrices in the Wannier gauge through
the definition $\overline{O}_\kk=U^\dagger_\kk O_\kk^\ww U_\kk$. Using
the band-diagonal elements of \equ{omega-h} in the expression for
$\Omega^{\rm tot}_{\alpha\beta}$ eventually leads to
%
\beq
\begin{split}
\label{eq:ahc-wanint}
\Omega^{\rm tot}_{\alpha\beta}=\sum_n^\nwann\, 
& f_n\overline{\Omega}_{\alpha\beta,nn}
+\sum_{mn}^\nwann\,(f_m-f_n)
\big(
        D_{\alpha,nm}\overline{A}_{\beta,mn} \\&
       -D_{\beta,nm}\overline{A}_{\alpha,mn}
       +iD_{\alpha,nm}D_{\beta,mn}
\big).
\end{split}
\eeq
%
This is the desired expression, which in the Wannier interpolation
scheme takes the place of the sum-over-states formula. In contrast to
\equ{berry-curv-sum}, note that the summations over bands now run over
the small set of Wannier-projected bands.

The basic ingredients going into \equ{ahc-wanint} are the Wannier
matrix elements of the Hamiltonian and of the position operator. From
a knowledge of $\bra{\0}H\ket{\RR}$ the energy eigenvalues and
occupation factors, as well as the matrices $U$ and $D_\alpha$, can be
found using band-structure interpolation
(Sec.~\ref{sec:band-interpolation}). The information about
$A^\ww_\alpha$ and $\Omega^\ww_{\alpha\beta}$ is instead encoded in
the matrix elements $\bra{\0}\rr\ket{\RR}$, as can be seen by
inverting \equ{rmatel},
%
\beq
\label{eq:A-w}
%\left[
A^\ww_{\alpha}
=\sum_\RR\,e^{i\kk\cdot\RR}\bra{\0}r_\alpha\ket{\RR}.
\eeq
%
As for $\Omega^\ww_{\alpha\beta}$, according to \equ{berry-curv-mat}
it is given by the curl of this expression, which can be taken
analytically.

The strategy outlined above was demonstrated by \textcite{wang-prb06}
in calculating the AHC of bcc Fe, using the spinor WFs of
Sec.~\ref{sec:fe-bands}. Both the $\kk$-space distribution of the Berry
curvature and the integrated AHC were found to be in excellent
agreement with the sum-over-states calculation of
\textcite{yao-prl04}.

Table~\ref{table:ahc} lists the AHC of the ferromagnetic transition
metal elements, calculated with the magnetization along the respective
easy axes. The magnetic anisotropy of the AHC was investigated by
\textcite{roman-prl09}.  While the AHC of the cubic metals Fe and Ni
is fairly isotropic, that of hcp Co was found to decrease by a factor
of four as the magnetization is rotated from the $c$-axis to the basal
plane.  The Wannier method has also been used to calculate the AHC in
FePt and FePd ordered alloys~\cite{seeman-prl10,zhang-prl11}, and the
spin-Hall conductivity in a number of metals~\cite{freimuth-prl10}.

\begin{table}
  \caption{
    Anomalous Hall conductivity in S/cm of the ferromagnetic
    transition metals, calculated from first-principles with the 
    magnetization along the respective easy axes.
    The first two rows show values obtained using the Wannier interpolation
    scheme
    to either integrate the Berry curvature over the Fermi 
    sea (see main text), or to
    evaluate the Berry phases of planar loops around the Fermi 
    surface~\cite{wang-prb07}. 
    Results obtained using the sum-over-states formula, \equ{berry-curv-sum},
    are included for comparison, as well as representative experimental
    values. Adapted from \textcite{wang-prb07}.}
\begin{ruledtabular}
\begin{tabular}{lrrr}
& bcc Fe & fcc Ni & hcp Co \cr
\hline
Berry curvature & 753\,\, & $-$2203\,\, & 477\,\, \cr
Berry phase & 750\,\, & $-$2275\,\, & 478\,\, \cr
Sum-over-states & 751\footnotemark[1]
    & $-$2073\footnotemark[2] & 492\footnotemark[2] \cr
Experiment& 1032  & $-$646 &480 \cr
\end{tabular}
\end{ruledtabular}
\footnotetext[1]{\textcite{yao-prl04}.}
\footnotetext[2]{Y. Yao, private communication.}
\label{table:ahc}
\end{table}

As already mentioned, for certain applications the Berry connection
matrix [\equ{berry-connection-matrix}] is the object of direct
interest. The interpolation procedure described above can be directly
applied to the off-diagonal elements describing vertical interband
transitions, and the magnetic circular dichroism spectrum of bcc Fe
has been determined in this way~\cite{yates-prb07}.
 
The treatment of the diagonal elements of the Berry connection matrix
is more subtle, as they are locally gauge-dependent. Nevertheless, the
Berry phase obtained by integrating over a closed loop in $\kk$-space,
$\varphi_n=\oint {\bf A}_{n\kk}\cdot d{\bf l}$, is
gauge-invariant~\cite{xiao-rmp10}.  Recalling that
${\bm{\mathcal{F}}}_{n\kk}={\nabla}_\kk\times {\bf
  A}_{n\kk}$ [\equ{curv-pot}] and using Stokes' theorem, \equ{ahc} for
the AHC can be recast in terms of the Berry phases of Fermi loops on
planar slices of the Fermi surface. This approach has been implemented by
\textcite{wang-prb07}, using Wannier interpolation to sample
efficiently the orbits with the very high density required near
band-crossings.  Table~\ref{table:ahc} lists values for the AHC
calculated using both the Berry curvature (``Fermi sea'') and
Berry-phase (``Fermi surface'') approaches. Similar convergence issues
that appear in calculating Knight shifts and orbital magnetization for
metals\cite{davezac-prb07} could be addressed very proficiently by Wannier interpolation.

\subsection{Electron-phonon coupling}
\label{sec:el-ph}

The electron-phonon interaction~\cite{grimvall-book81} plays a key
role in a number of phenomena, from superconductivity to the
resistivity of metals and the temperature dependence of the optical
spectra of semiconductors.  The matrix element for scattering an
electron from state $\psi_{n\kk}$ to state $\psi_{m,\kk+\qq}$ while
absorbing a phonon $\qq\nu$ is proportional to the electron-phonon
vertex
%
\beq
\label{eq:g_kq}
g_{\nu,mn}(\kk,\qq)=\bra{\psi_{m,\kk+\qq}}\partial_{\qq\nu} V\ket{\psi_{n\kk}}.
\eeq
%
Here $\partial_{\qq\nu} V$ is the derivative of the self-consistent
potential with respect to the amplitude of the phonon with branch
index $\nu$ and momentum $\qq$.  Evaluating this vertex is a key task
for a first-principles treatement of electron-phonon couplings.

State-of-the-art calculations using first-principles linear-response
techniques~\cite{baroni-rmp01} have been successfully applied to a
number of problems, starting with the works of
  \textcite{savrasov-prl94} and \textcite{mauri-prl96b}, who used
respectively the LMTO and planewave pseudopotential methods.  The cost of evaluating
\equ{g_kq} from first-principles over a large number of
$(\kk,\qq)$-points is quite high, however, and this has placed a
serious limitation on the scope and accuracy of first-principles
techniques for electron-phonon problems.

The similarity between the Wannier interpolation of energy bands and
the Fourier interpolation of phonon dispersions was already noted.  It
suggests the possibility of interpolating the electron-phonon vertex
in both the electron and the phonon momenta, once \equ{g_kq} has been
calculated on a relatively coarse uniform $(\kk,\qq)$-mesh.  Different
electron-phonon interpolation schemes have been put forth in the
literature~\cite{giustino-prl07,eiguren-prb08,calandra-prb10}
In the following we
describe the approach first developed by
\textcite{giustino-prb07} and implemented in the software package {\tt
  EPW}~\cite{noffsinger-cpc10}.
%
To begin, let us set the notation for lattice
dynamics~\cite{maradudin-rmp68}. We write the instantaneous nuclear
positions as
%
%\beq
$\RR+{\boldsymbol\tau}_s+{\bf u}_{\RR s}(t)$,
%\eeq
% 
where $\RR$ is the lattice vector, ${\boldsymbol\tau}_s$ is the
equilibrium intracell coordinate of ion $s=1,\ldots,p$, and ${\bf
  u}_{\RR s}(t)$ denotes the instantaneous displacement.  The normal
modes of vibration take the form
%
\beq
\label{eq:normal_mode}
{\bf u}^{\qq\nu}_{\RR s}(t)={\bf u}^{\qq\nu}_s
e^{i(\qq\cdot\RR-\omega_{\qq\nu}t)}.
\eeq
%
The eigenfrequencies $\omega_{\qq\nu}$ and mode amplitudes ${\bf
  u}^{\qq\nu}_s$ are obtained by diagonalizing the dynamical matrix
$\left[D_\qq^{\rm ph}\right]_{st}^{\alpha\beta}$, where $\alpha$ and
$\beta$ denote spatial directions.
%
It is expedient to introduce composite indices $\mu=(s,\alpha)$ and
$\nu=(t,\beta)$, and write $D^{\rm ph}_{\qq,\mu\nu}$ for the dynamical
matrix. With this notation, the eigenvalue equation becomes
%
\beq
\left[e_\qq^\dagger D^{\rm ph}_\qq e_\qq\right]_{\mu\nu}=
\delta_{\mu\nu}\omega^2_{\qq\nu},
\eeq
% 
where $e_\qq$ is a $3p\times 3p$ unitary matrix. In analogy with the
tight-binding eigenvectors $\vt \phi_{n\kk}\ra$ of
Sec.~\ref{sec:band-der}, we can view the columns of $e_{\qq,\mu\nu}$
as orthonormal phonon eigenvectors ${\bf e}^{\qq\nu}_s$.  They are
related to the complex phonon amplitudes by ${\bf
  u}^{\qq\nu}_s=(m_0/m_s)^{1/2}{\bf e}^{\qq\nu}_s$ ($m_0$ is a
reference mass), which we write in matrix form as $U^{\rm
  ph}_{\qq,\mu\nu}$.

Returning to the electron-phonon vertex, \equ{g_kq}, we can now write
explicitly the quantity $\partial_{\qq\nu} V$ therein as
%
\beq
\label{eq:derV-q}
\begin{split}
\partial_{\qq\nu} V(\rr)&=\frac{\partial}{\partial\eta}
V(\rr;\{\RR+{\boldsymbol\tau}_s
+\eta{\bf u}^{\qq\nu}_{\RR s} \})\\
&=\sum_{\RR,\mu}\,e^{i\qq\cdot\RR}
\partial_{\RR\mu}V(\rr)
U^{\rm ph}_{\qq,\mu\nu},
\end{split}
\eeq
%
where $\partial_{\RR\mu}V(\rr)$ is the derivative of the
self-consistent potential with respect to $u_{\RR s,\alpha}$.  As will
be discussed in Sec.~\ref{sec:other}, it is possible to view these
single-atom displacements as maximally-localized ``lattice Wannier
functions.''  With this interpretation in mind we define the
Wannier-gauge counterpart of $\partial_{\qq\nu} V(\rr)$
as
%
\beq
\label{eq:g-w}
\partial^\ww_{\qq\mu}V(\rr)=\sum_\RR\,e^{i\qq\cdot\RR}
\partial_{\RR\mu}V(\rr),
\eeq
%
related to the ``eigenmode-gauge''  quantity $\partial_{\qq\nu} V(\rr)$ by
%
\beq
\label{eq:latt-gauge-transf}
\partial_{\qq\nu}V(\rr)=\sum_\mu\,\partial^\ww_{\qq\mu}V(\rr)
U^{\rm ph}_{\qq,\mu\nu}.
\eeq
%
Next we introduce the Wannier-gauge vertex $g^\ww_\mu(\kk,\qq)=
\bra{\psi^\ww_{\kk+\qq}}\partial^\ww_{\qq\mu} V\ket{\psi^\ww_{\kk}}$,
which can be readily interpolated onto an arbitrary point
$(\kk',\qq')$ using Eqs.~(\ref{eq:bloch-w-gauge}) and (\ref{eq:g-w}),
%
\beq
\label{eq:g-intp-w}
g^\ww_\mu(\kk',\qq')=\sum_{\RR_e,\RR_p}\,
e^{i(\kk'\cdot\RR_e+\qq'\cdot\RR_p)}
\bra{\0_e}\partial_{\RR_p\mu}V\ket{\RR_e}.
\eeq
%
The object $\bra{\0_e}\partial_{\RR_p\mu}V\ket{\RR_e}$, the
electron-phonon vertex in the Wannier representation, is depicted
schematically in Fig.~\ref{fig:vertex-w}. Its localization in
real-space ensures that \equ{g-intp-w} smoothly interpolates $g^\ww$
in the electron and phonon momenta. Finally, we transform the
interpolated vertex back to the Hamiltonian/eigenmode gauge,
%
\beq
\label{eq:g-intp-h}
\begin{split}
g^\hh_\nu(\kk',\qq')&=
\bra{\psi^\hh_{\kk'+\qq'}}\partial_{\qq'\nu} V\ket{\psi^\hh_{\kk'}}\\
&=U^\dagger_{\kk'+\qq'}
\Big[
  \sum_\mu g^\ww_\mu(\kk',\qq')
  U^{\rm ph}_{\qq',\mu\nu}
\Big]
U_{\kk'},
\end{split}
\eeq
% 
where we used Eqs.~(\ref{eq:psik-h}) and (\ref{eq:latt-gauge-transf}).

Once the matrix elements $\bra{\0_e}\partial_{\RR_p\mu}V\ket{\RR_e}$
are known, the electron-phonon vertex can be evaluated at arbitrary
$(\kk',\qq')$ from the previous two equations.  The procedure for
evaluating those matrix elements is similar to that leading up to
\equ{ham-R} for $\bra{\0}H\ket{\RR}$: invert \equ{g-intp-w} over the
first-principles mesh, and then use Eqs.~(\ref{eq:psik-w}) and
(\ref{eq:latt-gauge-transf}).  
%
%\ism{Maybe cut the text in blue?\\
%DV: I guess yes...}
%
%\textcolor{blue}{The resulting expression~\cite{giustino-prb07}
%  contains the vertex (\ref{eq:g_kq}) calculated from first-principles
%  over the {\it ab initio} grid, and the set of matrices ${\cal
%    U}_\kk$ and $\widetilde{U}^{\rm ph}_\qq$ describing the
%  Wannierization of the electronic and ionic degrees of freedom.}

\begin{figure}
\begin{center}
\includegraphics[width=6.0cm]{figs/wannier-vertex-crop}
\caption{(Color online) Real-space representation of
  $\bra{\0_e}\partial_{\RR_p\mu}V\ket{\RR_e}$, the electron-pohnon
  vertex in the Wannier basis. The black squares denote the crystal
  lattice, the blue lines the electron Wannier functions $\ket{{\bf
      0}_e}$ and $\ket{\RR}_e$, and the red line the phonon
  perturbation in the lattice Wannier representation,
  $\partial_{\RR_p}V(\rr)$. Whenever two of these functions are
  centered on distant unit cells, the vertex becomes vanishingly
  small. From \textcite{giustino-prb07}.}
\label{fig:vertex-w}
\end{center}
\end{figure}

The above interpolation scheme has been applied to a number of
problems, including the estimation of $T_c$ in
superconductors~\cite{giustino-prl07,noffsinger-prb08,noffsinger-prb09};
the phonon renormalization of energy bands near the Fermi level in
graphene~\cite{park-prl07} and copper oxide
superconductors~\cite{giustino-n08}; the phonon renormalization of the
band gap of diamond~\cite{giustino-prl10}; and the vibrational 
lifetimes~\cite{park-nl08}
and electron linewidths~\cite{park-prl09}
arising from electron-phonon interactions in graphene.

We mention in closing the work of \textcite{calandra-prb10}, where the
linear-response calculation of the deformation potential
$\partial_{\qq\nu} V$ is carried out taking into account nonadiabatic
effects (that is, going beyond the usual approximation of static ionic
displacements). Using the electron-phonon interpolation scheme of
\textcite{giustino-prb07} to perform the BZ integrations, and a
Wannier interpolation of the dynamical matrix to capture Kohn anomalies, 
the authors
found significant nonadiabatic corrections to the phonon frequencies
and electron-phonon coupling matrix elements in MgB$_2$ and CaC$_6$.

