%% LyX 2.1.2 created this file.  For more info, see http://www.lyx.org/.
%% Do not edit unless you really know what you are doing.
\documentclass[10pt,english]{article}
\usepackage[latin9]{inputenc}
\usepackage[a4paper]{geometry}
\geometry{verbose,tmargin=2cm,bmargin=2cm,lmargin=2cm,rmargin=2cm}
\setlength{\parskip}{\medskipamount}
\setlength{\parindent}{0pt}
\usepackage{color}
\usepackage{babel}
\usepackage{float}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{esint}
\usepackage[unicode=true,pdfusetitle,
 bookmarks=true,bookmarksnumbered=true,bookmarksopen=true,bookmarksopenlevel=1,
 breaklinks=true,pdfborder={0 0 1},backref=section,colorlinks=true]
 {hyperref}
\hypersetup{
 linkcolor=darkred, citecolor=darkgreen, urlcolor=red, filecolor=magenta}
\usepackage{breakurl}

\makeatletter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% User specified LaTeX commands.
\usepackage{html}
\definecolor{darkgray}{rgb}{0.8,0.8,0.8}
\definecolor{darkblue}{rgb}{0,0,0.65}
\definecolor{darkred}{rgb}{0.65,0,0}

\makeatother

\begin{document}

\title{Wannier functions with FPLO}


\author{Klaus Koepernik and Marc H{\"o}ppner}


\date{\today}

\maketitle

\section{Intro}

Wannier functions can be defined in many ways since there is a gauge
freedom of chosing a phase transformation. The ideal way of fixing
the gauge is the requirement of maximum localization. This is a tedious
algorithm and we do not do this in FPLO. However, it turns out that
our definition leads to highly localized Wannier functions (WFs).
Should the Wannier functions turn out not to be localized, it in most
cases means that the WF is badly chosen. The main drawback of our
approach is, that the user has to decide, where the WFs shall sit
and which symmetry they shall have. However, this is at the same time
intended, since for modeling it is exactly what one wants to do.

The WF sitting in cell $R$ and being of type $\mu$ (which denotes
the WF center and its symmetry) is defined as
\begin{equation}
W_{R\mu}=\int\mathrm{e}^{-ikR}\sum_{n}\Psi_{n}^{k}U_{n\mu}^{k}\label{eq:def_WF}
\end{equation}
where $\Psi$ denotes the Kohns-Sham (KS) functions and $U$ is a
unitary matrix. If we define less WFs than KS functions $U$ is a
column unitary projector ($U^{+}U=1$, $UU^{+}=P$). It mapps all
the KS bands on some few WFs. The choice of $U$ is the choice of
the gauge. In FPLO we use a chemically motivated local orbital basis
$\Phi$ to construct 
\[
\Psi_{n}^{k}=\frac{1}{\sqrt{N}}\sum_{Rs\nu}\Phi_{Rs\nu}\mathrm{e}^{ik\left(R+s\right)}C_{s\nu,n}^{k}
\]
($s$ is the atom position and $\nu$ some quantum numbers specifying
the orbital). These orbitals although non-orthogonal are in a way
'optimally' localized by their construction. Hence, it is clear that
a WF centered at an atom and having a certain orbital symmetry has
the corresponding orbital as its main contribution. This allows the
following choice of $U$. We project the KS functions on a test function
$\chi$, which is an FPLO orbital in the simplest case. The resulting
number is the square root of the orbital character of the KS bands,
as plotted in the FAT bands. If we now select only the KS functions
whith a large such orbital character in Eq. (\ref{eq:def_WF}), we
will end up with a WF resembling $\chi$ the most and this is where
the localization of our WFs comes from. If we want WFs corresponding
to sub bands of a band complex with a character $\chi$ we have to
project onto a particular energy window as well. This happens if there
are bonding (B) and antibonding (AB) bands of character $\chi$ and
if one wants, say, only WFs of the AB bands. So, in total we define
a test function or WF projector $\chi$ and an energy window for each
WF. The WF projector $\chi$ can in principle be a linear combination
of FPLO orbitals. Example: suppose we have a cuprate plane. The bands
are formed of linear combinations of Cu $3d$ and ligand O-$2p$ orbitals.
Around certain $k$-points the upper bands are clearly AB and hence
correspond to a certain molecular orbital with a certain phase relation
between the central $3d$ orbital and the O-$2p$ orbitals. The lower
(B) bands are clearly formed of the same orbitals, but with a different
phase relation. If this is the case then defining WF-projectors using
the AB molecular orbital as $\chi$ will automatically project out
the antibonding bands, given that the band topology is dominated by
the clear character separation around the considered k-point. This
would make the energy window obsolete.

The Wannier transformation can be described in two steps. In the first
step the unitary projector is build from the users definitions of
$\chi$ and the energy windows. $U$ is applied to the KS functions
yielding the Bloch sums of the WFs
\begin{eqnarray}
W_{\mu}^{k}\left(r\right) & = & \sum_{n}\Psi_{n}^{k}\left(r\right)U_{n\mu}^{k}\label{eq:WF-Bloch-sums}\\
 & = & \frac{1}{\sqrt{N}}\sum_{R}\mathrm{e}^{ikR}W_{R\mu}\left(r\right)\nonumber 
\end{eqnarray}
\begin{equation}
W_{R\mu}\left(r\right)=\int\mathrm{e}^{-ikR}W_{\mu}^{k}\left(r\right)dk\label{eq:WF-Bloch-sum-Fouriertrans}
\end{equation}
(Note, that the $k$-integrations contain appropriate normalization
factors, left out in the formulas or being hidden in the integral
measure.) The orbtial Bloch functions are
\[
\Phi_{s\mu}^{k}=\frac{1}{\sqrt{N}}\sum_{R}\Phi_{Rs\mu}\mathrm{e}^{ik\left(R+s\right)}
\]
where the additional phase factor $\mathrm{e}^{iks}$ makes live easier
by removing the dependence on the coordinate origin. From the formulas
we get the representation of the Hamiltonian in orbital Bloch sums
\begin{eqnarray*}
\underline{H}_{s^{\prime}s}^{k} & = & \left\langle \Phi_{s^{\prime}}^{k}\hat{H}\Phi_{s}^{k}\right\rangle \\
 & = & \frac{1}{N}\sum_{R^{\prime}R}\left\langle \Phi\hat{H}\Phi\right\rangle _{R^{\prime}s^{\prime}Rs}\mathrm{e}^{ik\left(R+s-R^{\prime}-s^{\prime}\right)}\\
 & = & \sum_{R}\left\langle \Phi\hat{H}\Phi\right\rangle _{0s^{\prime},Rs}\mathrm{e}^{ik\left(R+s-s^{\prime}\right)}
\end{eqnarray*}
from which we get the KS eigenvalues
\begin{equation}
\varepsilon_{n}^{k}=\left(C^{k+}\underline{H}^{k}C^{k}\right)_{n}\label{eq:KS-bandstructure}
\end{equation}
This gives the WF Bloch representation
\begin{eqnarray}
\left\langle W_{\mu^{\prime}}^{q}\hat{H}W_{\mu}^{k}\right\rangle  & = & \sum_{n}U_{n\mu^{\prime}}^{q*}\varepsilon_{n}^{k}U_{n\mu}^{k}\delta_{qk}\label{eq:WF-Bloch_Ham}\\
 & = & \frac{1}{N}\sum_{RP}\left\langle W_{P\mu^{\prime}}\hat{H}W_{R\mu}\right\rangle \mathrm{e}^{ik\left(R-P\right)}\delta_{qk}\nonumber \\
 & = & \sum_{R}\left\langle W_{0\mu^{\prime}}\hat{H}W_{R\mu}\right\rangle \mathrm{e}^{ikR}\delta_{qk}\nonumber 
\end{eqnarray}
with 
\begin{eqnarray}
\varepsilon_{0\mu^{\prime},R\mu} & = & \left\langle W_{0\mu^{\prime}}\hat{H}W_{R\mu}\right\rangle \label{eq:WF-hoppings}\\
 & = & \int\mathrm{e}^{-ikR}\left\langle W_{\mu^{\prime}}^{k}\hat{H}W_{\mu}^{k}\right\rangle dk\nonumber 
\end{eqnarray}
which is the WF Hamiltonian in real space representation, which usualy
contains the model we are interested in. Its $k$-representation (Bloch
sums) Eq. (\ref{eq:WF-Bloch_Ham}) can be diagonalized and will give
the bandstructure corresponding to the model. If the WFs represent
the whole Hilbert space spanned by all $\Psi_{n}^{k}$ of course the
resulting WF-bandstructure coincides with the original bandstructure
Eq. (\ref{eq:KS-bandstructure}). There is a modification one can
do, which consists of restricting the matrix elements Eq. (\ref{eq:WF-hoppings})
by removing hoppings with distances above a certain cut off or hoppings,
which are smaller than a certain threshold. The resulting restricted
hopping matrix can be Bloch summed and diagonalized again. Note, however,
that this modified Hamiltonian does not strictly correspond to the
WFs calculated above.

The whole transformation can be summed up, by defining a transformation
from the FPLO basis into the WF basis
\[
W_{R\mu}\left(r\right)=W_{0\mu}\left(r-R\right)
\]
\begin{equation}
W_{0\mu}\left(r\right)=\sum_{R^{\prime}s\nu}\Phi_{R^{\prime}s\nu}D_{R^{\prime}s\nu,\mu}\label{eq:WF-from-Phi}
\end{equation}



\section{The FPLO WF module}

The Wannier function module in FPLO is currently a postprocessing
tool. First, one needs a converged calculation. Then, the relevant
information must be written to the hard disk in order to access it
conveniently later. This might take plenty of disk space! (See Sec.
\ref{sub:Saving-memory}) This information is then read and used in
a subsequent FPLO run to calculate the desired Wannier functions.
In order to use the module a particular file (\texttt{=.wandef}) must
be created by the user. If this file is found by a running FPLO process
the Wannier module is activated, if the keyword \textbf{doit} is found
on a single line in the file. To disable the module change the keyword
into something like e.g. \textbf{xdoit}. If FPLO finds the file \texttt{=.wandef}
with the keyword it will start dumping data after every Kohn-Sham
diagonalization process. If the FPLO process comes to convergence
(best to start with a converged calculation to begin with) it stops
like in normal mode. The file created after the initial data-dumping
is \texttt{+wancoeff} (besides all the usual files). Note: that the
band structure plot in the fedit submenu should be switched on.

On restart of FPLO (provided that \texttt{=.wandef} exist and the
\textbf{doit}-keyword is set therein) all the data are used and the
actual WF-module is executed. It reads WF definitions from \texttt{=.wandef}
and constructs WFs accordingly. The WFs can be produced in real space
for visualization and the WF Hamiltonian in Bloch representation Eq.
(\ref{eq:WF-Bloch_Ham}) can be written to the file \texttt{+hamongrid}
on a $k$-space grid. The WF hopping itegrals Eq. (\ref{eq:WF-hoppings})
are produced in the output and can be used to extract models. Beware
that the WF hamiltonian on the $k$-space grid is restricted by cutting
of the real-space Hamiltonian Eq. (\ref{eq:WF-hoppings}) according
to user input.

We will walk through an example step by step. The input-files are
in the example directory. The example is CaCuO2 in a ferro-magnetic
calculation (just to have that complication in). We want to create
a WF for the antibonding $3d_{x^{2}-y^{2}}$ band. In Figure \ref{fig:FPLO-fat-bands-CaCuO}
we show the spin-polarized 
\begin{figure}[H]
\begin{centering}
\begin{minipage}[c]{0.25\columnwidth}%
\begin{center}
\includegraphics[width=0.99\columnwidth]{pics/strcu}
\par\end{center}%
\end{minipage}%
\begin{minipage}[c]{0.75\columnwidth}%
\begin{center}
\includegraphics[width=0.99\columnwidth]{pics/fatbands}
\par\end{center}%
\end{minipage}
\par\end{centering}

\protect\caption{\label{fig:FPLO-fat-bands-CaCuO}FPLO fat-bands CaCuO$_{2}$}


\end{figure}
fat-bands as given by FPLO. The filled symbols denote the 3d majority
fat-bands and the open symbols the 3d minority bands. The $3d_{x^{2}-y^{2}}$
bands are called $3d+2$ according to the standard convention and
orientation of the global coordinate system and have magenta color.
Focus on the majority bands (filled symbols) , you see a bonding and
antibonding band (best seen at the M-point). We are interested in
a description of the antibonding band. There is only one Cu atom in
the unit cell and acoording to the list of sites in the output that
is site number 2. We now guess that the Cu, site 2, $3d_{x^{2}-y^{2}}$
orbital will have the largest contribution to a WF describing the
anti-bonding part of the filled magenta band(s). Thus, we chose this
orbital as a WF-projector $\chi$. We assume that the file \texttt{=.wandef}
exists and that the keyword \textbf{doit} is set and that the file
\texttt{+wancoeff} got already created.

We have to define, which linear combinations of local orbitals are
used as WF-projectors. Note, that we have two basic choices: the WF-projector
and the energy window. In our case of CaCuO$_{2}$ we already identified
the projector. The simplest projector consists of one single orbital
centered at on site. A WF-projector is defined using the keyword \textbf{wandef}.
A \textbf{wandef} can have many contributions \textbf{contrib}, each
denoting an FPLO orbital with a certain weight factor. Here comes
the example
\begin{quote}
\ldots{}\\
-------- Cu --------------\\
wandef \\
on \\

\begin{quote}
name Cu x2-y2\\
emin -4 -1 \\
emax -1 3 \\
de 1 1\\
contrib 
\begin{quote}
site 2 \\
difvec 0 0 0 \\
xaxis 1 0 0 \\
zaxis 0 0 1\\
orb 3d+2 \\
fac 1
\end{quote}
\end{quote}
\end{quote}
As can be seen there can be arbitrary lines in the file, and they
are considered comments as long as they do not start with a keyword.
In our case '--------- Cu ----------------' is such a comment. Then
comes the keyword \textbf{wandef} starting of a WF-projector. The
next thing is \textbf{on} or \textbf{off}. This allows to have several
\textbf{wandef}s in one file, where only some are used (playing around).
The \textbf{wandef} has a name, which can be anything after the keyword
\textbf{name} (actually it should not be more than 17 characters).
The energy window is discussed later. in our example the \textbf{wandef}
is made of one \textbf{contrib}. The \textbf{contrib} is an orbital
sitting at \textbf{site} 2, being \textbf{difvec} 0 0 0 away from
the (imagined) centerpoint of the WF, which in this case is the position
of site 2. The orbital is defined with respect to the local coordinate
system, whose \textbf{xaxis} is 1 0 0 and \textbf{zaxis} is 0 0 1.
In this local system the \textbf{orb} is 3d+2 and the \textbf{contrib}
enters the whole \textbf{wandef} with relative weight \textbf{fac}
1. If there are more than one \textbf{contrib}s the \textbf{fac} of
each \textbf{contrib} can be given. The whole thing is normalized
later on. If the \textbf{wandef} is made of several \textbf{contrib}s
at different sites it is important to correctly define the difference
vectors \textbf{difvec}. Choose a particular point in space as WF
center and express the sites of the \textbf{contrib}s by vectors pointing
from the WF center to the sites \ldots{} these are the \textbf{difvec}s.

If there is more than one site belonging to a crystalographic orbital
(all sites generated by the same Wyckoff position), the user has to
give every WF-projector separately in the moment. In this case it
is important that the user takes care that all the different WF-projectors
belonging to the cryst. orbit are symmetry related. This means setting
up proper \textbf{difvec}, \textbf{x/zaxis}, \textbf{orb}, \textbf{fac}
and so on values.

In our example the orbital used as WF-projector forms a bonding and
an anti-bonding part. There is only one Cu $3d_{x^{2}-y^{2}}$ orbital
in the simple unit cell and hence it can only form one band. However,
this orbital forms sigma-bonding hybrids with the plaquette oxygen
$2p$-orbitals and this allows for and creates ``two'' bands. Moreover,
it is clear from Figure \ref{fig:FPLO-fat-bands-CaCuO} that the anti-bonding
(AB) band is not clearly isolated due to other hybridizations. Hence,
a well defined Wannier function cannot exist. In order to have some
approximation for modeling we create a single WF of AB type. According
to the rules of how WFs are created it is clear that we need to project
away the part of the KS spectrum, which contains the bonding part
of the corresponding band. This is done by specifying an energy window.
All bands outside this window will be ignored. For the majority bands
we have the AB band between -5 and 0 eV and the B band between -8
and -5 eV. For the minority band it is shifted upwards accordingly.
The energy window could be defined as a sharp cutoff. This can lead
to weird effects, if the band character specified through the WF-projector
occurs ONLY outside this window for a particular k-point (such things
can actually happen). This would mean that for one k-point the weight
of the WF-projector in the considered Hilbert-subspace is zero and
this leads to an indefinite problem, when calculating the WFs. So,
it is better to make the window smooth in order to have the interesting
subspace in the main energy window but allowing to sample outside
of it, in case we have a character-run-away scenario. The energy window
is defined by a function being 1 between \textbf{emin} and \textbf{emax}
and falling of as a gaussian, with a width \textbf{de} outside this
window. Of course, the energy window applies to all \textbf{contrib}s
of one \textbf{wandef}. In spin polarized cases these three keywords
take two values one for the majority and one for the minority bands.
The user should also pay attention that the energy window is the same
for all \textbf{wandef}s belonging to the same cryst. orbit! 

Suppose that we defined the above described WF projector onto one
orbital, defining one WF per unit cell and that we did not specify
the energy window. This would lead to a WF whose corresponding band
will be pretty dispersionless situated at an energy between the B
and AB bands in the original band structure. This happens because
the projector $\chi$ makes one single band out of the B and AB part
of the KS-Hilbert space, which essentially forms two bands. In mathematical
terms this will lead to an average of the B and AB KS wavefunctions,
which will have non-bonding character. That should also make clear
what happens, if we extend the energy window more and more towards
the B energies. This will mix in more and more of the bonding KS functions,
which leads to a WF band whose energies get pulled down more and more.
The user can try to play with this to get a feeling.

Now, we have defined all the stuff in the file \texttt{=.wandef} and
we can run FPLO. 
\begin{enumerate}
\item The first thing in the WF module will be that the content of \texttt{=.wandef}
is copied to the output, followed by a more condensed printout of
the wannier parameters. The later output shows all parameters, also
the ones not explicitely set in \texttt{=.wandef} (which then will
have their default values). Then follows a section, in which the symmetry
of the WF-projectors $\chi$ is checked. If this check does not run
through properly there is a mistake in the symmetry relation between
\textbf{wandef}s of $\chi$-s belonging to the same cryst. orbit.
Check all keywords. In the moment the energy window is not checked
for proper symmetry setting. So, if the code runs through, but the
output seems weird, check the energy windows!
\item Now, if the file \texttt{+wancoeff} is found it will be loaded. If
not, the normal FPLO execution will continue. Reading the data is
not the fastest that is why there is an internal loop (see below).
See also Sec. \ref{sub:Saving-memory}.
\item Now, the Wannier bands get processed. As a result the hopping matrix
elements between WFs are printed to the output. This printout is controlled
by user defined restrictions to avoid enormous amounts of data. Example:

\begin{quotation}
spin 1: WF(Cu x2-y2) -\textgreater{} WF(Cu x2-y2) at relative 

T= 0.00000 0.00000 0.00000 hop= -2.093458893797475

T= 0.00000 0.00000 6.04712 hop= -0.061341124468540

T= 0.00000 0.00000 -6.04712 hop= -0.061341124468540

T= -7.29434 0.00000 0.00000 hop= -0.463606088774015

T= 0.00000 -7.29434 0.00000 hop= -0.463606088774014

T= 0.00000 7.29434 0.00000 hop= -0.463606088774014

T= 7.29434 0.00000 0.00000 hop= -0.463606088774015\\

\end{quotation}

The WFs are given by their name, and the symbol '-\textgreater{}'
means that the vectors, which follow ('T= \ldots{}') point from the
WF left of '-\textgreater{}' to the one right of '-\textgreater{}'.
After each vector 'T=' the hopping element in eV is printed. If T=0
the hopping element is the onsite matrix element. The information
written to the output is restricted such that only hoppings with $\left|t\right|>$\textbf{WF\_ham\_threshold}
are considered and only for $\left|T\right|\leq$\textbf{ham\_cutoff}.

\item After this the output on the reciprocial grid is performed. Here it
is important to realize that the Hamiltonian in $k$-space is written,
which corresponds to the Bloch sums of the hoppings Eq. (\ref{eq:WF-hoppings})
\textbf{after} the cut off of hoppings for which $\left|T\right|>$\textbf{ham\_cutoff}
and $\left|t\right|<$\textbf{WF\_ham\_threshold}. The resulting Hamiltonian
is written to \texttt{+hamongrid}. See the \texttt{wannier.f90} soure
code for the order of the data.
\item Now, \texttt{+WF\_coefficients} is written, which contains information
about the contributions of the FPLO orbitals to the WFs Eq. (\ref{eq:WF-from-Phi}).
This file shows only contributions, which are larger than \textbf{WF\_coeff\_threshold}
in \texttt{=.wandef}.
\item At the end of the module the WF are written to disk on a real space
grid. There will be data files \texttt{opendxWF.dx} and \texttt{wfdata...}
and control files \texttt{WF.net\textbar{}cfg}. The later can be use
to visualize the data using \texttt{opendx}:


\texttt{~~~~dx -image WF.net}


Beware that some standard installations of \texttt{opendx }are a bit
buggy. For instance, it might work better if num-lock is turned off.

\item After this FPLO pauses with the message


~~~~~~~CTRL\_C for abort, enter for next trial 


Here, one can type CTRL-C to stop, or one modifies the \textbf{wandef}s
(in another window/editor) and hits enter to re-run the module, with
re-scanning \texttt{=.wandef} but without re-reading \texttt{+wancoeff}
(which is a slow read, since large file).

\end{enumerate}
\begin{figure}[H]
\noindent \begin{centering}
\includegraphics[width=0.8\columnwidth]{pics/wanbands}
\par\end{centering}

\protect\caption{\label{fig:Wannier-fit-for-cacuo2}Wannier fit for CaCuO$_{2}$}


\end{figure}


Figure \ref{fig:Wannier-fit-for-cacuo2} shows the result of the CaCuO3
example. Dark yellow and magenta are the $3d_{x^{2}-y^{2}}$ fat bands
as produced by FPLO. Red shows the full Wannier function transformation
(without cutoffs) from the file +wanband and green shows the cutoff
model WF bands including cutoffs. One can nicely observe how the WF
bands interpolate the hybridization gap between $\Gamma$ and X. and
M and $\Gamma$.


\subsection{Files}

There are three band structure files produced.
\begin{description}
\item [{\texttt{+wanband}}] The Bloch-sums of the WF Eq. (\ref{eq:WF-Bloch-sums})
are obtained by doing the unitary transformation in \textbf{$k$}-space
($U$). The Hamiltonian in WF-Bloch basis in $k$-space Eq. (\ref{eq:WF-Bloch_Ham})
is directly related to the WF-Bloch sums and hence is obtained straight
from the $U$-transformed KS functions and Hamiltonian. The resulting
$k$-dependent Hamiltonian has the dimension of the WF-basis defined
by the user. It can be diagonalized to get the band structure belonging
to the WF model. The result of this is written to \texttt{+wanband}.
If the WFs describe an isolated band complex (one needs as many WFs
as there are KS bands in the band complex) then the \texttt{+wanband}
band structure must coincide with the bandstructure of this band complex
in the full FPLO band structure plot. Deviations are possible, if
there are not enough $k$-points in the FPLO-SCF calculation. This
comes about since a discrete approximation to the $k$-integral in
Eq. (\ref{eq:def_WF}) defines WFs, which have periodic images with
a period of the Born von Karman torus.
\item [{\texttt{+wanbweights}}] This are the fat-bands corresponding to
\texttt{+wanband}. The band weights are determined with respect to
the WF character in the corresponding bands. Let's make that clearer.
In a local orbital basis like FPLO the overlap of the orbitals with
the KS functions determines the fat-bands or how much of a certain
orbital a band is made of. This concept can be applied to the WFs
as well thinking of the WFs as a basis. So, how much of a particular
WF makes up a band? This is given by the information in this file.
Note, that the FPLO fat-bands and WF fat-bands do not coincide. That
is the main reason why we do a WF analysis, see Figure \ref{fig:Comparison-of-FPLO-and-WF-fatbands}
\begin{figure}[H]
\begin{minipage}[c]{0.25\columnwidth}%
\includegraphics[width=1.2\columnwidth]{pics/strnio2}%
\end{minipage}%
\begin{minipage}[t]{0.75\columnwidth}%
\begin{minipage}[t]{1\columnwidth}%
\begin{center}
\includegraphics[width=0.8\columnwidth]{pics/fplofatbands}
\par\end{center}%
\end{minipage}\\
\begin{minipage}[t]{1\columnwidth}%
\begin{center}
\includegraphics[width=0.8\columnwidth]{pics/niowanbands}
\par\end{center}%
\end{minipage}%
\end{minipage}

\protect\caption{\label{fig:Comparison-of-FPLO-and-WF-fatbands}Comparison of FPLO
fat-bands (upper panel) with WF fat-bands (lower panel) for a NiO$_{2}$-chain.
Some WFs are linear combination of several orbitals, such that the
WFs are centered between two orbitals. One can clearly see that in
the WF basis each band has a pure character, i.e. the bands are completely
decoupled by chosing a symmetry adapted basis (WFs).}


\end{figure}
.
\item [{\texttt{+wanbandtb}}] Now, we can go a step further and Fourier
transform the WF-Bloch sums Eq. (\ref{eq:WF-Bloch-sums}) into real
space Eq. (\ref{eq:WF-Bloch-sum-Fouriertrans}). This gives the actual
WFs. In real space these WFs overlap and form hopping integrals Eq.
(\ref{eq:WF-hoppings}). These hopping integrals can be cut off at
some distance or if they are smaller than a threshold. This defines
a modified Hamiltonian in real space. We transform this modified (model)
Hamiltonain back into $k$-space, diagonalize it and get the band
structure in this file. If the cut off and threshold is moderate the
result should equal \texttt{+wanband} otherwise it can differ. Especially,
for non isolated bands we can get non-analyticities in \texttt{+wanband},
which are un-avoidable in some cases. It turns out that the real-space
cutoff reduces these unwanted effects. 
\end{description}
The band structure data files can be visualized using the bandplot-tools
of fplo. Besides these files there are a few more:
\begin{description}
\item [{\texttt{+WFstat\ldots{}}}] contains the absolute values of the
coefficients in Eq. (\ref{eq:WF-from-Phi}) as a function of the distance
of the corresponding orbitals from the Wannier function center. Plot
this file (e.g. \texttt{xfbp ./+WFstat\ldots{}}) to have an idea
of the localization of the WFs.
\item [{\texttt{+T\ldots{}}}] contains the hoppings from a given WF to
the neighbouring WFs corresponding to Eq. (\ref{eq:WF-hoppings})
as a function of the distance between the WF centers.
\end{description}

\subsection{Keywords}

The keywords in \texttt{=.wandef }are explained in the following.
\begin{description}
\item [{doit}] switch the WF module on.
\item [{wandef}] start a definition of a single Wannier function.

\begin{description}
\item [{on/off}] make the corresponding \textbf{wandef} active/inactive.
\item [{emin}] the lower bound of the energy window (two numbers if spin-polarized).
\item [{emax}] the upper bound of the energy window (two numbers if spin-polarized).
\item [{de}] the width of the Gaussian tail above and below \textbf{emax}/\textbf{emin
}defining a smooth energy window (two numbers if spin-polarized).
\item [{delower/deupper}] like \textbf{de}, but for the lower/upper end
of the energy interval. Note, that \textbf{de} will overwrites \textbf{delower}/\textbf{deupper}
if specified after the latter!
\item [{ubands/lbands}] band indices, which specify the maximum and mimimum
band to be contained in the projector (energy window). (two numbers
if spin-polarized) 
\item [{contrib}] add an FPLO orbital to the \textbf{wandef}. There can
be several \textbf{contrib}s in one \textbf{wandef}, one after the
other.

\begin{description}
\item [{site}] the site number of the orbital according to FPLO output.
\item [{name}] an arbitrary name to identify the WF.
\item [{difvec}] the distance of the contributing orbital from th WF center.
The \textbf{difvec} of the first contrib to a \textbf{wandef} together
with the \textbf{site} of the first \textbf{wandef} implicitly define
the WF center. The choice of the WF centers decide whether the WF
is real or not. WF must be defined such that they are real. E.g. a
single orbital/contrib WF should always have \textbf{difvec} 0 0 0. 
\item [{xaxis}] the local x-axis expressed in global coordinates
\item [{zaxis}] the local z-axis expressed in global coordinates
\item [{orb}] the orbital (e.g. 2s+0 or 3d-1). In full-relativistic mode
the orbitals can either be pseudo non-relativistic projections denoted
by ``3d-1 up'' or ``3d-1 dn'' or spherical spinors denoted by
3d3/2-1/2 or 3d5/2+3/2.
\item [{fac}] a weight factor, determining the relative weight of the orbital/contrib
in cases of multiple-contrib wandefs. The weights need not be normalized,
this is taken care of. 
\end{description}
\end{description}
\item [{coefficients\_format}] \label{coefficientsformat}can be \textbf{bin}
or something else. Since version 14.00, the file \texttt{+wancoeff}
can be converted into binary format for faster loading. After the
first FPLO run with a valid \texttt{=.wandef} present, \texttt{+wancoeff}
will have been created. A rerun will start the WF creation process.
If this option is set to \textbf{bin} the data file will be converted
into binary format, if not already done. On any further run the binary
file \texttt{+wancoeffbin} will be read instead of \texttt{+wancoeff}.
This is faster. The user has to take care of deleteing \texttt{+wancoeffbin}
whenever \texttt{+wancoeff} got changed due to settings changes by
the user.
\item [{ham\_cutoff}] restricts output of the real-space WF Hamiltonian
in standard output and \texttt{+T\ldots{}. }Also restricts the matrix
elements used in creating \texttt{+hamongrid}.
\item [{WF\_coeff\_threshold}] restricts the output of coefficients in
\texttt{+WF\_coefficients}.
\item [{WF\_ham\_threshold}] restricts the output of real-space Hamiltonian
in standard output and \texttt{+T\ldots{} }and restricts the hoppings
used in creating \texttt{+wanbandtb }and \texttt{+hamongrid}.
\item [{WF\_write\_coeff\_stats}] can be \textbf{on}/\textbf{off} and triggers
the output of the files \texttt{+WFstat\ldots{}}
\item [{ham\_write\_t\_stats}] can be \textbf{on}/\textbf{off} and triggers
the output of the files \texttt{+T\ldots{}}
\end{description}
For the output of the WFs on the real space grid (visualization) we
have to define a grid
\begin{description}
\item [{WF\_grid\_basis}] can be \textbf{conv}/\textbf{prim}. This defines
the basis $\underline{B}=\left(\begin{array}{c}
\vec{b}_{1}\\
\vec{b}_{2}\\
\vec{b}_{3}
\end{array}\right)$ of the box-like grid. We can use the conventional or the primitive
basis vectors.
\item [{WF\_grid\_directions}] in terms of the basis we can define three
vectors forming the rows of $\underline{V}=\underline{D}\underline{B}$.
These three directions span the grid-box. The input here is the matrix
$\underline{D}$.
\item [{WF\_grid\_subdivision}] subdivide the box along the directions
$\underline{V}$ accordingly.
\item [{WF\_grid\_origin}] put the origin of the box here. If this keyword
is commented out (e.g. \textbf{xWF\_grid\_origin}) the box will be
centered around the WF center.
\end{description}
For the output of the WF Hamiltonian on the $k$-space grid (file
\texttt{+hamongrid}) we have to define a grid
\begin{description}
\item [{k\_grid\_basis}] the basis $\underline{B}$ in reciproical space.
(Note the reciprocial relations between bcc and fcc and the like.)
\item [{k\_grid\_directions}] the directions $\underline{V}=\underline{D}\underline{B}$
defined here by giving $\underline{D}$.
\item [{k\_grid\_subdivision}] subdivisions, see above.
\item [{k\_grid\_incl\_periodic\_points}] can be \textbf{on}/\textbf{off}.
The Hamiltonian in $k$-space is periodic. We can include or exclude
the periodically equivalent points at the boundary of the box.
\end{description}

\subsection{Definition of real spherical harmonics}

We define real spherical harmonics $Y_{lm}$ with magnetic qn. numbers
$m=-l,\ldots,l$ as 
\[
Y_{lm}\left(x,y,z\right)\propto P_{l}^{\left|m\right|}\left(\frac{z}{r}\right)\begin{cases}
\sin\left(\left|m\right|\varphi\right) & \; m<0\\
\cos\left(\left|m\right|\varphi\right) & \; m\geq0
\end{cases}
\]
The sin/cos can be expanded according to the addition theorems, e.g.
$\sin\left(2\varphi\right)=2\sin\varphi\cos\varphi$. Using $\sin\varphi\propto y$
and $\cos\varphi\propto x$ we get $\sin\left(2\varphi\right)\propto xy$.
Additionally we need $P_{l}^{\left|m\right|}\left(\frac{z}{r}\right)\propto\mbox{polynomial of degree }l-\underline{m}\mbox{ in }z$.
Thus,

\[
Y_{2,-2}\propto P_{2}^{2}\left(\frac{z}{r}\right)\sin\left(2\varphi\right)\propto xy
\]
\[
Y_{2,-1}\propto P_{2}^{1}\left(\frac{z}{r}\right)\sin\left(\varphi\right)\propto zy
\]
\[
Y_{2,0}\propto P_{2}^{0}\left(\frac{z}{r}\right)\cos\left(0\varphi\right)\propto z^{2}
\]
\[
Y_{2,1}\propto P_{2}^{1}\left(\frac{z}{r}\right)\cos\left(\varphi\right)\propto zx
\]
\[
Y_{2,2}\propto P_{2}^{2}\left(\frac{z}{r}\right)\cos\left(2\varphi\right)\propto x^{2}-y^{2}
\]
For more specific information on the polynomial in $z$ one has to
look up the associated Legendre polynomials $P_{l}^{\left|m\right|}\left(z\right)$.


\subsection{Problems }

The number of $\mathbf{k}$-points used in the SCF calculation influences
the Wannier function quality. If the $\mathbf{k}$-mesh is not fine
enough, the band structure determined by diagonalizing the WF Bloch
Hamiltonain Eq. (\ref{eq:WF-Bloch_Ham}) (\texttt{+wanband}) will
not coincide with the corresponding bands of the full FPLO band structure
(\texttt{+band}). (Coincidence can of course only happen anyway, if
the WFs describe an isolated band complex.) 

If the real space Hamiltonian cutoff is smaller than the extend of
the WFs the fitted Wannier bandstructure (\texttt{+wanbandtb}) will
not conincide with the full WF band structure. This may also be due
to a bad Wannier definition, which is not localizing the WFs sufficiently
(check energy window, use molecular orbitals (several \textbf{contrib}s
per \textbf{wandef}) instead of simple orbitals).

There might be unnatural spikes at certain $k$-points in the WF band
structure. This usually means that the energy window is to narrow
and that the character of $\chi$ is only large outside the window
at the corresponding k-points.

If slabs or chains are calculated, there is an artificial periodicity
in the irrelevant directions. If not enough SCF $\mathbf{k}$-points
are used in these directions, we get influences of the artificial
periodic replica. Hence, either the vacuum spacing is very large,
which is detrimental for other parts of the code, or one has to use
sufficiently many points in the irreevant directions. (Only some experimenting
can tell.)


\subsection{Saving memory (and time)\label{sub:Saving-memory}}

In order to reduce the size of +wancoeff the user can set the FEDIT
option ``restrict bands to window'' in the bandplot submenu. Then
lower and upper energy bounds in the same submenu defines the bands,
which will be written to \texttt{+wancoeff} for all data contained
in this file.

Use the \htmlref{\textbf{coefficients\_format}}{coefficientsformat}
option in \texttt{=.wandef} in order to convert \texttt{+wancoeff}
to binary format. 


\section{Examples}


\subsection{Hexagonal: $sp^{2}$, Graphene, MgB$_{2}$}

Here we analyse the Wannier function choices for a hexagonal lattice
with essentially two atoms per unit cell. The main issue is to understand
the symmetry considerations.

Have a look at the accompanying example directories for MgB$_{2}$.

\noindent \begin{center}
\begin{figure}[H]
\noindent \begin{centering}
\includegraphics[width=0.6\columnwidth]{pics/wanhex1}
\par\end{centering}

\protect\caption{\label{fig:real-and-reciprocal}Real and reciprocal cell of graphene,
MgB$_{2}$}
\end{figure}

\par\end{center}

Figure \ref{fig:real-and-reciprocal} shows the basic lattice structure.
We have a unit cell given by
\[
\boldsymbol{A}=\left(\begin{array}{c}
\boldsymbol{a}_{1}^{T}\\
\boldsymbol{a}_{2}^{T}\\
\boldsymbol{a}_{3}^{T}
\end{array}\right)=\underline{A}\boldsymbol{E},\quad\boldsymbol{E}=\left(\begin{array}{c}
\boldsymbol{e}_{x}^{T}\\
\boldsymbol{e}_{y}^{T}\\
\boldsymbol{e}_{z}^{T}
\end{array}\right),\quad\underline{A}=\left(\begin{array}{ccc}
a_{H} & 0 & 0\\
0 & a_{H} & 0\\
0 & 0 & c_{H}
\end{array}\right)\left(\begin{array}{ccc}
\frac{\sqrt{3}}{2} & -\frac{1}{2} & 0\\
0 & 1 & 0\\
0 & 0 & 1
\end{array}\right)
\]
and two sites $\boldsymbol{s}=\underline{s}^{T}\underline{A}$, $\underline{s}_{1}=\left(\frac{1}{3}\frac{2}{3}0\right)$,
$\underline{s}_{2}=\left(\frac{2}{3}\frac{1}{3}0\right)$. The symmorphic
space group is $191$ with a maximal point group $D_{6h}$ at $\left(000\right)$
with generators $C_{6}\left(z\right)$, $C_{2}\left(0\right)$ and
$I$ (indicated in Figure \ref{fig:real-and-reciprocal}). Clearly
$C_{6}\boldsymbol{s}_{1}=\boldsymbol{s}_{2}$. The local pointgroup
at the two sites is $D_{3h}$ with generators $C_{3}$, $C_{2}\left(0\right)$
and $\overline{6}=IC_{6}$. For the sake of definitness let's define
operations. We understand the operation $g$ to actively transform
the basis vectors $\boldsymbol{e}_{x,y,z}$ like 
\[
g\boldsymbol{E}^{T}=\boldsymbol{E}^{T}R\left(g\right)
\]
where $R\left(g\right)$ is the 3D representation matrix. For example
\[
g=C_{6}=\left(\begin{array}{ccc}
\frac{1}{2} & -\frac{\sqrt{3}}{2} & 0\\
\frac{\sqrt{3}}{2} & \frac{1}{2} & 0\\
0 & 0 & 1
\end{array}\right)
\]
The coordinates $\underline{r}$ of a vector $\boldsymbol{r}=\boldsymbol{E}^{T}\underline{r}$
then transform with the transposed matrix: $\underline{r}^{\prime T}=\underline{r}^{T}R\left(g\right)^{T}$
(proof: $\boldsymbol{E}^{T}\underline{r}^{\prime}=g\boldsymbol{r}=\boldsymbol{E}^{T}R\left(g\right)\underline{r}$).
A set of functions of the coordinates (orbitals) transforms as 
\[
gf_{m}\left(\underline{r}\right)=f_{m}\left(g^{-1}\underline{r}\right)=f_{m}\left(\underline{r}^{T}R\left(g\right)\right)=\sum_{n}f_{n}\left(\underline{r}\right)D_{nm}^{f}\left(g\right)
\]
where $D_{nm}^{f}\left(g\right)$ are the representation matrices
of the group in the space spanned by $f_{m}$. In other words the
orbitals $gf_{m}\left(\underline{r}\right)$ are also actively transformes
like the cartesian basis vectors $\boldsymbol{e}_{x,y,z}$.

In our examples we have a light $p$ element siting at the two sites
of lower than maximal symmetry. The essential basis atoms are one
$s$ and three $p$ orbitals. Figure \ref{fig:Bandweights-of-boron}
shows the band structure of MgB$_{2}$, where the boron bands are
highlighted. (The Mg bands do not play a big role in this energy window.)
\begin{figure}[H]
\noindent \begin{centering}
\includegraphics[width=0.6\columnwidth]{pics/MgB2/bwfplo}
\par\end{centering}

\protect\caption{\label{fig:Bandweights-of-boron}Bandweights of boron in MgB$_{2}$}


\end{figure}


The irreducible representations of the site symmetry $D_{3h}$ make
$p_{x,y}$ an $E^{\prime\prime}$ doublet and $p_{z}$ a $A_{2}^{\prime\prime}$
singlet, while the $s$-orbital represents an $A_{1}^{\prime}$ singlet.
In total we have $8$ oribitals and hence bands. Due to the symmetry
the $p_{z}$ bands largely decouple from the other orbitals along
high symmetry lines in the $k_{z}=0$ and $k_{z}=\frac{\pi}{c_{H}}$
planes. At general points they off course couple. Furthermore, these
bands cross the Fermi level, while the $s$ and $p_{x,y}$ bands form
bonding and anti-bonding band complexes of three bands each separated
by a gap. The bonding band complex has exactly three bands, while
the AB complex shows a high degree of band entangling wiht other bands.
The two $p_{z}$ bands cannot be separated into two distinct bands
because of their Fermi surface and symmetry. However the planar-orbital
bands can be represented by different kinds of Wannier functions.
The simplest way of defining WFs is to use the atomic $s$, $p_{x,y}$
orbitals sitting at the two sites. This more or less reproduces the
FPLO band characters. Another alternative is to try to find WFs for
the bonding bands (anti-bonding) bands only). There are essentially
three bands in each band complex and it is imediately clear that atom
centered functions cannot fulfill the crystall symmetry due to the
number of sites (two). The bands must be combinations of orbitals
from both sites and hence an odd (3) number of orbitals from an even
number of sites (2) cannot be symmetric. We discussed the irreducible
representations of the orbitals in the site symmetry above. We have
a singlet $\left(s\right)$ and a doublet ($p_{x,y}$) (as is reflected
in the degeneracies at the $\Gamma$-point) but all three orbitals
are mixed in the bonding bands (no further decoupling). 
\begin{figure}[H]
\noindent \begin{centering}
\includegraphics[width=0.6\columnwidth]{pics/MgB2/sp2-orbitals}
\par\end{centering}

\protect\caption{\label{fig:sp2-hybrid-basis}$sp^{2}$-hybrid basis. The left panel
shows one $sp^{2}$-orbital at both sites (filled) and the remaining
two $120^{\circ}/240^{\circ}$-rotated orbitals at $s_{1}$ (open).
The right panel show a bonding bond-centered combination (full) and
two symmetry related combinations (open).}
\end{figure}
 Unfortunately, there is hardly a simple linear combination of say
the $p_{x,y}$ orbitals from both sites, which will fulfill the full
crystal symmetry. Hence, we have to resort to reducible representations
aiming at constructing as highly symmetric orbitals as possible. The
well known answer is of course $sp^{2}$ hybrids. We introduce the
following orbitals at each site 
\[
\left(\begin{array}{c}
\Phi_{1}\\
\Phi_{2}\\
\Phi_{3}
\end{array}\right)=\frac{1}{\sqrt{3}}\left(\begin{array}{c}
s+\sqrt{2}p_{x}\\
s+\sqrt{2}C_{3}p_{x}\\
s+\sqrt{2}C_{3}^{2}p_{x}
\end{array}\right)
\]
which are depicted in Figure \ref{fig:sp2-hybrid-basis}. 

Obviously, $C_{3}p_{x}$ is not parallel to $p_{x}$ and hence contains
some $p_{y}$. On the other hand $\left(1+C_{3}+C_{3}^{2}\right)=3P_{\parallel}$
where $P_{\parallel}$ is a projector onto the rotation axis. For
in-plane objects $P_{\parallel}p_{x,y}=0$. This consideration shows
that the $p$-parts of $\Phi_{1,2,3}$ are linearly dependent, however,
with the $s$-admixture they span the same space as $s$ and $p_{x,y}$.
The factors are choosen such that $\left|\Phi_{i}\right|^{2}=1$ and
that $\sum_{i}\left|\Phi_{i}\right|^{2}$ contains one $\left|s\right|^{2}$
density and two $\left|p\right|^{2}$ densities. 

Now, we can put three $sp^{2}$ hybrids onto each site (rotated by
$180^{\circ}$ at the second site, which is not really necessary)
and clearly they posses maximum crystal symmetry. (Figure \ref{fig:sp2-hybrid-basis},
left panel). The symmetry relations are
\[
C_{6}\Phi_{s_{2},1}=\Phi_{s_{1},3},\ldots
\]
\[
C_{3}\Phi_{s,1}=\Phi_{s,2},\ldots
\]
\[
C_{2}\left(0\right)\Phi_{s_{2},1}=\Phi_{s_{2},1},\quad C_{2}\left(0\right)\Phi_{s_{2},2}=\Phi_{s_{2},3},\ldots
\]
\[
I\Phi_{s_{2},1}=\Phi_{s_{1},2},\ldots
\]
The resulting Wannier function fat bands are shown in Figure \ref{fig:3-orbital--Wannier-function}.

\begin{figure}
\noindent \begin{centering}
\includegraphics[width=0.8\columnwidth]{pics/MgB2/wanbweights_6sp2}\includegraphics[width=0.2\columnwidth]{pics/MgB2/sp21_at_2.epsf}
\par\end{centering}

\protect\caption{\label{fig:3-orbital--Wannier-function}Wannier function fit with
3 $sp^{2}$-orbital per site.}


\end{figure}
The advantage of this high symmetry is that we now can create linear
combinations of one $sp^{2}$ from each site at the bond centers and
by construction this generates symmetry related additional two similar
combinations (Figure \ref{fig:sp2-hybrid-basis}, right panel). This
were not possible if we used the irreducible $s$ and $p_{xy}$ orbitals.
In that way one can separate the bonding from the anti-bonding orbitals.
This is an examples for the fact that in many covalent $p$-electron
systems the bonding/anti-bonding bands must be described by bond-centered
Wannier functions. (Carbon chain: $sp^{1}$, diamond $sp^{3}$). Once
we have separated the bonding Wannier functions, we can of course
produce either a separate fit for the bonding or anti-bonding bands
respectively (Figure \ref{fig:Bonding/anti-bonding-bond-center})
or a fit for bonding and anti-bonding bands together (however with
clear energy separation of the orbital).
\begin{figure}[H]
\noindent \begin{centering}
\includegraphics[width=0.8\columnwidth]{pics/MgB2/wanbweights_sp2_B}\includegraphics[width=0.2\columnwidth]{pics/MgB2/image_B.epsf}
\par\end{centering}

\noindent \begin{centering}
\includegraphics[width=0.8\columnwidth]{pics/MgB2/wanbweights_sp2_AB}\includegraphics[width=0.2\columnwidth]{pics/MgB2/image_AB.epsf}
\par\end{centering}

\protect\caption{\label{fig:Bonding/anti-bonding-bond-center}Bonding/anti-bonding
bond centered Wannier function fits.}
\end{figure}


One last consideration concerns the choice of energy windows. For
the entangled anti-bonding bands to be fitted better we put a narrow
energy window onto the occupied bonding bands but allow for a large
Gaussian tail at the upper side of the enrgy window. This pulls the
majority weight from the well separated bonding bands and magically
adds the missing anti-bonding weights. If we extend the main energy
window to encompass the anti-bonding bands it would pull in weights
from higher lying bands, which also show $s$ and $p_{x,y}$ character,
which in turn pulls the fitted bands up in energy.


\end{document}
