The coupling between the electronic and the atomic degrees of freedom plays a key role in several physical phenomena. For example it affects the temperature dependence of carriers mobility in organic devices or the position and intensity of Raman peaks. The electron-phonon coupling is also the driving force that causes excitons dissociation at the donor/acceptor interface in organic photovoltaic and the transition to a superconducting phase in molecular solids. From the theoretical point of view the role of the atomic vibrations has often been treated in a semi-empirical manner. Such approaches, based on model hamiltonians, rely on parameters that are difficult to extract from experiments. In contrast the ab-initio methods describe, and in some cases predict, in a quantitative manner the optical and electronic properties of many different systems, without resorting to any external parameter.
In this tutorial we will describe how to calculate the electron-phonon induced renormalization of the band-gap of bulk silicon and diamond. We will
also investigate the electronic line-widths induced by the electron-phonon coupling and how these can be used to interpret the finite
temperature absorption spectrum of bulk Silicon.
>yambo_ph -J 02_OMS
>yambo_ph -J 02_OMS -F INPUTS/02_OMS
Silicon is an indirect gap seminconductor. If you want to calculate the correction to the indirect gap you have to locate the Conduction Band Minimum (CBM) and Valence Band Maximum (VBM).
In the k-grid provided with this tutorial the Γ point is point #1, while the X point (the nearest to the CBM) is at point #7.
We start by entering in the Electron-Phonon/Section_A folder relative to the folder YAMBO_TUTORIALS. The list of folders gives
Section_A>ls
DFT/ DFT_commons/ DFT_dVscf/ DFT_el-ph/ INPUTS/ REFERENCE/ SAVE/
The folders starting with DFT contains relative scripts to run the Quantum-Espresso code in order to calculate the electron-phonon matrix elements and phonon frequencies. We will review shortly this procedure in section E.
We start, therefore, with the usual SAVE already containing the electron-phonon databases imported from QE:
Section_A>ls SAVE/
ndb.elph_gkkp ndb.elph_gkkp_fragment_2 ndb.elph_gkkp_fragment_4
ndb.elph_gkkp_fragment_6 ndb.elph_gkkp_fragment_8 ns.wf
ndb.elph_gkkp_fragment_1 ndb.elph_gkkp_fragment_3
ndb.elph_gkkp_fragment_5 ndb.elph_gkkp_fragment_7 ns.db1
The ndb.elph_gkkp* files contain the and matrix elements and the phonon frequencies needed to calculate the Fan and Debye Waller self-energies.
In the INPUTS folder all input files are provided. However we prefer to ask you to create yourself the input files. So let's start with the initialization that we do by typing
>yambo_ph -i -V kpt
then, in the input file, uncomment the flag
#MinusQ # [KPT] Use -{q} grid
in order to force yambo to append a (-1) to the automatic provided q-grid. This is because of the different convention between Yambo and QE in defining the matrix elements. More precisely yambo adopts the
|k> -> |k-q> (Yambo)
|k> -> |k+q> (QE)
By using the MinusQ flag Yambo and QE conventions will be the same. You will better understand this difference in the next section. Now we can run Yambo being careful in using a proper JOB naming convention to spot more easily the reference files. So let's use
>yambo_ph -J 01_init
To calculate the QP energies in the pure HAC we must use the combination of command options
>yambo_ph -g n -s p -V qp
The -V qp option is needed to activate some additional flags/variables that are switched off by default. The one we need now is
#OnMassShell # [F GW] On mass shell approximation
This flag is needed to perform a On-The-Mass-Shell calculation, that is a static limit of the dynamical theory. So, now, please edit your input file
(that, if you did not provide any -I
>yambo_ph -J 02_OMS
In this run Yambo will evaluate the Fan and DW self-energies at zero temperature by imposing the On-The-Mass-Shell approximation in order to get the QP energies defined as
The result of the calculation is reported in the report file and in the output file o-02_OMS.qp. The gap in Silicon is indirect and, in the k-points provided with this database, is given by transitions between the 1^{st} and the 7^{th} point:
#
# K-point Band Eo Eqp E-Eo Z Width[meV]
#
1.00000 4.00000 0.00000 0.03942 0.03942 1.00000 12.16324
7.00000 5.00000 0.60780 0.57831 -0.02949 1.00000 -9.01723
We notice that: (i) the gap correction is negative and of about 50 meV, (ii) the lifetimes sign changes, (iii) the renormalization factor Z is equal to one.
EXERCISE: Why is the sign of the lifetime different. Will it be the same at finite temperature ? Why is Z=1?
The electron-phoon matrix elements diverge as q → 0. This divergenve is compensated by the volume factor coming from the d^{3}q integration prefactor.
We can now create anothet input file or edit the previous one and enter a non zero value for the input file variable
RandQpts=100000 # [RIM] Number of random q-points in the BZ
and run
>yambo_ph -J 03_OMS_RIM
This variable controls the number of random q-points to be used to integrate via a Montecarlo technique the integral of the divergence.
#
# K-point Band Eo Eqp E-Eo Z Width[meV]
#
1.00000 4.00000 0.00000 0.03260 0.03260 1.00000 -0.05072
7.00000 5.00000 0.60780 0.57854 -0.02926 1.00000 -7.74586
In this case the effect of the Random Integration is small. This scenario changes drammatically in system with a compressed Brillouin zone, like nano--structures or surfaces.
The generalized Eliashberg functions, defined as,
with the residuals defined as
The generalized Eliashberg functions allows to visualize which phonon modes contribute to the electronic energy renormalzation. Indeed it is easy to realize that
In order to calculate the generalized Eliashberg functions we need to use the ypp_ph post-processor by using the command line
>ypp_ph -s e
then edit the ypp.in input file by changing the variable values in order to be the same of the example contained in INPUTS/04_eliashberg and run
>ln -s 02_OMS 04_eliashberg
>ypp_ph -J 04_eliashberg
The first command is to give to ypp_ph the input databases where to read from the eliashberg functions. At this point you can plot the two eliashberg functions corresponding to the conduction band minimum (CBM) amd valence band maximum (VBM). Use the columns 1:3 (Fan), 1:4 (DW) and 1:5 (SUM).
From these plots it is evident the cancellation between the DW and Fan contribution in the low-energy part of the spectrum. This is linked to the role played by the two contributions in order to satisfy the translational invariance of the theory.
Altough this tutorial is mainly based on the Heine-Allen-Cardona approach (HAC) it is istructive to see whether a dynamical approach agrees or not with the HAC approach. Therefore in this part of the tutorial we go beyond the HAC approach by simply redo the same procedure as before commenting the OnMassShell. The reference input file to use/create is INPUTS/05_QP
#OnMassShell # [F GW] On mass shell approximation
We, then, run
>yambo_ph -J 05_QP
The ouput is, now, more structured. Indeed there are several columns
#
# K-point Band Eo Eqp E-Eo Sc(Eo) Z Width[meV] Width[fs]
#
1.000 1.000 -12.08 -12.11 -.2396E-1 .2429E-1 0.9861 0.3104E-1 0.2121E+5
1.000 2.000 0.000 0.3136E-1 0.3136E-1 .3268E-1 0.9596 -.1237 -5321.
1.000 3.000 0.000 0.3136E-1 0.3136E-1 .3268E-1 0.9596 -.1237 -5321.
1.000 4.000 0.000 0.3136E-1 0.3136E-1 .3268E-1 0.9596 -.1237 -5321.
1.000 5.000 2.557 2.568 0.1124E-1 .9461E-2 0.9522 -23.69 -27.78
1.000 6.000 2.557 2.568 0.1124E-1 .9461E-2 0.9522 -23.69 -27.78
1.000 7.000 2.557 2.568 0.1124E-1 .9461E-2 0.9522 -23.69 -27.78
1.00000 8.00000 3.53155 3.46545 -0.06610 -0.06234 1.06837 -17.75930 -37.06294
1.000 9.000 7.759 7.742 -.1691E-1 .1773E-1 0.9556 -2.574 -255.7
1.000 10.00 7.896 7.906 0.1003E-1 .1060E-1 0.9511 -3.135 -210.0
7.000 1.000 -7.855 -7.879 -.2402E-1 .1887E-1 1.134 41.78 15.75
7.000 2.000 -7.855 -7.879 -.2402E-1 .1887E-1 1.134 41.78 15.75
7.000 3.000 -2.944 -2.938 0.6718E-2 .6969E-2 0.9631 2.213 297.4
7.000 4.000 -2.944 -2.938 0.6718E-2 .6969E-2 0.9631 2.213 297.4
7.000 5.000 0.6078 0.5770 -.3079E-1 .2933E-1 1.047 -7.859 -83.75
7.000 6.000 0.6078 0.5770 -.3079E-1 .2933E-1 1.047 -7.859 -83.75
7.000 7.000 10.22 10.13 -.8477E-1 .5134E-1 1.480 -57.70 -11.41
7.000 8.000 10.22 10.13 -.8477E-1 .5134E-1 1.480 -57.70 -11.41
7.000 9.000 11.22 11.22 -.7947E-4 .6968E-2 1.120 -54.64 -12.05
7.000 10.00 11.22 11.22 -.7947E-4 .6968E-2 1.120 -54.64 -12.05
We notice immediately the some renormalization factors are even unphysical.
EXERCISE: Why ? How can you explain using simple arguments the appearance of these anomalous Z factors ?EXERCISE: Repeat the calculation of the point A.1: at room temperature. In practice this means to add the field
BoseTemp= 300 K # Bosonic Temperature
to the INPUTS/02_OMS input file.
EXERCISE: Does the gap correction increase or decrease at room temperature ? And why ? [Tip. Plot the eliashberg functions to get a deeper understanding]
This section of the tutorial contains a crucial difference with the previous one. The transferred momenta are not relative to a uniform grid but are, instead, generated randomly in the Brillouin zone. This procedure allows a better convergence with the number of points used. From the numerical point of view the difference is that in the definition of the the |k+q> is not a point of the k grid.
This procudure is particularly important for the calculation of the el-ph induced lifetimes, that is the subject of the present part of the tutorial.
We start by entering in the Electron-Phonon/Section_B folder relative to the folder YAMBO_TUTORIALS. The list of
folders gives
Section_B>ls
DFT/ DFT_commons/ DFT_dVscf/ DFT_el-ph/ INPUTS/ REFERENCE/ SAVE/
In this case the SAVE folder contains 50 ndb.elph_gkkp_fragment databases corresponding to 50 random q-points. This is different from the previous case as the given grid of k-points produces a total of 8 q-points. And this is indeed the number of ndb.elph_gkkp_fragment databases provided in the section A tutorial.
Please proceed in the same way of section A.1
Please proceed in the same way of section A.2 and A.3 by also testing the effect of the Random Integration around the gamma point.
Also in this case the procedure is the same of section A.4 but the use of random points makes the results much smoother and the cancellation of the Fan and DW terms more evident. We plot, then, again the two eliashberg functions corresponding to the conduction band minimum (CBM) amd valence band maximum (VBM). Use the columns 1:3 (Fan), 1:4 (DW) and 1:5 (SUM).
Here we can use the HAC command line already used in Section A:
>yambo_ph -g n -s p -V qp
and uncomment the #OnMassShell option. We remind again that the input files correspinding to this section are in INPUTS named 05_OMS_RIM_100K, 07_OMS_RIM_297K, 08_OMS_RIM_400K, 09_OMS_RIM_465K, 10_OMS_RIM_600K. To reproduce these inputs by hand remember to add the temperature option
BoseTemp= 300 K # Bosonic Temperature
corresponding to the temperature you are interested in. Note that the BoseTemp appears by using the -V gen verbosity option. Now take the values of the linewidth at one band/k-point from the o-02....qp files and create a table using as a first column the temperature in Kelvin. Use this table to plot the linewidths as a function of the temperature. You should obtain a plot like the following where we used the band 5^{th} at the gamma point.
The aim of this part of the tutorial is to reproduce the finite temperature absorption spectra of bulk Silicon previously published by A. Marini
To this end we enter the Electron-Phonon/Section_C folder after having performed the Section_B part of this tutorial as the calculated finite temperature linewidths will be used as input.
So let's start with the initialization that we do by typing
>yambo -i
>yambo
We first rename the precompiled screened interaction database to speed up the calculation
>cd SAVE
>mv EM1S ndb.em1s
Now we calculate the absorption spectrum in the standard approach by using
>yambo -o b -k sex -y i
>yambo -J 02_bse
The next step is to use a randomly generated list of k-points (produced using PWscf and the input files contained in the DFT_E_RIM/ folder) to make the energy dependence of the Green's functions smoother. The idea is the following. The BSE is a linear equation of the electron-hole propagator written in the electron-hole basis
where
When the RIM is used the definition of the bare propagator is
The random k-points enter in the sum by turning the spectrum in a much smoother function of the frequency. To run the BSE by using the RIM we can use the kernel prevously calculated and we have to give as input the ndb.E_RIM database:
>cp 02_bse/ndb.BS_Q1 SAVE
>cp RIM/ndb.E_RIM SAVE
>yambo -J 03_bse_rim
The comparison of the output of the two calculations yields:
And now the last step. We want to make the absorption spectrum temperature dependent. The most important contribution comes from the poles of the non-interacting electron-hole Green's functtion. We make Yambo, now, to read the polaronic energies and widths calculated in Section B. To this end we perform a standard BSE calculation, like in the C.3 case. But with to important modifications. We add, indeed, to the input file the lines (changing the damping section)
KfnQPdb = "E W < ../Section_B/02_OMS/ndb.QP" # [EXTQP BSK BSS] Database"
% BDmRange
0.01000 | 0.01000 | eV # [BSS] Damping range
%
The first line makes yambo to read the QP energies and linewidths and the second puts (almost) to zero the internal damping so that the one observed is the real electron-phonon induced one. Now we can run the code
>yambo -J 04_bse_rim_0K
We can now perform the same steps for the other temperatures used in the Section B. Each time it is enough to include in the input file the path to the correspoding QP database.
Keeping in mind that this tutorial represents an highly uncoverged calculation and ALL parameters should be carefully checked and converged the agreement is really good!
In this section we calculate the QP energies within the HAC approximation and beyond (including dynamical effects) in the case of Bulk Diamond. We will consider, in the following, the correction to the optical gap of Diamond that is located at the Γ point, that is band #4-#5 at k-point #1. We start by entering in the Electron-Phonon/Section_D folder relative to the folder YAMBO_TUTORIALS. At this point we ask you to go trough the same steps of Section A by repeating EXACTLY the same calculations.
The result you should find by using/creating the INPUTS/03_OMS_RIM for your el-ph induced correction is
#
# K-point Band Eo Eqp E-Eo Z Width[meV]
#
1.000 4.000 0.000 0.1449 0.1449 1.000 0.1273E-4
1.000 5.000 5.653 5.28799 -0.36467 1.00000 -64.61224
This corresponds to a huge gap correction of 600 meV (!!!). This is much larger than the bulk silicon case and in reasonable agreement with the number publised by E. Cannuccia and F. Giustino results.
A deeper insight in this huge gap correction can be obtained by using the genaralized Eliashberg functions that in the case of Diamond look like
We notice immediately two maior points: (i) the Eliashberg functions extend on a wider energy range thanks to the larger Debye Energy of Diamond, (ii) the intensity of the functions is much larger than in the case of Silicon.
First of all download the Quantum Espresso 4.0.5 source and the files that have been changed/added in order to calculate the electron-phonon matrix elements. Please note that in the near future these changes will be embodied in the standard quantum-espresso distribution.
After having unpacked the source
>ls
espresso-4.0.5/
enter the espresso-4.0.5/ folder and unzip the patched files
>unzip NEWFILES_ELPH_4.0.5.zip
Archive: NEWFILES_ELPH_4.0.5.zip
replace Modules/input_parameters.f90? [y]es, [n]o, [A]ll, [N]one, [r]ename: A
inflating: Modules/input_parameters.f90
inflating: Modules/read_namelists.f90
inflating: PW/sgam_at.f90
inflating: Modules/Makefile
inflating: PH/allocate_phq.f90
inflating: PH/bcast_ph_input.f90
inflating: PH/deallocate_phq.f90
inflating: PH/elphon.f90
inflating: PH/Makefile
inflating: PH/phonon.f90
inflating: PH/phq_init.f90
inflating: PH/phq_readin.f90
inflating: Modules/yambo.f90
inflating: PH/debye_waller.f90
then follow the standard compilation procedure.
The procedure to get the el-ph databases ndb.elph_gkkp_fragment will be described in this section shortly