\documentclass[aps,preprint]{revtex4}

\usepackage[dvips]{graphicx}
\usepackage{textcomp}
\graphicspath{figs}

\unitlength 1cm

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}

\title {Abinit electron phonon interaction calculations for geniuses}
\author{Matthieu Verstraete}
\affiliation{Unit\'e de Physico-Chimie et de Physique des Mat\'eriaux (PCPM), and
CERMIN, Universit\'e Catholique de Louvain, 1 Croix du Sud, B-1348
Louvain-la-Neuve, Belgium}

\maketitle


\section{INTRODUCTION}
This document is designed as a rudimentary tutorial to learn to use the
electron phonon capabilities implemented in the ABINIT software package. The
author supposes that the reader has knowledge of the theory involved in basic
Density Functional Theory (DFT)\cite{ABINIT_short, payne_1992_review_dft,
hohenberg_1964_DFT, kohn_1965_DFT_LDA}, Density Functional Perturbation theory
(DFPT)\cite{baroni_1987_elast_const_lin_resp, baroni_2001_phonon_review,
gonze_1989_DFT_nonlinear_response, gonze_1997_higher_order_DFTPT,
gonze_1995_perturbation_variational}, as well as that involved in the standard
descrption of superconductivity and the electron-phonon interaction
(EPI)\cite{allen_1982_superconducting_tc}, including some of its
implementations in DFT\cite{dacorogna_1985_EPI, lam_1986_EPI,
savrasov_1996_EPI_implementation}. Some references to the latter may be added
in the text, but this is in no event an introduction to any of these fields.
It is at best a rough description of the usage of ABINIT for standard
run-of-the-mill calculations on simple systems. It follows the elphon response
function tutorial in ABINIT, and explains some of the operations and concepts.

Other standard abbreviations include Brillouin Zone (BZ), and Fermi Surface
(FS).

\section{CALCULATING THE RAW ELECTRON-PHONON MATRIX ELEMENTS}

The actual calculation of the EPI matrix elements is performed in the course of
the secondorder-perturbation-theory run, routinely used to calculate phonon
frequencies, Born effective charges, phonon polarization vectors etc... The
matrix elements are by-products of the Sternheimer equation used to calculate
the first-order perturbed wavefunctions $\psi^{(1)}$. In this equation, the
ground state hamiltonian $H^{(0)}$ is applied to $\psi^{(1)}$ and the first
order hamiltonian $H^{(1))}$ is applied to $\psi^{(0)}$. The calculation of the
scalar product of $\langle \psi^{(0)}_{k,n} \mid$ with $H^{(1)} \mid
\psi^{(0)}_{k+q,n'} \rangle$ yields the matrix element we are looking for. $H^{(1)}$
denotes only the perturbation of the hamiltonian for a given q vector with
respect to one atomic position $R_\tau$ in one reduced direction $\alpha$, i.e.
$H^{q \tau \alpha}$.

The first step is a ground-state calculation of the density and wavefunctions
on a sufficiently fine grid of kpoints. Kpoint convergence is an important
issue in the EPI as only electrons on or near the Fermi Surface will
contribute. As the EPI is in general very small except at a few points in the
BZ, and the final integration contains a double delta function selecting
kpoints on the FS, the kpoint mesh must be exceptionaly fine. A Monkhorst-Pack
16x16x16 grid appears well converged for simple FCC systems like Al or Pb. See
the appendices for sample input files at each step of the calculation (GS,
response function, mrggkk and anaddb) for FCC aluminium using the HGH
pseudopotential Al.hgh provided with abinit.


With the ground state data in hand, following the standard ABINIT tutorial, do
a response-function calculation of the phonons on a fine enough qpoint grid in
the BZ. The density of the qpoint mesh is less essential than that of the
kpoint one, as we will see, and a grid converged for the calculation of phonon
frequencies appears to be enough. Here we use a 4$\times$4$\times$4 grid of qpoints.

First limitation and warning: the qpoint mesh must be a sub-mesh of the kpoint
mesh, and must contain $\Gamma$ !: qpt 2$\times$2$\times$2 with kpt
10$\times$10$\times$10 but not qpt 3$\times$3$\times$3 with kpt
10$\times$10$\times$10. This constraint could be relaxed, but a bit more work
would be needed.

Here comes the second limitation and warning: by default abinit only does the
minimal number of irreducible perturbations (atoms, directions and qpoints) it
needs to, and presumes the rest can be re-generated by symmetry. The
completion of the EPI matrix elements over qpoints is implemented in anaddb,
but NOT the completion over atoms or reduced directions. You must therefore
force abinit to do all calculations for all atoms in all directions, e.g. by
specifiying separately

\begin{verbatim}
rfatpol1 1 1
rfdir1 1 0 0

rfatpol2 1 1
rfdir2 0 1 0

rfatpol3 1 1
rfdir3 0 0 1

rfatpol4 2 2
rfdir4 1 0 0
 ...
\end{verbatim}

for each irreducible qpoint. It is not yet clear whether this restriction will
ever be lifted. Each response-function calculation will produce a 1WF file
containing the first order wavefunctions and our precious first-order matrix
elements. At the end of this section you should have $n_{qpt} \times 3 \times
n_{atom}$ of these files. In our example case that makes 24 (8 irred qpoints).
Merge all the DDB files into one Al total ddb, as in the standard tutorials,
with mrgddb.


\section{EXTRACTING AND MERGING THE MATRIX ELEMENTS INTO ONE FILE}

A small utility has been added to the package to merge the necessary bits of
the 1WF files together into a single file, containing the needed header data
and the raw first-order matrix elements, abusively called a GKK file. The
utility is called mrggkk in reference to mrgddb which does the same thing for
DDB files. Copy all the 1WF files into the present directory, with
corresponding names for each qpoint, modify the names in telphon\_3.in, and run
as
\begin{verbatim}
mrggkk < telphon_3.in
\end{verbatim}

As detailed in the example input file, the first line gives the name of the
final output file, the second should be kept to 0 (binary output). The third
line gives the name of the GS wavefunction file.  The fourth line contains the
number of 1WF files, and 0 0 (the last two numbers allow you to re-merge GKK
files but this won't be needed). Follow the names of all the 1WF files,
ordered by qpt.

This is the third important important limitation and warning: the code
presumes the perturbations for one qpoint are grouped, and throws them out
once it finds a new qpoint. If all the perturbations were not present, anaddb
will exit when it realises information is missing.

\section{RUNNING ANADDB}

The main external parameter for the moment is the elphflag variable in the
telphon\_4.in file, which should be set to 1. Future variables will include the
$\mu^{\star}$ parameter used in the determination of T$_c$.

\subsection{Calculation of $g_{kn,k'n'}$}
The matrix elements which come out of ABINIT are not exactly those used in
electron-phonon theory, because ABINIT works with simple perturbations of one
atom in one crystalline direction. The normal definition of the ``GKK" matrix
element is:

\begin{eqnarray}
g^{q,j}_{k',n';k,n} &=& \sqrt{\frac{1}{2 M_\tau \omega_{q,j}}} \langle \psi_{k',n'} \mid H^{q,j}_{k',k} \mid \psi_{k,n}\rangle \\
&=& \sqrt{\frac{1}{2 M_\tau \omega_{q,j}}} \sum_{\tau, \alpha} {e}^{q,j}_{\tau, \alpha} \cdot \langle \psi_{k',n'} \mid H^{\tau, \alpha}_{k',k} \mid \psi_{k,n}\rangle
\end{eqnarray}

where $\psi_{k,n}$ is the wavefunction at kpoint k and band n, $M_\tau$ is the
mass of the atomic species, $\omega_{q,j}$ and ${e}^{q,j}_{\tau, \alpha}$
are the frequency and eigenvector of the phonon mode j, and $H^{\tau,
\alpha}_{k',k}$ is the first order perturbing hamiltonian, moving atom $\tau$
in direction $\alpha$, which is actually applied in ABINIT.

We see that the $g^{q,j}_{k',n';k,n}$ is phonon-mode specific, instead of atom and
direction. As its calculation requires all $\tau$, $\alpha$ perturbations for a given q,
the first step is to read in all the matrix elements, calculate the phonon
frequencies, and perform the scalar products with ${e}^{q,j}$.

\subsection{Integration over the Fermi Surface and isotropic constants}
At this point, we can calculate all isotropic constants (T$_c$, $\lambda$,
$\omega_{log}$) and the FS-averaged phonon linewidths. The
$g^{q,j}_{k',n';k,n}$ are summed first over $n$ and $n'$, using a weighting factor
$exp -((\epsilon_{k,n} - \epsilon_F )/\sigma)^2$ for each. The $\sigma$ is an input variable.

Further, the $g^{q,j}_{k';k}$ are integrated over k, to give $g_{q,j}$. The
phonon linewidth is just $\gamma_{q,j} = 2 \pi \omega_{q,j} \cdot g^{q,j}$ and
is output over a path in reciprocal space corresponding to the FCC special
points $\Gamma$ -X-$\Gamma$ -L-X-W -L (corresponding to the qpath and nqpath
input variables). The values of the linewidths on the path are output to a file
appended LWD.

The $g^{q,j}$ are interpolated wrt q, as described in the next section, on a much
finer grid (the kpoint grid), and integrated over q and j to give the
Eliashberg spectral function:

\begin{equation}
\alpha^2 F(\omega) = \frac{1}{N(\epsilon_F)} \sum_{q,j} g^{q,j} \delta(\omega - \omega_{q,j})
\end{equation}
output to a file appended A2F.

\subsection{Interpolation with respect to q}
Once the actual GKK have been calculated (on the given kpoint and qpoint
grids), anaddb does a Fourier interpolation of the matrix elements, to obtain
them on a fine grid of qpoints, identical to the kpoint grid. In this way, $k$,
$q$, and $k' = k + q$ all span the whole grid of points in the BZ. The
$g^{q,j}_{k',n';k,n}$ are FTed to $g^{R,j}_{k',n';k,n}$ on a set of points R in real space
chosen in the same way as for phonon interpolation in anaddb, with weights to
account for their belonging (or not) to the surface of the Wigner-Seitz cell in
real space. Then the $g^{R,j}_{k',n';k,n}$ are FTed back for all q on the fine mesh of
kpoints. For the moment these are not used, but work goes on.

\section{DISCLAIMER}

This piece of code is not waranteed suitable for any purpose save my own. As
of today September 9 2004 it has only been tested of FCC Al and Pb: no
multi-component or noncubic systems. If anyone has successes or bugs to
report, please do so - I will appreciate it.

%%
%%  Bibliography
%%
\begin{thebibliography}{13}
\expandafter\ifx\csname natexlab\endcsname\relax\def\natexlab#1{#1}\fi
\expandafter\ifx\csname bibnamefont\endcsname\relax
  \def\bibnamefont#1{#1}\fi
\expandafter\ifx\csname bibfnamefont\endcsname\relax
  \def\bibfnamefont#1{#1}\fi
\expandafter\ifx\csname citenamefont\endcsname\relax
  \def\citenamefont#1{#1}\fi
\expandafter\ifx\csname url\endcsname\relax
  \def\url#1{\texttt{#1}}\fi
\expandafter\ifx\csname urlprefix\endcsname\relax\def\urlprefix{URL }\fi
\providecommand{\bibinfo}[2]{#2}
\providecommand{\eprint}[2][]{\url{#2}}

\bibitem[{\citenamefont{Gonze et~al.}(2002)\citenamefont{Gonze, Beuken,
  Caracas, Detraux, Fuchs, Rignanese, Sindic, Verstraete, Zerah, Jollet
  et~al.}}]{ABINIT_short}
\bibinfo{author}{\bibfnamefont{X.}~\bibnamefont{Gonze}},
  \bibinfo{author}{\bibfnamefont{J.-M.} \bibnamefont{Beuken}},
  \bibinfo{author}{\bibfnamefont{R.}~\bibnamefont{Caracas}},
  \bibinfo{author}{\bibfnamefont{F.}~\bibnamefont{Detraux}},
  \bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Fuchs}},
  \bibinfo{author}{\bibfnamefont{G.-M.} \bibnamefont{Rignanese}},
  \bibinfo{author}{\bibfnamefont{L.}~\bibnamefont{Sindic}},
  \bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Verstraete}},
  \bibinfo{author}{\bibfnamefont{G.}~\bibnamefont{Zerah}},
  \bibinfo{author}{\bibfnamefont{F.}~\bibnamefont{Jollet}},
  \bibnamefont{et~al.}, \bibinfo{journal}{Comp. Mat. Sci.}
  \textbf{\bibinfo{volume}{25}}, \bibinfo{pages}{478} (\bibinfo{year}{2002}).

\bibitem[{\citenamefont{Payne et~al.}(1992)\citenamefont{Payne, Teter, Allan,
  Arias, and Joannopoulos}}]{payne_1992_review_dft}
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Payne}},
  \bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Teter}},
  \bibinfo{author}{\bibfnamefont{D.}~\bibnamefont{Allan}},
  \bibinfo{author}{\bibfnamefont{T.}~\bibnamefont{Arias}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Joannopoulos}},
  \bibinfo{journal}{Rev. Mod. Phys.} \textbf{\bibinfo{volume}{64}},
  \bibinfo{pages}{1045} (\bibinfo{year}{1992}).

\bibitem[{\citenamefont{Hohenberg and Kohn}(1964)}]{hohenberg_1964_DFT}
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Hohenberg}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{W.}~\bibnamefont{Kohn}},
  \bibinfo{journal}{Phys. Rev.} \textbf{\bibinfo{volume}{136}},
  \bibinfo{pages}{864} (\bibinfo{year}{1964}).

\bibitem[{\citenamefont{Kohn and Sham}(1965)}]{kohn_1965_DFT_LDA}
\bibinfo{author}{\bibfnamefont{W.}~\bibnamefont{Kohn}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{L.}~\bibnamefont{Sham}},
  \bibinfo{journal}{Phys. Rev.} \textbf{\bibinfo{volume}{140}},
  \bibinfo{pages}{A 1133 } (\bibinfo{year}{1965}).

\bibitem[{\citenamefont{Baroni et~al.}(1987)\citenamefont{Baroni, Giannozzi,
  and Testa}}]{baroni_1987_elast_const_lin_resp}
\bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Baroni}},
  \bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Giannozzi}},
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Testa}},
  \bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{59}},
  \bibinfo{pages}{2662} (\bibinfo{year}{1987}).

\bibitem[{\citenamefont{Baroni et~al.}(2001)\citenamefont{Baroni, de~Gironcoli,
  Dal~Corso, and Giannozzi}}]{baroni_2001_phonon_review}
\bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Baroni}},
  \bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{de~Gironcoli}},
  \bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Dal~Corso}},
  \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Giannozzi}},
  \bibinfo{journal}{Rev. Mod. Phys.} \textbf{\bibinfo{volume}{73}},
  \bibinfo{pages}{515} (\bibinfo{year}{2001}).

\bibitem[{\citenamefont{Gonze and
  Vigneron}(1989)}]{gonze_1989_DFT_nonlinear_response}
\bibinfo{author}{\bibfnamefont{X.}~\bibnamefont{Gonze}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{J.-P.} \bibnamefont{Vigneron}},
  \bibinfo{journal}{Phys. Rev. B} \textbf{\bibinfo{volume}{39}},
  \bibinfo{pages}{13120} (\bibinfo{year}{1989}).

\bibitem[{\citenamefont{Gonze}(1997)}]{gonze_1997_higher_order_DFTPT}
\bibinfo{author}{\bibfnamefont{X.}~\bibnamefont{Gonze}},
  \bibinfo{journal}{Phys. Rev. B} \textbf{\bibinfo{volume}{55}},
  \bibinfo{pages}{10337} (\bibinfo{year}{1997}).

\bibitem[{\citenamefont{Gonze}(1995)}]{gonze_1995_perturbation_variational}
\bibinfo{author}{\bibfnamefont{X.}~\bibnamefont{Gonze}},
  \bibinfo{journal}{Phys. Rev. B} \textbf{\bibinfo{volume}{52}},
  \bibinfo{pages}{1086} (\bibinfo{year}{1995}).

\bibitem[{\citenamefont{Allen and
  Mitrovi\'{c}}(1982)}]{allen_1982_superconducting_tc}
\bibinfo{author}{\bibfnamefont{P.~B.} \bibnamefont{Allen}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{B.}~\bibnamefont{Mitrovi\'{c}}},
  \emph{\bibinfo{title}{Theory of Superconducting T$_{c}$}},
  vol.~\bibinfo{volume}{37} of \emph{\bibinfo{series}{Solid State Phys.}}
  (\bibinfo{publisher}{Academic Press}, \bibinfo{address}{New York},
  \bibinfo{year}{1982}).

\bibitem[{\citenamefont{Dacorogna et~al.}(1985)\citenamefont{Dacorogna, Cohen,
  and Lam}}]{dacorogna_1985_EPI}
\bibinfo{author}{\bibfnamefont{M.~M.} \bibnamefont{Dacorogna}},
  \bibinfo{author}{\bibfnamefont{M.~L.} \bibnamefont{Cohen}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{P.~K.} \bibnamefont{Lam}},
  \bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{55}},
  \bibinfo{pages}{837} (\bibinfo{year}{1985}).

\bibitem[{\citenamefont{Lam et~al.}(1986)\citenamefont{Lam, Dacorogna, and
  Cohen}}]{lam_1986_EPI}
\bibinfo{author}{\bibfnamefont{P.~K.} \bibnamefont{Lam}},
  \bibinfo{author}{\bibfnamefont{M.~M.} \bibnamefont{Dacorogna}},
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{M.~L.} \bibnamefont{Cohen}},
  \bibinfo{journal}{Phys. Rev. B} \textbf{\bibinfo{volume}{34}},
  \bibinfo{pages}{5065} (\bibinfo{year}{1986}).

\bibitem[{\citenamefont{Savrasov and
  Savrasov}(1996)}]{savrasov_1996_EPI_implementation}
\bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Savrasov}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{D.}~\bibnamefont{Savrasov}},
  \bibinfo{journal}{Phys. Rev. B} \textbf{\bibinfo{volume}{54}},
  \bibinfo{pages}{16487} (\bibinfo{year}{1996}).

\end{thebibliography}

\end{document}
