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[1] 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 calculation 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 for bulk silicon. All input file are availabe in the following tgz file: si.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 if you are working with 2D systems, not our case, we have to add the flag assume_isolated="2D" in such a way correct 2D phonons.
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/si.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 4x4x4 grid and 12 bands.
In the nscf output you will find the list of k-points corresponding to the 4x4x2 grid:
.....
number of k points= 8
cart. coord. in units 2pi/alat
k( 1) = ( 0.0000000 0.0000000 0.0000000), wk = 0.0312500
k( 2) = ( 0.1250000 0.1250000 -0.1250000), wk = 0.2500000
k( 3) = ( -0.2500000 -0.2500000 0.2500000), wk = 0.1250000
k( 4) = ( 0.2500000 0.0000000 0.0000000), wk = 0.1875000
k( 5) = ( -0.1250000 -0.3750000 0.3750000), wk = 0.7500000
k( 6) = ( 0.0000000 -0.2500000 0.2500000), wk = 0.3750000
k( 7) = ( -0.5000000 0.0000000 0.0000000), wk = 0.0937500
k( 8) = ( -0.5000000 0.2500000 0.0000000), wk = 0.1875000
.....
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 use 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. You can use this simple python script minuq_kpoints.py to read k-points from the QE output and multiply them by minus one. The final input will be:
&inputph
verbosity = 'high'
tr2_ph = 1e-12h
prefix = 'si'
fildvscf = 'si-dvscf'
fildyn = 'si.dyn',
electron_phonon = 'dvscf',
epsil = .true.
trans = .true.
ldisp = .false.
qplot = .true.
/
8
0.000000000 0.000000000 0.000000000 1
-0.125000000 -0.125000000 0.125000000 1
0.250000000 0.250000000 -0.250000000 1
-0.250000000 0.000000000 0.000000000 1
-0.375000000 -0.125000000 0.125000000 1
0.000000000 0.250000000 -0.250000000 1
0.500000000 0.000000000 0.000000000 1
0.500000000 -0.250000000 0.000000000 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/si.save ./ and cp -r ../phonon/_ph0 ../phonon/*.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 = 'si'
fildvscf = 'si-dvscf'
fildyn = 'si.dyn'
electron_phonon = 'yambo',
epsil = .true.
trans = .false.
ldisp = .false.
qplot = .true.
/
8
0.000000000 0.000000000 0.000000000 1
-0.125000000 -0.125000000 0.125000000 1
0.250000000 0.250000000 -0.250000000 1
-0.250000000 0.000000000 0.000000000 1
-0.375000000 -0.125000000 0.125000000 1
0.000000000 0.250000000 -0.250000000 1
0.500000000 0.000000000 0.000000000 1
0.500000000 -0.250000000 0.000000000 1
this run will generate a new folder elph_dir with all electron-phonon matrix elements in a format compatible with Yambo.
Nota bene: use the same parallelization in the phonon and dvscf calculation to avoid problem reading wave-functions.
5. Now you have to read the wave-functions. Go in the dvscf/si.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 8 Q-grid compatible <---> ELPH databases (WRITE) |########################################| [100%] --(E) --(X) <---> Modes : 6 <---> Bands range : 12 .....
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). In the input we require the correction only at gamma point:
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= 12 # [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|1|4|5|
%
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= 12 # [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|1|4|5| %
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. If you repeat the calculation for different temperature you can obtain the trend of the gap vs temperature, see figure below:
If you uncomment the flag WRgFsq the code saves information about the Eliashberg functions that can be plotted using the ypp_ph, see below.
Finally if you add -V qp in the input generation a new flag appears OnMassShell, if you un-comment this flag calculation will be performed in the "on mass shell" approximation, namely the static limit the Quasi-Particle approximation, for a discussion see reference [2]
Optical properties
Now you repeat the previous calculation but including all k-points, the last 3 valence and the first 3 conduction bands:
..... %QPkrange # [GW] QP generalized Kpoint/Band indices 1|8|2|7| % ....
and save the results of the 0K and 300K temperature in two separate folder with the -J option. Now you can use the correction to the energy levels and the induced width to calculate the optical absorption at finite temperature. Generate the input with the command yambo_ph -o c -V qp
optics # [R] Linear Response optical properties
chi # [R][CHI] Dyson equation for Chi.
dipoles # [R] Oscillator strenghts (or dipoles)
Chimod= "IP" # [X] IP/Hartree/ALDA/LRC/PF/BSfxc
% QpntsRXd
1 | 1 | # [Xd] Transferred momenta
%
% BndsRnXd
2 | 7 | # [Xd] Polarization function bands
%
% EnRngeXd
0.000000 | 5.000000 | eV # [Xd] Energy range
%
% DmRngeXd
0.0100000 | 0.0100000 | eV # [Xd] Damping range
%
ETStpsXd= 500 # [Xd] Total Energy steps
% LongDrXd
1.000000 | 0.000000 | 0.000000 | # [Xd] [cc] Electric Field
%
XfnQPdb= "E W < T300/ndb.QP" # [EXTQP Xd] Database action
set the path of the ndb.QP you want to read and perform the calculations. Notice that from the QP database we read two quantities the correction to the energy levels E and the width W. In this calculation we also included a small smearing 0.01eV to mimic the electronic smearing. Hereafter the result without and with electron-phonon coupling at two different temperatures:
The same procedure can be applied to the Bethe-Salpeter calculations. However in this case results can be slightly wrong. If fact electron-phonon matrix elements should be rotate in the excitonic space, as it was done in ref. [3], while at present Yambo calculate the exciton width from the correction to the single particle bands, see ref. [4] . 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. This is due to the poor parameters used in this tutorial to make the calculations fast. In order to have converged results, first of all you have to be sure to have converged phonons, in order to do that increase plane-wave cutoff, number of k-points and if necessary reduce tr2_ph.
Then change the parameters 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 phonons. In the following section we will describe a smart way to accelerate convergence.
Nota bene: At present Yambo neither implements the Fröhlich at q=0[5] nor the quadupolar correction,[6] therefore convergence in polar material can be very slow with the number of q-points.
Double-grid method for the electron-phonon coupling
To be done
Miscellaneous and post-processing
There are a list of utilities to analyze electron-phonon coupling results:
Phonon density of states
You can plot phonon density of states with the command: ypp_ph -p d. In order to get a nice plot set in the input
phonons # [R] Phononic properties dos # [R] DOS PhBroad= 0.0005000 eV # Phonon broadening (Eliashberg & DOS) PhStps=1000 # Energy steps
notice that this is an easy quantity to check for the convergence in q-points.
Eliashberg Functions
You can plot Eliashberg functions[7] for both electrons and excitons. In order to plot Eliashberg functions you must have calculated Quasi-Particle correction with the flag WRgFsq, see above.
The command ypp_ph -s e generate the input for the electronic Eliashberg functions:
electrons # [R] Electronic properties eliashberg # [R] Eliashberg PhBroad= 0.0010000 eV # Phonon broadening (Eliashberg & DOS) PhStps= 200 # Energy steps %QPkrange # generalized Kpoint/Band indices 1|1|4|5| %
in this example we plot Eliashberg functions for the top valence and bottom conduction band at Gamma point:
For the excitonic Eliashberg functions use the command ypp_ph -e e, an example of these functions can be found in ref.[4].
Atomic displacement amplitudes
Running ypp_ph -p a will plot the atomic displacement for each atom in the cell each direction.
Other variables in the input files
In the input files of the present tutorial there are other variable not used in this tutorial.
In particular GkkpExpand is used for other calculations not present in the GPL, while GkkpConvert is used to checking purpose the IO for different versions
of the databases.
References
- ↑ This tutorial is based on Elena Cannuccia blog
- ↑ H. Kawai et al. Phys. Rev. B 89, 085202 (2014)
- ↑ H. Chen, D. Sangalli, and M. Bernardi, Phys. Rev. Lett. 125 107401 (2020)
- ↑ 4.0 4.1 A. Marini, Phys. Rev. Lett. 101 106405 (2008)
- ↑ Carla Verdi and Feliciano Giustino Phys. Rev. Lett. 115, 176401
- ↑ G. Brunin et al. Phys. Rev. Lett. 125, 136601
- ↑ F. Marsiglio, J.P. Carbotte Electron - Phonon Superconductivity