abinit/doc/guide/elphon_manual.tex

\documentclass[aps,preprint]{revtex4-1}
%\documentclass{article}

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

\unitlength 1cm

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

\title {Abinit electron phonon interaction calculations for geniuses}
\author{Matthieu Verstraete, Bin Xu, and Momar Diakhate}
\affiliation{Unit\'e Nanomat, Dept of Physics, Universit\'e de Li\`ege, 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
authors suppose 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 electron-phonon response
function (eph) 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 second order 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 k-points. K-point 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
k-points on the FS, the k-point mesh must be exceptionaly fine. A Monkhorst-Pack
16x16x16 grid appears well converged for simple FCC systems like Al or Pb, but
in many cases 24x24x24 or 36x36x36 grids are not rare. 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 as essential as that of the
k-point one, since both k and k+q must be on the FS.
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 k-point
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. Similarly a $\Gamma$ centered grid for the electron
k-points is better - a shifted grid appears to work, but the related
convergence properties are not ensured.

Here comes the second limitation and warning: by default abinit only does the
minimal number of irreducible perturbations (atoms, directions and qpoints),
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.

We will calculate the matrix elements (potentially on a different, denser,
k-point grid) in a second step. These are contained in files suffixed GKK. 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, 1 atom, and 3
directions). In the abinit input file, the datasets are divided into
\begin{itemize}
\item the ground state with a normal number of k-points
\item the phonon perturbations, with the same grid of k-points, sufficient to
converge good phonon frequencies. From these we will recover the perturbed
electronic density and the DDB files (second derivatives of the total energy)
\item the calculation of the DDK perturbation, which will be used to calculate
the Fermi velocities, when they are needed for transport calculations
\item a non self-consistent calculation of the GS wavefunctions on a dense grid
of k
\item a non self-consistent calculation of the matrix elements of the EPI, on
the dense grid. The hamiltonian is reconstructed using the density converged
above. Here the prtgkk flag is set to 1 to output the GKK files, and only 1
step is necessary, as there is nothing to converge (we care only about the
matrix elements, not the perturbed wave functions)
\item a non self-consistent calculation of the DDK matrix elements (idem for transport)
\end{itemize}

Important Nota Bene: the last 3 steps are new to abinit version 7.6. The prtgkk flag now
disables some symmetry reduction of k-points, and will make your self
consistent phonon calculations much slower (unnecessarily). The new scheme
decouples the calculation of the GKK and the convergence of the phonons, which
do not need a super dense k-point grid. The latter are now complete (full
k-point grid with kptopt 3), and calculated for all of the perturbations. This
was important for certain degenerate electronic bands, whose correspondence
between equivalent k-points can contain a random rotation in the degenerate
subspace.


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

We will merge all the DDB files into one Al total ddb, as in the
standard tutorials, with mrgddb.
A small utility has been added to the package to merge the necessary bits of
the matrix element 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 GKK files into the present directory, with
corresponding names for each qpoint, modify the names in teph\_3.in, and run as
\begin{verbatim}
mrggkk < teph_3.in
\end{verbatim}

As detailed in the example input file, the first line gives the name of the
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 (usually 0 now, but these can also be used to extract the matrix elements),
and 24 24 (the difference between the last two numbers allows you to re-merge GKK
files with several perturbations per file, but this won't be needed here). 
Follow the names of all the 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
teph\_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 k-point 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 k-point 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 k-point and qpoint
grids), anaddb does a Fourier interpolation of the matrix elements, to obtain
them on a fine grid of qpoints, identical to the k-point 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
k-points, which can then be integrated over the Fermi Surface.


%%
%%  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}