Electron Phonon Coupling

Here we show step-by-step how to use Quantum Espresso to calculate phonons and electron-phonon matrix-elements on a regular q-grid, with the final aim to allow Yambo to read these databases and calculate the temperature-dependent correction to the electronic states. This tutorial is quite complicated, take your time to follow all the steps

Electron-phonon matrix elements

In this first section we will explain how to generate electron-phonon matrix elements in Quantum-Espresso, and then import them in Yambo. Calculations will be divided in different folders:

• pseudo the pseudo potential folder
• scf for the self-consistent calculation
• nscf for the non-self-consistent calcaultion with a larger number of bands
• phonon for the phonons calculations
• dvscf for the calculation of electron-phonon matrix elements

In this tutorial we will show how to calculate electron-phonon induced corrections to the bands and optical properties of 2D hexagonal boron nitride. All input file are availabe in the following tgz file: hBN.epc.tgz

1. In scf we run a standard scf calculation choosing the a large k-grid in such a way to converge density. Do not forget to set force_symmorphic=.true., because not symmorphic symmetries are not supported yet in Yambo. Notice that because the present system is two-dimensional we added the flag assume_isolated="2D" in such a way correct phonons in 2D, remove this flag is you have a system with a different dimensionality (a bulk, a molecule etc...)

2. Go in the nscf folder, and then copy the ${PREFIX}.save folder from scf to nscf, in the present example just do cp -r ../scf/bn.save ./. In the nscf input you have to choose the number of k-points and bands you will use for the electron-phonon coupling and Yambo calculations, in our case we will a 9x9x1 grid and 8 bands.

.....   
 number of k points=    12
                      cart. coord. in units 2pi/alat
       k(    1) = (   0.0000000   0.0000000   0.0000000), wk =   0.0246914
       k(    2) = (   0.0000000   0.1283001   0.0000000), wk =   0.1481481
       k(    3) = (   0.0000000   0.2566001   0.0000000), wk =   0.1481481
       k(    4) = (   0.0000000   0.3849002   0.0000000), wk =   0.1481481
       k(    5) = (   0.0000000   0.5132002   0.0000000), wk =   0.1481481
       k(    6) = (   0.1111111   0.1924501   0.0000000), wk =   0.1481481
       k(    7) = (   0.1111111   0.3207501   0.0000000), wk =   0.2962963
       k(    8) = (   0.1111111   0.4490502   0.0000000), wk =   0.2962963
       k(    9) = (   0.1111111   0.5773503   0.0000000), wk =   0.1481481
       k(   10) = (   0.2222222   0.3849002   0.0000000), wk =   0.1481481
       k(   11) = (   0.2222222   0.5132002   0.0000000), wk =   0.2962963
       k(   12) = (   0.3333333   0.5773503   0.0000000), wk =   0.0493827
.....

3. Go in the phonon directory. You have to copy the k-points list in the 2pi/alat units from the previous nscf run and provide it as q-grid for the phonon calculations. Notice that due to an internal convection of Yambo the q-points have to be multiplied for -1 before add them to the phonons input. Then use the same k-point grid of the nscf calculation for the phonons, the final input will be:

&inputph
           verbosity = 'high'
              tr2_ph = 1e-12
              prefix = 'bn'
            fildvscf = 'bn-dvscf'
              fildyn = 'bn.dyn',
     electron_phonon = 'dvscf',
               epsil = .true.
               trans = .true.
               ldisp = .false.
               qplot = .true.
     nk1=9, nk2=9, nk3 = 1
/
12
   0.0000000    0.0000000   0.0000000  1
   0.0000000   -0.1283001   0.0000000  1
   0.0000000   -0.2566001   0.0000000  1
   0.0000000   -0.3849002   0.0000000  1
   0.0000000   -0.5132002   0.0000000  1
  -0.1111111   -0.1924501   0.0000000  1
  -0.1111111   -0.3207501   0.0000000  1
  -0.1111111   -0.4490502   0.0000000  1
  -0.1111111   -0.5773503   0.0000000  1
  -0.2222222   -0.3849002   0.0000000  1
  -0.2222222   -0.5132002   0.0000000  1
  -0.3333333   -0.5773503   0.0000000  1

4. Go in the dvscf folder. From the nscf folder copy density and wave-functions and from the phonon folder copy the dynamical matrices and the dvscf matrix elements: cp -r ../nscf/bn.save ./ and cp -r ../phonon/_ph0 ../phonon/bn*.dyn* ./. Then modify the phonon input to generate electron-phonon matrix elements for Yambo, by changing trans = .false. and electron_phonon = yambo. You do not need to specify the k-point grid because it is read from the nscf wave-functions. In this way we will not recalculate phonons, but only the electron-phonon matrix elements for the number of bands present in the nscf input file, in our case 8 bands:

&inputph
          verbosity = 'high'
             tr2_ph = 1e-12
             prefix = 'bn'
           fildvscf = 'bn-dvscf'
             fildyn = 'bn.dyn'
    electron_phonon = 'yambo',
              epsil = .true.
              trans = .false.
              ldisp = .false.
              qplot = .true.
/
12
   0.0000000    0.0000000   0.0000000  1
   0.0000000   -0.1283001   0.0000000  1
   0.0000000   -0.2566001   0.0000000  1
   0.0000000   -0.3849002   0.0000000  1
   0.0000000   -0.5132002   0.0000000  1
  -0.1111111   -0.1924501   0.0000000  1
  -0.1111111   -0.3207501   0.0000000  1
  -0.1111111   -0.4490502   0.0000000  1
  -0.1111111   -0.5773503   0.0000000  1
  -0.2222222   -0.3849002   0.0000000  1
  -0.2222222   -0.5132002   0.0000000  1
  -0.3333333   -0.5773503   0.0000000  1

this run will generate a new folder elph_dir with all electron-phonon matrix elements in a format compatible with Yambo.

5. Now you have to read the wave-functions. Go in the dvscf/bn.save folder, do p2y, then yambo_ph -i -V kpt and turn on the flag BSEscatt, then run the setup. Now you use ypp_ph to import electron-phonon coupling, by doing ypp_ph -g' and the the PW el-ph folder:

gkkp                             # [R] gkkp databases
DBsPATH= "../elph_dir/"                     # Path to the PW el-ph databases
PHfreqF= "none"                  # PWscf format file containing the phonon frequencies
PHmodeF= "none"                  # PWscf format file containing the phonon modes
#GkkpExpand                    # Expand the gkkp in the whole BZ
#GkkpConvert                   # Convert the gkkp to new I/O format

run ypp_ph, and if everything went well you will get in output

.....
-> [06] == Electron-Phonon Interface: PW->Yambo Databases ==
<---> PW(ELPH) databases ...[PHONON] ...found 12 Q-grid compatible
<---> ELPH databases (WRITE) |########################################| [100%] --(E) --(X)
<--->  Modes           :   6
<--->  Bands range     :   8
.....

Quasi-particle band structure

The first quantity we will calculate is the correction to the band structure induced by the electron-phonon coupling. In order to generate the corresponding input do yambo_ph -g n -p fan -c p -V gen (in the newer versions of Yambo this line will be changed in yambo_ph -g n -p ep -c p -V gen):

dyson                            # [R] Dyson Equation solver
gw0                              # [R] GW approximation
el_ph_corr                       # [R] Electron-Phonon Correlation
Nelectro= 8.000000               # Electrons number
ElecTemp= 0.000000         eV    # Electronic Temperature
BoseTemp= 0.000000         eV    # Bosonic Temperature
OccTresh= 0.100000E-4            # Occupation treshold (metallic bands)
SE_Threads=0                     # [OPENMP/GW] Number of threads for self-energy
DysSolver= "n"                   # [GW] Dyson Equation solver ("n","s","g")
GphBRnge=  8                     # [ELPH] G[W] bands range
% ElPhModes
  1 |  6 |                           # [ELPH] Phonon modes included
%
RandQpts=0                       # [RIM] Number of random q-points in the BZ
#WRgFsq                        # [ELPH] Dump on file gFsq coefficients
%QPkrange                        # [GW] QP generalized Kpoint/Band indices
1|12|1|8|
%

If you run yambo_ph -J T0 you will get correction at zero kelvin to the band structure. You can change the temperature to 300 K by doing:

dyson                            # [R] Dyson Equation solver
gw0                              # [R] GW approximation
el_ph_corr                       # [R] Electron-Phonon Correlation
Nelectro= 8.000000               # Electrons number
ElecTemp= 0.000000         eV    # Electronic Temperature
BoseTemp= 300.0            K    # Bosonic Temperature
OccTresh= 0.100000E-4            # Occupation treshold (metallic bands)
SE_Threads=0                     # [OPENMP/GW] Number of threads for self-energy
DysSolver= "n"                   # [GW] Dyson Equation solver ("n","s","g")
GphBRnge=  8                     # [ELPH] G[W] bands range
% ElPhModes
  1 |  6 |                           # [ELPH] Phonon modes included
%
RandQpts=0                       # [RIM] Number of random q-points in the BZ
#WRgFsq                        # [ELPH] Dump on file gFsq coefficients
%QPkrange                        # [GW] QP generalized Kpoint/Band indices
1|12|1|8|
%

run again yambo_ph -J T300 and compare the two output file o-T0.qp and o-T300.qp, to find how much the gap change with the temperature.

Optical properties

You can use the previous calculate quasi-particle band structure to get optical properties, by doing yambo_ph -o b -k sex -y d -V qp -X s and then reading the corresponding ndb.QP database KfnQPdb= "E < SAVE/ndb.QP"

Nota bene: Calculation of exciton width is not completely correct. If fact electron-phonon matrix elements should be rotate in the excitonic space, as it was done in Ref. Phys. Rev. Lett. 125 107401 (2020), while at present Yambo calculate the exciton width from the correction to the single particle bands, see ref. Phys. Rev. Lett. 101 106405 (2008). This approximation is correct in case of not-strong bounded excitons, and unfortunately this is not the case.

Convergence

The results of this tutorial are not converged, and if you check output there are also some phonon mode with negative energies. This is due to the poor parameters used in this tutorial. In order to have converged results, first of all be sure to have converged phonons, increase the plane-wave cutoff, the number of k-points and if necessary reduce tr2_ph. Then change the paramters for the Yambo calculations, increasing the number of k and q points and the number of bands. If you want to increase only the number of bands, just repeat the nscf and dvscf calculations without recalculate the phonons.