1\documentclass[12pt,a4paper]{article} 2\def\version{6.6} 3\def\qe{{\sc Quantum ESPRESSO}} 4 5\usepackage{html} 6 7% BEWARE: don't revert from graphicx for epsfig, because latex2html 8% doesn't handle epsfig commands !!! 9\usepackage{graphicx} 10 11\textwidth = 17cm 12\textheight = 24cm 13\topmargin =-1 cm 14\oddsidemargin = 0 cm 15 16\def\pwx{\texttt{pw.x}} 17\def\phx{\texttt{ph.x}} 18\def\configure{\texttt{configure}} 19\def\PWscf{\texttt{PWscf}} 20\def\PHonon{\texttt{PHonon}} 21\def\make{\texttt{make}} 22 23\begin{document} 24\author{} 25\date{} 26 27\def\qeImage{../../Doc/quantum_espresso} 28 29\title{ 30 \includegraphics[width=5cm]{\qeImage} \\ 31 % title 32 \Huge \PHonon\ User's Guide (v. \version) 33 \\ \Large (only partially updated) 34} 35 36\maketitle 37 38\tableofcontents 39 40\section{Introduction} 41 42This guide covers the usage of the \PHonon\ package, a 43part of the \qe\ distribution. 44Further documentation, beyond what is provided 45in this guide, can be found in the directory 46\texttt{PHonon/Doc/}, containing a copy of this guide. 47 48{\em Important notice: due to the lack of time and of manpower, this 49manual is only partially updated and may contain outdated information.} 50 51This guide assumes that you know the contents of 52the general User's Guide for \qe\ and of the User's 53Guide for \PWscf. It also assumes that you have 54already installed \qe\ (\PHonon\ is not a stand-alone 55package: it requires \PWscf\ to be compiled and used). 56If not, please locate the general User's Guide in directory 57\texttt{Doc/} two levels above the one containing this guide, 58and the User's Guide for \PWscf\ in \texttt{PW/Doc/}; 59or consult the web site:\\ 60\texttt{http://www.quantum-espresso.org}. 61It is also assumed that you know the physics behind \qe, 62the methods it implements, and in particular the physics 63and the methods of \PHonon. 64 65% People who want to modify or contribute to \PHonon\ should read 66% the Developer Manual: \texttt{Doc/developer\_man.pdf}. 67 68\PHonon\ has the following directory structure, 69contained in a subdirectory \texttt{PHonon/} 70of the main \qe\ tree: 71 72\begin{tabular}{ll} 73\texttt{Doc/} & : contains the user\_guide and input data description \\ 74\texttt{examples/} & : some running examples \\ 75\texttt{PH/} & : source files for phonon calculations 76 and analysis\\ 77\texttt{Gamma/} & : source files for Gamma-only phonon calculation\\ 78\end{tabular}\\ 79{\em Important Notice:} since v.5.4, many modules and routines that were 80common to all linear-response \qe\ codes are moved into the new 81\texttt{LR\_Modules} subdirectory of the main tree. Since v.6.0, the 82\texttt{D3} code for anharmonic force constant calculations has been 83superseded by the \texttt{D3Q} coda, available on 84\texttt{http://www.qe-forge.org/gf/project/d3q/}. 85 86The codes available in the \PHonon\ package can perform the following 87types of calculations: 88\begin{itemize} 89 \item phonon frequencies and eigenvectors at a generic wave vector, 90 using Density-Functional Perturbation Theory; 91 \item effective charges and dielectric tensors; 92 \item electron-phonon interaction coefficients for metals; 93 \item interatomic force constants in real space; 94 \item Infrared and Raman (nonresonant) cross section. 95\end{itemize} 96 97Phonons can be plotted using the \texttt{PlotPhon} package. 98Calculations of the vibrational free energy in the Quasi-Harmonic 99approximations can be performed using the \texttt{QHA} package. 100{\em Note:} since v.5.4, these two packages are separately distributed 101and no longer bundled with \PHonon. Their latest version can be 102found in the tarballs of \PHonon\ v.5.3. 103 104\section{People} 105The \PHonon\ package 106was originally developed by Stefano Baroni, Stefano 107de Gironcoli, Andrea Dal Corso (SISSA), Paolo Giannozzi (Univ. Udine), 108and many others. 109We quote in particular: 110\begin{itemize} 111 \item Michele Lazzeri (Univ.Paris VI) for the 2n+1 code and Raman 112 cross section calculation with 2nd-order response; 113 \item Andrea Dal Corso for the implementation of Ultrasoft, PAW, 114 noncolinear, spin-orbit extensions to \PHonon; 115 \item Mitsuaki Kawamura (U.Tokyo) for implementation of the optimized 116 tetrahedron method in phonon and electron-phonon calculations; 117 \item Thibault Sohier, Matteo Calandra, Francesco Mauri, for 118 phonons in two-dimensional systems; 119 \item Andrea Floris, Iurii Timrov, Burak Himmetoglu, Nicola Marzari, 120 Stefano de Gironcoli, and Matteo Cococcioni, for phonons using 121 Hubbard-corrected DFPT (DFPT+$U$). 122\end{itemize} 123 124The \texttt{PlotPhon} and \texttt{QHA} packages were contribute by the 125late Prof. Eyvaz Isaev. 126 127Other contributors include: Lorenzo Paulatto (Univ. Paris VI) for 128PAW, 2n+1 code; William Parker (Argonne) for phonon terms in dielectric 129tensor; Tobias Wassmann (Univ. Paris VI) for third-order derivatives of GGA 130potential. Weicheng Bao (Nanjing University) and Norge Cruz Hernandez 131(Universidad de Sevilla) helped debugging the phonon code. 132 133%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 134\input ../../Doc/quote.tex 135%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 136 137\section{Installation} 138 139\PHonon\ is a package tightly bound to \qe. 140For instruction on how to download and compile \qe, please refer 141to the general Users' Guide, available in file \texttt{Doc/user\_guide.pdf} 142under the main \qe\ directory, or in web site 143\texttt{http://www.quantum-espresso.org}. 144 145Once \qe\ is correctly configured, \PHonon\ can be automatically 146downloaded, unpacked and compiled by 147just typing \texttt{make ph}, from the main \qe\ directory. 148 149\subsection{Compilation} 150 151Typing \texttt{make ph} from the root \qe\ directory, or \texttt{make} 152from the \PHonon\ directory, produces the following codes: 153\begin{itemize} 154 \item \texttt{PH/ph.x}: Calculates phonon frequencies and displacement patterns, 155 dielectric tensors, effective charges (uses data produced by \pwx). 156 \item \texttt{PH/dynmat.x}: applies various kinds of Acoustic Sum Rule (ASR), 157 calculates LO-TO splitting at ${\bf q} = 0$ in insulators, IR and Raman 158 cross sections (if the coefficients have been properly calculated), 159 from the dynamical matrix produced by \phx 160 \item \texttt{PH/q2r.x}: calculates Interatomic Force Constants (IFC) in real space 161 from dynamical matrices produced by \phx\ on a regular {\bf q}-grid 162 \item \texttt{PH/matdyn.x}: produces phonon frequencies at a generic wave vector 163 using the IFC file calculated by \texttt{q2r.x}; may also calculate phonon DOS, 164 the electron-phonon coefficient $\lambda$, the function $\alpha^2F(\omega)$ 165\item \texttt{PH/lambda.x}: also calculates $\lambda$ and $\alpha^2F(\omega)$, 166 plus $T_c$ for superconductivity using the McMillan formula 167\item \texttt{PH/alpha2f.x}: also calculates $\lambda$ and $\alpha^2F(\omega)$. 168 It is used together with the optimized tetrahedron method 169 (\verb|occupations="tetrahedra_opt"| in \pwx) and shifted {\it q}-grid 170 (\verb|lshiftq=.true.| in \phx). 171\item \texttt{PH/fqha.x}: a simple code to calculate vibrational entropy with 172 the quasi-harmonic approximation 173\item \texttt{PH/dvscf\_q2r.x}: performs inverse Fourier transformation of phonon 174 potential from a regular {\bf q} grid to real space. 175\item \texttt{Gamma/phcg.x}: 176 a version of \phx\ that calculates phonons at ${\bf q} = 0$ using 177 conjugate-gradient minimization of the density functional expanded to 178 second-order. Only the $\Gamma$ (${\bf k} = 0$) point is used for 179 Brillouin zone integration. It is faster and takes less memory than 180 \phx, but does not support spin polarization, USPP and PAW. 181\end{itemize} 182Links to the main \qe\ \texttt{bin/} directory are automatically generated. 183 184\section{Using \PHonon} 185 186Phonon calculation is presently a two-step process. 187First, you have to find the ground-state atomic and electronic configuration; 188Second, you can calculate phonons using Density-Functional Perturbation Theory. 189Further processing to calculate Interatomic Force Constants, to add macroscopic 190electric field and impose Acoustic Sum Rules at ${\bf q}=0$ may be needed. 191In the following, we will indicate by ${\bf q}$ the phonon wavevectors, 192while ${\bf k}$ will indicate Bloch vectors used for summing over the 193Brillouin Zone. 194 195The main code \phx\ can be used whenever \PWscf\ can be used, 196with the exceptions of DFT+U, hybrid functionals, external electric fields, 197constraints on magnetization, nonperiodic boundary conditions. 198USPP and PAW are not implemented for higher-order response calculations. 199See the header of file \texttt{PHonon/PH/phonon.f90} for a complete and 200updated list of what \PHonon\ can and cannot do. 201 202Since version 4.0 it is possible to safely stop execution of \phx\ code using 203the same mechanism of the \pwx\ code, i.e. by creating a file 204\texttt{prefix.EXIT} in the working directory. Execution can be resumed by 205setting \texttt{recover=.true.} in the subsequent input data. 206Moreover the execution can be (cleanly) stopped after a given time is elapsed, 207using variable \texttt{max\_seconds}. See \texttt{example/Recover\_example/}. 208 209\subsection{Single-{\bf q} calculation} 210 211The phonon code \phx\ calculates normal modes at a given {\bf q}-vector, 212starting from data files produced by \pwx\ with a simple SCF calculation. 213NOTE: the alternative procedure in which a band-structure calculation 214with \texttt{calculation='phonon'} was performed as an intermediate step is no 215longer implemented since version 4.1. It is also no longer needed to 216specify \texttt{lnscf=.true.} for ${\bf q}\ne 0$. 217 218The output data files appear in the directory specified by the 219variable {\tt outdir}, with names specified by the variable 220{\tt prefix}. After the output file(s) has been produced (do not remove 221any of the files, unless you know which are used and which are not), 222you can run \phx. 223 224The first input line of \phx\ is a job identifier. At the second line the 225namelist {\tt \&INPUTPH} starts. The meaning of the variables in the namelist 226(most of them having a default value) is described in file 227\texttt{Doc/INPUT\_PH.*}. Variables \texttt{outdir} and \texttt{prefix} 228must be the same as in the input data of \pwx. Presently 229you can specify \texttt{amass(i)} (a real variable) the atomic mass 230of atomic type $i$ or you can use the default one deduced from the 231periodic table on the basis of the element name. If 232{\tt amass(i)} is not given as input of \phx, the one given as 233input in \pwx\ is used. When this is {\tt 0} the default one is used. 234 235After the namelist you must specify the {\bf q}-vector of the phonon mode, 236in Cartesian coordinates and in units of $2\pi/a$. 237 238Notice that the dynamical matrix calculated by \phx\ at ${\bf q}=0$ does not 239contain the non-analytic term occurring in polar materials, i.e. there is no 240LO-TO splitting in insulators. Moreover no Acoustic Sum Rule (ASR) is 241applied. In order to have the complete dynamical matrix at ${\bf q}=0$ 242including the non-analytic terms, you need to calculate effective charges 243by specifying option \texttt{epsil=.true.} to \phx. This is however not 244possible (because not physical!) for metals (i.e. any system subject to 245a broadening). 246 247At ${\bf q}=0$, use program \texttt{dynmat.x} to calculate the correct LO-TO 248splitting, IR cross sections, and to impose various forms of ASR. 249If \phx\ was instructed to calculate Raman coefficients, 250\texttt{dynmat.x} will also calculate Raman cross sections 251for a typical experimental setup. 252Input documentation in the header of \texttt{PHonon/PH/dynmat.f90}. 253 254See Example 01 for a simple phonon calculations in Si, Example 06 for 255fully-relativistic calculations (LDA) on Pt, Example 07 for 256fully-relativistic GGA calculations. 257 258\subsection{Calculation of interatomic force constants in real space} 259 260First, dynamical matrices are calculated and saved for a suitable uniform 261grid of {\bf q}-vectors (only those in the Irreducible Brillouin Zone of the 262crystal are needed). Although this can be done one {\bf q}-vector at the 263time, a 264simpler procedure is to specify variable \texttt{ldisp=.true.} and to set 265variables \texttt{nq1}, \texttt{nq2}, \texttt{nq3} to some suitable 266Monkhorst-Pack grid, that will be automatically generated, centered at 267${\bf q}=0$. 268 269Second, code \texttt{q2r.x} reads the dynamical matrices produced in the 270preceding step and Fourier-transform them, writing a file of Interatomic Force 271Constants in real space, up to a distance that depends on the size of the grid 272of {\bf q}-vectors. Input documentation in the header of \texttt{PHonon/PH/q2r.f90}. 273 274Program \texttt{matdyn.x} may be used to produce phonon modes and 275frequencies at any {\bf q} using the Interatomic Force Constants file as input. 276Input documentation in the header of \texttt{PHonon/PH/matdyn.f90}. 277 278See Example 02 for a complete calculation of phonon dispersions in AlAs. 279 280\subsection{Calculation of electron-phonon interaction coefficients} 281 282Since v.5.0, there are two ways of calculating electron-phonon 283coefficients, distinguished according to the value of variable 284\texttt{electron\_phonon}. The following holds for the case 285\texttt{electron\_phonon=} {\tt'interpolated'} (see also Example 03). 286 287The calculation of electron-phonon coefficients in metals is made difficult 288by the slow convergence of the sum at the Fermi energy. It is convenient to 289use a coarse {\bf k}-point grid to calculate phonons on a suitable 290wavevector grid; 291a dense {\bf k}-point grid to calculate the sum at the Fermi energy. 292The calculation 293proceeds in this way: 294\begin{enumerate} 295\item a scf calculation for the dense ${\bf k}$-point grid (or a scf calculation 296followed by a non-scf one on the dense ${\bf k}$-point grid); specify 297option \texttt{la2f=.true.} to \pwx\ in order to save a file with 298the eigenvalues on the dense {\bf k}-point grid. The latter MUST contain 299all ${\bf k}$ and ${\bf k}+{\bf q}$ grid points used in the subsequent 300electron-phonon 301calculation. All grids MUST be unshifted, i.e. include ${\bf k}=0$. 302\item a normal scf + phonon dispersion calculation on the coarse {\bf k}-point 303grid, specifying option \texttt{electron\_phonon='interpolated'}, and 304the file name where 305the self-consistent first-order variation of the potential is to be 306stored: variable \texttt{fildvscf}). 307The electron-phonon coefficients are calculated using several 308values of Gaussian broadening (see \texttt{PHonon/PH/elphon.f90}) 309because this quickly 310shows whether results are converged or not with respect to the 311{\bf k}-point grid and Gaussian broadening. 312\item Finally, you can use \texttt{matdyn.x} and \texttt{lambda.x} 313(input documentation in the header of \texttt{PHonon/PH/lambda.f90}) 314to get the $\alpha^2F(\omega)$ function, the electron-phonon coefficient 315$\lambda$, and an estimate of the critical temperature $T_c$. 316\end{enumerate} 317 318See the appendix for the relevant formulae. 319{\bf Important notice}: the $q\rightarrow 0$ limit of the contribution 320to the electron-phonon coefficient diverges for optical modes! please 321be very careful, consult the relevant literature. 322 323\subsection{DFPT with the tetrahedron method} 324 325In order to use the tetrahedron method for phonon calculations, 326you should run \pwx\ and \phx\ as follows: 327\begin{enumerate} 328 \item Run \pwx\ with \verb|occupation = "tetrahedra_opt"| and \verb|K_POINT automatic|. 329 \item Run \phx. 330\end{enumerate} 331 332There is an example in \verb|PHonon/example/tetra_example/|. 333 334\subsection{Calculation of electron-phonon interaction coefficients with the tetrahedron method} 335 336When you perform a calculation of electron-phonon interaction coefficients 337with the tetrahedron method, 338you have to use an offset $q$-point grid in order to avoid a singularity 339at $q=\Gamma$; you can perform this calculation as follows: 340 341\begin{enumerate} 342 \item Run \pwx\ with \verb|occupation = "tetrahedra_opt"| and \verb|K_POINT automatic|. 343 \item Run \phx\ with \verb|lshift_q = .true.| and \verb|electron_phonon = ""| (or unset it) 344 to generate the dynamical matrix and 345 the deformation potential (in \verb|_ph*/{prefix}_q*/|) of each $q$. 346 \item Run \phx\ with \verb|electron_phonon = "lambda_tetra"|. 347 You should use a denser $k$ grid by setting \verb|nk1|, \verb|nk2|, and \verb|nk3|. 348 Then \verb|lambda*.dat| are generated; they contain $\lambda_{q \nu}$. 349 \item Run \verb|alpha2f.x| with an input file as follows: 350\begin{verbatim} 351&INPUTPH 352! The same as that for the electron-phonon calculation with ph.x 353 : 354/ 355&INPUTA2F 356 nfreq = Number of frequency-points for a2F(omega), 357/ 358\end{verbatim} 359Then $\lambda$, and $\omega_{\ln}$ are computed and they are printed to the standard output. 360$\alpha^2F(\omega)$ and (partial) phonon-DOS are also computed; 361they are printed to a file \textit{prefix}\verb|.a2F.dat|. 362\end{enumerate} 363 364There is an example in \verb|PHonon/example/tetra_example/|. 365 366\subsection{Phonons for two-dimensional crystals} 367 368The extension of DFPT to two dimensional crystals, 369in particular gated two-dimensional heterostructure, 370s described in the following paper: 371 372T. Sohier, M. Calandra, and F. Mauri, Phys. Rev. B {\bf 96}, 075448 (2017), 373https://doi.org/10.1103/PhysRevB.96.075448 374 375See example \verb|PHonon/example/example17/|. 376 377\subsection{Phonons from DFPT+$U$} 378 379The extension of DFPT to inlcude Hubbard $U$ correction is described 380in the following papers: 381 382A.~Floris, S.~de~Gironcoli, E.~K.~U.~Gross, M.~Cococcioni, Phys. Rev. B {\bf 84}, 383161102(R), (2011); 384 385A.~Floris, I.~Timrov, B.~Himmetoglu, N.~Marzari, S.~de~Gironcoli, and M.~Cococcioni, 386Phys. Rev. B {\bf 101}, 064305 (2020). 387 388See example \verb|PHonon/example/example18/|. 389 390\subsection{Fourier interpolation of phonon potential} 391 392The potential perturbation caused by the displacement of a single atom is 393spatially localized. Hence, the phonon potential can be interpolated from 394$q$ points on a coarse grid to other $q$ points using Fourier interpolation. 395To use this functionality, first, one shoud run a \phx\ calculation for $q$ 396points on a regular coarse grid. 397Then, a \texttt{dvscf\_q2r.x} run performs an inverse Fourier 398transformation of the phonon potentials from a $q$ grid to a real-space 399supercell. Finally, by specifying \texttt{ldvscf\_interpolation=.true.} in 400\phx, the phonon potentials are Fourier transformed to given $q$ points. 401 402For insulators, the nonanalytic long-ranged dipole part of the potential needs 403to be subtracted and added before and after the interpolation, respectively. 404This treatment is activated by specifying \texttt{do\_long\_range=.true.} in the 405input files of \texttt{dvscf\_q2r.x} and \phx. 406 407Due to numerical inaccuracies, the calculated Born effective charges may not 408add up to zero, violating the charge neutrality condition. This error may lead 409to nonphysical polar divergence of the phonon potential for $q$ points close 410to $\Gamma$, even in IR-inactive materials. To avoid this problem, one can 411specify \texttt{do\_charge\_neutral=.true.} in the input files of 412\texttt{dvscf\_q2r.x} and \phx. Then, the phonon potentials and the Born 413effective charges are renormalized by enforcing the charge neutrality condition, 414following the scheme of S.~ Ponce et al, J. Chem. Phys. (2015). 415 416The Fourier interpolation of phonon potential is proposed and described in the 417following papers: 418 419A.~Eiguren and C.~Ambrosch-Draxl, Phys. Rev. B {\bf 78}, 045124 (2008); 420 421S.~ Ponce et al, J. Chem. Phys. {\bf 143}, 102813 (2015); 422 423X.~Gonze et al, Comput. Phys. Commun., {\bf 248}, 107042 (2020). 424 425\subsection{Calculation of phonon-renormalization of electron bands} 426The phonon-induced renormalization of electron bands can be computed using \PHonon. 427After SCF, PHONON, and NSCF calculations, one can run \phx\ with the 428\texttt{electron\_phonon=`ahc'} option, which generates binary files containing 429quantities required for the calculation of electron self-energy. 430Then, a \texttt{postahc.x} run reads these binary files and compute the 431phonon-induced electron self-energy at a given temperature. 432 433For more details, see the \verb|PHonon/Doc/dfpt_self_energy.pdf| file. 434 435Also, there is an example in \verb|PHonon/example/example19/|. 436 437Implementation of this functionality in \qe\ is described in the following 438paper: 439 440J.-M.~Lihm and C.-H.~Park, Phys. Rev. B {\bf 101}, 121102 (2020). 441 442 443\section{Parallelism} 444\label{Sec:para} 445 446We refer to the corresponding section of the \PWscf\ guide for 447an explanation of how parallelism works. 448 449\phx\ may take advantage of MPI parallelization on images, plane waves (PW) 450and on {\bf k}-points (``pools''). Currently all other MPI and explicit 451OpenMP parallelizations have very limited to nonexistent implementation. 452\texttt{phcg.x} implements only PW parallelization. 453All other codes may be launched in parallel, but will execute 454on a single processor. 455 456In ``image'' parallelization, processors can be divided into different 457``images", corresponding to one (or more than one) ``irrep'' or {\bf q} 458vectors. Images are loosely coupled: processors communicate 459between different images only once in a while, so image parallelization 460is suitable for cheap communication hardware (e.g. Gigabit Ethernet). 461Image parallelization is activated by specifying the option 462\texttt{-nimage N} to \phx. Inside an image, PW and {\bf k}-point 463parallelization can be performed: for instance, 464\begin{verbatim} 465 mpirun -np 64 ph.x -ni 8 -nk 2 ... 466\end{verbatim} 467will run $8$ images on $8$ processors each, subdivided into $2$ pools 468of $4$ processors for {\bf k}-point parallelization. In order 469to run the \phx\ code with these flags the \pwx\ run has to be run with: 470\begin{verbatim} 471 mpirun -np 8 pw.x -nk 2 ... 472\end{verbatim} 473without any {\tt -nimage} flag. 474After the phonon calculation with images the dynmical matrices of 475{\bf q}-vectors calculated in different images are not present in the 476working directory. To obtain them you need to run 477\phx\ again with: 478\begin{verbatim} 479 mpirun -np 8 ph.x -nk 2 ... 480\end{verbatim} 481and the {\tt recover=.true.} flag. This scheme is quite automatic and 482does not require any additional work by the user, but it wastes some 483CPU time because all images stops when the image that requires the 484largest amount of time finishes the calculation. Load balancing 485between images is still at 486an experimental stage. You can look into the routine {\tt image\_q\_irr} 487inside {\tt PHonon/PH/check\_initial\_status} to see the present 488algorithm for work distribution and modify it if you think that 489you can improve the load balancing. 490 491A different paradigm is the usage of the GRID concept, instead of MPI, 492to achieve parallelization over irreps and {\bf q} vectors. 493Complete phonon dispersion calculation can be quite long and 494expensive, but it can be split into a number of semi-independent 495calculations, using options \texttt{start\_q}, \texttt{last\_q}, 496\texttt{start\_irr}, \texttt{last\_irr}. An example on how to 497distribute the calculations and collect the results can be found 498in \texttt{examples/GRID\_example}. Reference:\\ 499{\it Calculation of Phonon Dispersions on the GRID using Quantum 500 ESPRESSO}, 501 R. di Meo, A. Dal Corso, P. Giannozzi, and S. Cozzini, in 502 {\it Chemistry and Material Science Applications on Grid Infrastructures}, 503 editors: S. Cozzini, A. Lagan\`a, ICTP Lecture Notes Series, 504 Vol. 24, pp.165-183 (2009). 505 506 507\section{Troubleshooting} 508 509\paragraph{ph.x stops with {\em error reading file}} 510The data file produced by \pwx\ is bad or incomplete or produced 511by an incompatible version of the code. 512 513\paragraph{ph.x mumbles something like {\em cannot recover} or {\em error 514 reading recover file}} 515You have a bad restart file from a preceding failed execution. 516Remove all files \texttt{recover*} in \texttt{outdir}. 517 518\paragraph{ph.x says {\em occupation numbers probably wrong} and 519 continues} You have a 520metallic or spin-polarized system but occupations are not set to 521\texttt{`smearing'}. 522 523\paragraph{ph.x does not yield acoustic modes with zero frequency at 524${\bf q}=0$} 525This may not be an error: the Acoustic Sum Rule (ASR) is never exactly 526verified, because the system is never exactly translationally 527invariant as it should be. The calculated frequency of the acoustic 528mode is typically less than 10 cm$^{-1}$, but in some cases it may be 529much higher, up to 100 cm$^{-1}$. The ultimate test is to diagonalize 530the dynamical matrix with program \texttt{dynmat.x}, imposing the ASR. If you 531obtain an acoustic mode with a much smaller $\omega$ (let us say 532$< 1 \mbox{cm}^{-1}$ ) 533with all other modes virtually unchanged, you can trust your results. 534 535``The problem is [...] in the fact that the XC 536energy is computed in real space on a discrete grid and hence the 537total energy is invariant (...) only for translation in the FFT 538grid. Increasing the charge density cutoff increases the grid density 539thus making the integral more exact thus reducing the problem, 540unfortunately rather slowly...This problem is usually more severe for 541GGA than with LDA because the GGA functionals have functional forms 542that vary more strongly with the position; particularly so for 543isolated molecules or system with significant portions of ``vacuum'' 544because in the exponential tail of the charge density a) the finite 545cutoff (hence there is an effect due to cutoff) induces oscillations 546in rho and b) the reduced gradient is diverging.''(info by Stefano de 547Gironcoli, June 2008) 548 549\paragraph{ph.x yields really lousy phonons, with bad or ``negative'' 550 frequencies or wrong symmetries or gross ASR violations} 551Possible reasons: 552\begin{itemize} 553\item if this happens only for acoustic modes at ${\bf q}=0$ that should 554 have $\omega=0$: Acoustic Sum Rule violation, see the item before 555 this one. 556\item wrong data file read. 557\item wrong atomic masses given in input will yield wrong frequencies 558 (but the content of file fildyn should be valid, since the force 559 constants, not the dynamical matrix, are written to file). 560\item convergence threshold for either SCF (\texttt{conv\_thr}) or phonon 561 calculation (\texttt{tr2\_ph}) too large: try to reduce them. 562\item maybe your system does have negative or strange phonon 563 frequencies, with the approximations you used. A negative frequency 564 signals a mechanical instability of the chosen structure. Check that 565 the structure is reasonable, and check the following parameters: 566\begin{itemize} 567\item The cutoff for wavefunctions, \texttt{ecutwfc} 568\item For USPP and PAW: the cutoff for the charge density, \texttt{ecutrho} 569\item The {\bf k}-point grid, especially for metallic systems. 570\end{itemize} 571\item For metallic systems: it has been observed that the convergence with 572 respect to the k-point grid and smearing is very slow in presence of 573 semicore states, and for phonon wave-vectors that are not commensurate i 574 with the k-point grid (that is, ${\bf q}\ne {\bf k}_i-{\bf k}_j$) 575\end{itemize} 576Note that ``negative'' frequencies are actually imaginary: the negative 577sign flags eigenvalues of the dynamical matrix for which $\omega^2 < 5780$. 579 580\paragraph{{\em Wrong degeneracy} error in star\_q} 581Verify the {\bf q}-vector for which you are calculating phonons. In order to 582check whether a symmetry operation belongs to the small group of ${\bf q}$, 583the code compares ${\bf q}$ and the rotated ${\bf q}$, with an acceptance tolerance of 584$10^{-5}$ (set in routine \texttt{PW/src/eqvect.f90}). You may run into trouble if 585your {\bf q}-vector differs from a high-symmetry point by an amount in that 586order of magnitude. 587 588\paragraph{Mysterious symmetry-related errors} 589Symmetry-related errors like {\em symmetry operation is non orthogonal}, 590or {\em Wrong representation}, or {\em Wrong degeneracy}, are almost 591invariably a consequence of atomic positions that are close to, 592but not sufficiently close to, symmetry positions. If such errors occur, 593set the Bravais lattice using the correct \texttt{ibrav} value (i.e. do 594not use \texttt{ibrav=0}), use Wyckoff positions if known. This must be 595done in the self-consistent calculation. 596 597\appendix 598\section{Appendix: Electron-phonon coefficients} 599 600\def\r{{\bf r}} 601\def\d{{\bf d}} 602\def\k{{\bf k}} 603\def\q{{\bf q}} 604\def\G{{\bf G}} 605\def\R{{\bf R}} 606 607\noindent The electron-phonon coefficients $g$ 608are defined as 609\begin{equation} 610g_{\q\nu}(\k,i,j) =\left({\hbar\over 2M\omega_{\q\nu}}\right)^{1/2} 611\langle\psi_{i,\k}| {dV_{SCF}\over d {\hat u}_{\q\nu} }\cdot 612 \hat \epsilon_{\q\nu}|\psi_{j,\k+\q}\rangle. 613\end{equation} 614The phonon linewidth $\gamma_{\q\nu}$ is defined by 615\begin{equation} 616\gamma_{\q\nu} = 2\pi\omega_{\q\nu} \sum_{ij} 617 \int {d^3k\over \Omega_{BZ}} |g_{\q\nu}(\k,i,j)|^2 618 \delta(e_{\q,i} - e_F) \delta(e_{\k+\q,j} - e_F), 619\end{equation} 620while the electron-phonon coupling constant $\lambda_{\q\nu}$ for 621mode $\nu$ at wavevector $\q$ is defined as 622\begin{equation} 623\lambda_{\q\nu} ={\gamma_{\q\nu} \over \pi\hbar N(e_F)\omega^2_{\q\nu}} 624\end{equation} 625where $N(e_F)$ is the DOS at the Fermi level. 626The spectral function is defined as 627\begin{equation} 628\alpha^2F(\omega) = {1\over 2\pi N(e_F)}\sum_{\q\nu} 629 \delta(\omega-\omega_{\q\nu}) 630 {\gamma_{\q\nu}\over\hbar\omega_{\q\nu}}. 631\end{equation} 632The electron-phonon mass enhancement parameter $\lambda$ 633can also be defined as the first reciprocal momentum of 634the spectral function: 635\begin{equation} 636\lambda = \sum_{\q\nu} \lambda_{\q\nu} = 6372 \int {\alpha^2F(\omega) \over \omega} d\omega. 638\end{equation} 639 640Note that a factor $M^{-1/2}$ is hidden in the definition of 641normal modes as used in the code. 642 643McMillan: 644\begin{equation} 645T_c = {\Theta_D \over 1.45} \mbox{exp} \left [ 646 {-1.04(1+\lambda)\over \lambda(1-0.62\mu^*)-\mu^*}\right ] 647\end{equation} 648or (better?) 649\begin{equation} 650T_c = {\omega_{log}\over 1.2} \mbox{exp} \left [ 651 {-1.04(1+\lambda)\over \lambda(1-0.62\mu^*)-\mu^*}\right ] 652\end{equation} 653where 654\begin{equation} 655\omega_{log} = \mbox{exp} \left [ {2\over\lambda} \int {d\omega\over\omega} 656 \alpha^2F(\omega) \mbox{log}\omega \right ] 657\end{equation} 658 659 660\end{document} 661