The Yambo Project - User contributions [en]

Phonon-assisted luminescence by finite atomic displacements

2022-04-26T14:06:45Z

Pierre: /* Calculate phonons at momentum q of the indirect transition */

[[File:Representative-instantaneous-atomic-configurations-of-H2O-top-and-H3O-bottom-Also.gif|thumb| 200px|Finite atomic displacements]]

In recent years, several articles have been published with Yambo, who calculated the phonon-assisted luminescence by means of finite atomic displacements<ref name='cann'>[https://arxiv.org/abs/1807.11797 Theory of phonon-assisted luminescence in solids: Application to hexagonal boron nitride], E. Cannuccia, B. Monserrat and C. Attaccalite, Phys. Rev. B '''99''', 081109(R) (2019)</ref><ref name='fprl'>[https://arxiv.org/abs/1810.08976 Exciton-Phonon Coupling in the Ultraviolet Absorption and Emission Spectra of Bulk Hexagonal Boron Nitride], F. Paleari et al. PRL '''122''', 187401(2019) </ref><ref name=strain>[https://arxiv.org/abs/2112.12417 Excitons under strain: light absorption and emission in strained hexagonal boron nitride],P. Lechifflart, F. Paleari and C. Attaccalite (2022). </ref><ref name='mprl'>[https://arxiv.org/abs/2111.03518 Phonon-Assisted Luminescence in Defect Centers from Many-Body Perturbation Theory], F. Libbi, P. M. M. C. de Melo, Z. Zanolli, M. J. Verstraete, and N. Marzari, Phys. Rev. Lett. '''128''', 167401 (2022) </ref>.
In this tutorial we will show how to use [https://github.com/yambo-code/yambopy YamboPy] in combination with Yambo to perform these kind of calculations. For a general review of finite displacements calculations of vibrational properties we recommend this article<ref>Electron-phonon coupling from finite differences, Bartomeu Monserrat, J. Phys.: Condens. Matter '''30''' 083001(2018)</ref>. 
Notice that calculation of phonon-assisted luminescence by finite atomic displacements are quite involved and expensive from a computational point of view, for this reason they are limited to relatively simple systems. 
In this tutorial we will follow the approach of Ref.<ref name=fprl></ref>, and calculate luminescence using the Roosbroeck–Shockley (RS) relation applied to the excitonic case. 
If you use this approach to calculate luminescence or phonon-assisted absorption please cite references <ref name=cann></ref> and <ref name=fprl></ref>. 
'''Important:''' This tutorial requires that you are familiar with Bethe-Salpeter calculations with Yambo, and phonons with Quantum-Espresso, and YamboPy. If it is not the case please study the corresponding tutorial before continuing.

== Locate the lowest indirect exciton (or indirect band gap) ==

In this tutorial as example we will consider the hexagonal-BN, that is an indirect insulator. The first step to calculate the phonon-assisted luminescence is to locate the lowest indirect exciton. The cleanest strategy to identify the lowest excitation would be to calculate the excitons along the whole Irreducible Brillouin zone(IBZ), then interpolate them as explained in [http://www.yambo-code.org/wiki/index.php?title=How_to_analyse_excitons#Interpolate_exciton_dispersion_(only_in_Yambo_5.x) Interpolate exciton dispersion] and find the '''q''' point that correspond to the minimum. 
A simpler approach is to use the GW (or Kohn-Sham) band structure to find the momentum responsible for the indirect band gap, in the major part of cases this correspond to the lowest exciton. 
In the h-BN case we have:
[[File:Bn band structure.png|1000px|center| Yambo tutorial image]]

The indirect gap is between a point close to '''K''' and the '''M''' point, see panel (a) of the figure above. We will approximate the momentum responsible for the indirect emission with '''q=K-M'''. Notice that then we should build supercells that contain this moment, so it is good to approximate it with fractions of integers that are not too large so as not to have giant supercells. In this case we get '''q=(1/3,-1/6, 0)'''. 
Notice that if you have a direct band gap material, the momentum '''q''' corresponding to the lowest transition is '''q=0''' and you do not need to construct special supercells to sample the corresponding vibration, it will be sufficient to displace atoms in the primitive cell.

== Calculate phonons at momentum '''q''' of the indirect transition ==

Now we need to calculate all atomic vibrations compatible with the momentum of the lowest indirect transition. In order to do so we first calculate the phonons with QE, here the inputs [http://www.yambo-code.org/educational/tutorials/files/hBN_phonons.tgz hBN_phonons.tgz], then we calculate the force constants with q2r.x and now we can interpolate phonons along all the BZ as shown in panel (b) of Figure 1. Then we use matdyn.x to calculate phonons at momentum '''q=K-M = (1/3, -1/6, 0)''':

&input
asr='crystal',
flfrc='bn.fc',
flfrq='bn.freq',
deltaE=1.d0,
q_in_cryst_coord=.true.
fleig='bn.eig'
/
2
0 0 0
0.3333333333333333 -0.16666666666666666666 0.0

The obtained phonon modes correspond phonons at the vertical red line in panel(b) of Fig. 1. They are written in the ''matdyn.modes'' file with the corresponding eigenvectors and will be later used to displace atoms along the phonon modes.
If you open the 'matdyn.modes' file you will find phonon modes for the '''q=0''' and '''q=(1/3, -1/6, 0)''' that in cartesian coordinates is (1/3,0,0). We copy the modes corresponding to the second '''q=(1/3, -1/6, 0)''' in a separate file, called ''Q2.modes'' that will look like:

diagonalizing the dynamical matrix ...

q = 0.3333 0.0000 0.0000
**************************************************************************
freq ( 1) = 5.194139 [THz] = 173.257815 [cm-1]
( 0.000000 -0.000000 -0.000000 -0.000000 -0.274372 -0.475224 )
( -0.000000 0.000000 -0.000000 -0.000000 -0.445962 0.000001 )
( -0.000000 -0.000000 -0.000000 -0.000000 -0.222981 -0.386214 )
( 0.000000 0.000000 -0.000000 0.000000 -0.548742 0.000000 )
...........................

== Generate the supercell==

Now we use YamboPy to generate the supercell corresponding to the momentum '''q=(1/3, -1/6, 0)'''. A tutorial is already present in the file ''tutorial/supercell/generate_supercells.py''.
You can edit the file to specify the momentum '''q''', the name of the SCF initial file (uc_filnm) and the name of the phonon mode file (mode_file).

.....
uc_filnm = 'hBN.scf.in' #pw input
Q = [[1,-1,0],[3,6,0]] #qpoint to fold in nondiagonal supercell: Q=(1/3,-1/6,0)
kpoints = [12,12,4] #NB: these fractional crystal coordinates must exactly divide the kpoint mesh!
modes_file = 'Q2.modes' #file with phonon eigenmodes
.....

run the script to generate the supercell without displace atoms ''python generate_supercells.py -N''. The script will generate a non-diagonal supercell<ref>[https://arxiv.org/abs/1510.04418 Lattice dynamics and electron-phonon coupling calculations using nondiagonal supercells], J. H. Lloyd-Williams and B. Monserrat, Phys. Rev. B '''92''', 184301(2015) </ref> in the file "sc_nondiagonal.bn.scf' that is shown below:

[[File:Supercell.png|500px|center| Supercell for the momentum q=(1/3, -1/6, 0) and the corresponding reciprocal space]]

The YamboPy script automatically reduce the number of k-point in such a way to have the same sampling of the corresponding primitive cell, however we advice you to check if the new k-point grid is what you expect, and modify it in case of incorrect grid. 
In the same way you can generate the input file for non-self-consistent calculation in the supercell.

==Calculate the optical spectra in the supercell==

Now you have to calculate the optical spectra and exciton in the new supercell. Do not forget to increase the number of conduction bands in both the screening the Bethe-Salpeter equation. For example in the supercell shown above there are 24 atoms, 6 times more than in the original cell, so you will have to multiply the number of bands by 6, if in the base cell you have used 40 bands here you will have to put 240 to have the same convergence parameters. Then print exciton according to their energy with the command ypp -e s. 

'''N. B. ''': in order to not recalculate W for each atomic displacement you can perform calculations without symmetries in the pristine supercell and then copy the corresponding dielectric constant in the calculations with displaced atoms, this is usually a safe approximation described in Ref.<ref name="cann"></ref>. 

After solving the BSE in the supecell you will find the excitons that were already present in the primary cell (within the numerical noise) and new lower energy excitons, with zero dipole matrix elements, that correspond to the indirect excitons that have now been mapped at '''q=0''', something like:

# Maximum Residual Value = 0.38659E+01
#
# E [ev] Strength Index
#
5.7165937 0.47290171E-11 1.0000000 <-- Indirect exciton mapped at q=0
5.7166047 0.58239158E-09 2.0000000 <-- Indirect exciton mapped at q=0
5.7367945 0.31830746E-09 3.0000000 <-- Indirect exciton mapped at q=0
5.7368207 0.33674005E-10 4.0000000 <-- Indirect exciton mapped at q=0
5.80542803 0.100763231E-6 5.00000000 <-- Direct exciton
5.80598068 0.384604526E-8 6.00000000 <-- Direct exciton
5.88944197 1.00000000 7.00000000 <-- Direct exciton
5.88967037 0.102969212E-1 8.00000000 <-- Direct exciton
................

This result shows that indirect exciton without the coupling with phonon do not contribute the optical response. 

==Displace atoms and calculate dipoles derivatives and luminescence==
Now you can generate the displaced supercell for all the phonon modes compatible with the indirect band gap '''q''' point by doing: ''python generate_supercells.py -S''. The script will produce a set of input files for each of the 12 phonon modes at '''q=(1/3,-1/6,0)'''. You have to perform BSE calculation for each one of these modes and get exciton dipoles. 
For each phonon mode the atoms are displace along the corresponding eigenvector read from the "Q2.modes" file multiplied for 0.1 Borh (the variable Temp=0.1), see the following line in "generate_supercell.py":

...
sc.displace(eivs,nd_atom_positions,Temp=0.1) # Generate list of displaced supercells as PwIn objects called 'modes_qe'
...

this is usually a good value to get correct derivatives. 
Once you get dipoles for a give phonon mode you can calculate the second order dipole derivatives using a finite different formula and the fact the first order derivatives are zero by construction. 
Now you have all the information necessary to construct the emission spectra. Using the formula, derived in Ref.<ref name=fprl></ref> you obtain the phonon-assisted luminescence as:
[[File:Light emission.png|center|Light emission formula]]

in this formula:
* nr(ω) is the refractive index, that usually can be approximated with a constant in the luminescence energy range.
* ES'''q''' are the excitonic energies in the undistorted supercell, including indirect excitons
* B(ES,Texc) is the Bose occupation of excitons, that can be approximated with a Boltzman function, see Ref.<ref name=fprl></ref><ref name=strain></ref>
* Texc is the excitonic temperature that should be obtained from experimental measurements or used as a parameters, for a discussion see Ref.<ref name=fprl></ref>
* d2|TS|2/dR2λ'''q''' are the derivatives of the exciton dipole matrix elements for a given phonon mode
* Ωλ,'''q''' are the phonon frequencies and the sum on '''q''' is on the different supercells obtained from phonon modes responsible for the luminescence, in this tutorial there is just a single '''q''' point.
* η is the exciton line-width that is given as input parameter to define the broadening of the spectra

In the h-BN example shown here, the Texc was obtained from experimental measurements, only the lowest two excitons were included in the sum and the Bose function was replaced by a Botlzman distribution with the energy zero was fixed at the lowest excitonic energy. An appropriate scissor operator independent of the atomic displacements was added in the calculation. The final result is shown in the figure above and compared with experimental measurement at 10 Kelvin
[[File:Lum spectra.png|center|800px|Luminescence spectra]]

'''Some additional considerations on luminescence calculations'''

* In general due to the '''Bose''' (or '''Boltzmann''') factor only the lower excitons (often only the lowest one) contribute to the luminescence signal
* The lowest excitons are '''degenerate''' all of them have to be included in the sum, as in the h-BN example
* In the present formulation of luminescence spectra excitons are treated as Bosons (or Boltzman particles), in principle it is possible to define a '''excitonic occupation''' starting from a quasi-equilibrium distribution of excited electron-hole, but this requires a rotation of the electron-hole occupations in the excitonic base for a discussion see Ref.<ref name='cann'></ref><ref name='mprl'></ref>
* In order to '''speed up calculations''' you can perform all calculation without symmetries and copy the dielectric constant (''ndb.em1s*'' or ''ndb.pp*'' files) from the undistorted supercell to the displaced ones. In this case if you do not need to calculate GW correction, for example using a scissor, you can also reduce the number of bands in the the non-self-consistent calculation to the one required by the BSE, for a discussion see Ref<ref name=strain></ref>
* Notice that the '''dipole matrix elements''' that will enter in optical absorption and luminescence depend form the direction of the electric field. In this example show before the luminescence has been averaged on ''x'' and ''y'' directions because we were interested in the in-place luminescence. In general is it possible to select specific directions to reproduce particular experimental setups<ref name=cann></ref><ref name=fprl></ref>. 

== References ==

Yambopy tutorial: band structures

2022-04-08T07:37:58Z

Pierre:

This tutorial will show you how to visualise wave-function and band-structure related information following a DFT Quantum Espresso calculation using the "qepy" module of Yambopy.

'''Concerning the ICTP school: recall that for the yambopy tutorials you have to load the yambo, quantum-espresso and anaconda3 modules:'''
spack load yambo
spack load quantum-espresso
spack load anaconda3
'''You can copy the tutorial databases as explained in the virtual machine setup page: '''
cp -r /media/ictpuser/smr3694/ictptutor/YAMBOPY_TUTORIALS ~/

The full tutorial, including the Quantum espresso and Yambo databases that we will read, can be downloaded and extracted from the yambo website:
$wget www.yambo-code.org/educational/tutorials/files/databases_qepy.tar
$tar -xvf databases_qepy.tar
$cd databases_qepy

==Tutorial 1. BN (semiconductor). Band structure==

First enter in the folder
cd databases_qepy/bn-semiconductor
and have a look to the
vim plot-qe-bands.py
The qepy classes are useful both to execute Quantum Espresso and to analyze the results. The full qepy library is imported by simply doing:

from qepy import *

===Plot Band structure===

The qepy class '''PwXML''' reads the data file generated by Quantum Espresso and post-processes the data. The class is instanced by doing:

xml = PwXML(prefix='bn', path='bands')

The variable prefix corresponds to the same variable of the QE input. The folder location is indicated by variable path. In order to plot the bands, we also define the k-points path (in crystal coordinates) using the function Path:

npoints = 50
path_kpoints = Path([ [[0.0, 0.0, 0.0],'$\Gamma$'],
[[0.5, 0.0, 0.0],'M'],
[[1./3,1./3,0.0],'K'],
[[0.0, 0.0, 0.0],'$\Gamma$']], [int(npoints*2),int(npoints),int(sqrt(5)*npoints)])

It is worth to note that the path should coincide with the selected path for the QE band calculations.

In order to show the plot we call the '''plot_eigen''' method of the '''PwXML''' class:

xml.plot_eigen(path_kpoints)

This function will automatically plot the bands as shown below running the script:

python plot-qe-bands.py

[[File:Bands BN 1.png| 400 px | center |Band structure of BN calculated at the level of DFT-LDA]]

Alternatively, we can use the function '''plot_eigen_ax'''. This functions requires as input a matplotlib '''figure''' object with given axes, as you will see in the next example.

===Plot atomic orbital projected Band structure===

In addition to the band structure, useful information regarding the atomic orbital nature of the electronic wave functions can be displayed using the class '''ProjwfcXML'''.
In order to make quantum espresso calculate the relevant data, we need to use the QE executable '''projwfc.x''', which will create the file '''atomic_proj.xml'''. This executable projects the Kohn-Sham wavefunctions onto orthogonalized atomic orbitals, among others functionalities. The orbital indexing and ordering are explained in the input documentation of the projwfc.x function which you are invited to check (https://www.quantum-espresso.org/Doc/INPUT_PROJWFC.html#idm94). We can run '''projwf.x''' directly from the python script using the qepy class '''ProjwfcIn''':

proj = ProjwfcIn(prefix='bn')
proj.run(folder='bands')

Be aware that this can take a while in a large system with many k-points. As in the class '''PwXML''', we then define the path_kpoints and also the selected atomic orbitals to project our bands. We have chosen to represent the projection weight onto the nitrogen (1) and boron (2) orbitals, which according to the projection orderings is given by

atom_1 = list(range(16))
atom_2 = list(range(16,32))

The full list of orbitals is written in the file '''bands/prowfc.log'''. We have also defined the figure box

import matplotlib.pyplot as plt
fig = plt.figure(figsize=(5,7))
ax = fig.add_axes( [ 0.12, 0.10, 0.70, 0.80 ])

The class '''ProjwfcXML''' then runs as in this example:

band = ProjwfcXML(prefix='bn',path='bands',qe_version='7.0')
band.plot_eigen(ax,path_kpoints=path_kpoints,cmap='viridis',selected_orbitals=atom_1,selected_orbitals_2=atom_2)

We can run now the file:

python plot-qe-orbitals.py

[[File:Bands AOP BN 1.png| 400 px | center | Atomic orbital projected band structure of monolayer BN]]

We have chosen the colormap 'viridis' (variable cmap). You see that the colormap goes from maximum '''selected_orbitals''' content (in this case nitrogen) to the maximum '''selected_orbitals_2''' content (in this case boron).
The colormap can be represented in many ways and qepy does not include a particular function for this. An example is:

import matplotlib as mpl
cmap=plt.get_cmap('viridis')
bx = fig.add_axes( [ 0.88, 0.10, 0.05, 0.80 ])
norm = mpl.colors.Normalize(vmin=0.,vmax=1.)
cb1 = mpl.colorbar.ColorbarBase(bx, cmap=cmap, norm=norm,orientation='vertical',ticks=[0,1])
cb1.set_ticklabels(['B', 'N'])

We can also plot the band orbital charater with size variation instead of a color scale. In this case we have to pass only the variable selected_orbitals (see the next tutorial).

==Tutorial 2. Iron. Ferromagnetic metalic material==

In the case of spin-polarized calculations we can plot the spin up and down band structures. We have included in the tutorial a small workflow example to run quantum espresso calculations from scratch. This is done using the classes Tasks and Flows developed in yambopy. You can check the file flow-iron.py for an example.
Yambopy includes basic predefined workflows to run the common QE and Yambo calculations. In this case we are using the class '''PwBandsTasks'''.

python flow-iron.py

In order to plot the spin-polarized bands. After doing all the calculations from scratch with the flows (flow-iron.py), we can make the band plot by running the script:

python plot-qe-bands.py

The class PwXML automatically detects the spin polarized case (nspin=2 in the QE input file). The spin up channel is displayed with red and the spin down channel with blue. In the case of iron we have selected this k-point path:

npoints = 50
path_kpoints = Path([ [[0.0, 0.0, 0.0 ],'G'],
[[0.0, 0.0, 1.0 ],'H'],
[[1./2,0.0,1./2.],'N'],
[[0.0, 0.0, 0.0 ],'G'],
[[1./2, 1./2, 1./2 ],'P'],
[[1./2,0.0,1./2. ],'N']], [npoints,npoints,npoints,npoints,npoints])

xml = PwXML(prefix='pw',path='bands/t0')
xml.plot_eigen(path_kpoints)

[[File:Figure 3-iron-bands.png| 600 px | center |Spin polarized band structure of iron calculated by DFT]]

The analysis of the projected atomic orbitals is also implemented. In this case the results are more cumbersome because the projection is separated in spin up and down channels.
Let us look first at the file '''plot-qe-orbitals-size'''.

'''NB: If you generated the quantum espresso databases using the flows instead of relying on the precomputed databases, you also need to rerun the quantum espresso executable projwfc.x to recompute the orbital projections. In this case, please uncomment the following lines in the script:'''
#proj = ProjwfcIn(prefix='pw')
#proj.run(folder='bands/t0')

Now, we can use the dot size as a function of the weight of the orbitals
atom_s = [8]
atom_p = [0,1,2]
atom_d = [3,4,5,6,7]

and the plots are done with
band = ProjwfcXML(prefix='pw',path='bands/t0',qe_version='7.0')
band.plot_eigen(ax,path_kpoints=path_kpoints,selected_orbitals=atom_s,color='pink',color_2='black')
band.plot_eigen(ax,path_kpoints=path_kpoints,selected_orbitals=atom_p,color='green',color_2='orange')
band.plot_eigen(ax,path_kpoints=path_kpoints,selected_orbitals=atom_d,,color='red',color_2='blue')

As an example, we can select just the ''d'' orbitals by commenting the first two plots and then running the file:

plot-qe-orbitals-size.py

[[File:Figure 4-iron-bands-size-d-orbitals.png|400px|center|Iron band structure. Size is proportional to the weights of the projection on atomic d-orbitals. Red (blue) is up (down) spin polarization.]]

From the plot is clear that ''d'' orbitals are mainly localized around the Fermi energy. The plot above is in red and blue, while the default choice in your script should be pink and black. You can experiment with the colors and other ''matplotlib'' plot options and also plot the other orbital types.

Another option is to plot the orbital composition as a colormap running the file:

plot-qe-orbitals-colormap.py

[[File:Figure 5-iron-bands-colormap.png|400px|center]]

Here we have added the p and d orbitals in the analysis:

atom_p = [0,1,2]
atom_d = [3,4,5,6,7]
band = ProjwfcXML(prefix='pw',path='bands/t0',qe_version='7.0')
band.plot_eigen(ax,path_kpoints=path_kpoints,cmap='viridis',cmap2='rainbow',selected_orbitals=atom_p,selected_orbitals_2=atom_d)

The colormap bar is added in the same way as in Tutorial 1 (see script), but this time we have a different colormap for each spin polarisation.

==Tutorial 3==

Yambopy can be used either to run Yambo and QE calculations, or to analyse the results of QE and Yambo by dealing with their generated databases. This is done with a variety of classes included in the qepy (for QE) or yambopy (for Yambo) modules.
In the case of Yambo GW quasi-particle calculations, we can use the yambopy class '''YamboQPDB''' to read the database produced by the simulation. We can use this to find the scissor operator, plot the GW bands and to interpolate the GW bands on a smoother k-path. The example runs by typing:

python plot-qp.py

See below for an explanation of the tutorial. As usual, we can import the qepy and yambopy libraries:

from qepy import *
from yambopy import *
import matplotlib.pyplot as plt

We define the k-points path.
npoints = 10
path = Path([ [[ 0.0, 0.0, 0.0],'$\Gamma$'],
[[ 0.5, 0.0, 0.0],'M'],
[[1./3.,1./3., 0.0],'K'],
[[ 0.0, 0.0, 0.0],'$\Gamma$']], [int(npoints*2),int(npoints),int(sqrt(5)*npoints)] )

Importantly, the number of points is a free choice. We can increase the variable npoints as we wish, it just means that the interpolation step will take more time. In order to analyse GW results we need to have the file related to the basic data of our Yambo calculation ('''SAVE/ns.db1''') and the netcdf file with the quasi-particle results ('''ndb.QP'''). We load the data calling the respective classes:

# Read Lattice information from SAVE
lat = YamboSaveDB.from_db_file(folder='SAVE',filename='ns.db1')
# Read QP database
ydb = YamboQPDB.from_db(filename='ndb.QP',folder='qp-gw')
(in the yambopy module, each class is specialised to read a specific Yambo database)

The first option is to plot the energies and calculate the ideal Kohn-Sham to GW scissor operator. We need to select the index of the top valence band:

n_top_vb = 4
ydb.plot_scissor_ax(ax,n_top_vb)

Yambopy displays the fitting and also the data of the slope of each fitting. Notice that this is also a test if the GW calculations are running well. '''If the dependence is not linear you should double-check your results!'''

[[File:Figure 6-slope-scissor.png|400px|center]]

In this case the slope is:

valence bands:
slope: 1.0515886598785766
conduction bands:
slope: 1.026524081134514
scissor list (shift,c,v) [eV,adim,adim]: [1.8985204833551723, 1.026524081134514, 1.0515886598785766]

In addition to the scissor operator, we can plot the GW (and DFT) band structure along the path. The first choice would be to plot the actual GW calculations, without interpolation, to check that the results are meaningful (or not). This plot is independent of the number of k-points ('''npoints''') that we want to put in the interpolation. The class YamboQPDB finds the calculated points that belong to the path and plots them. Be aware that if we use coarse grids the class would not find any point and the function will not work.

ks_bs_0, qp_bs_0 = ydb.get_bs_path(lat,path)
ks_bs_0.plot_ax(ax,legend=True,color_bands='r',label='KS')
qp_bs_0.plot_ax(ax,legend=True,color_bands='b',label='QP-GW')

[[File:Figure 7-GW-band-structure-non-interpolated.png|400px|center]]

The interpolation of the DFT and GW band structures looks similar:

ks_bs, qp_bs = ydb.interpolate(lat,path,what='QP+KS',lpratio=20)
ks_bs.plot_ax(ax,legend=True,color_bands='r',label='KS')
qp_bs.plot_ax(ax,legend=True,color_bands='b',label='QP-GW')

The '''lpratio''' can be increased if the interpolation does not work as well as intended. The interpolation is the same one implemented in abipy and currently requires abipy installed in order to work.

[[File:Figure 8-GW-band-structure-interpolated.png|400px|center]]

Finally, we can compare the calculated GW eigenvalues with the interpolation.

[[File:Figure 8-GW-band-structure-comparison.png|400px|center]]

==Links==
* Back to [[ICTP 2022#Tutorials]]

Yambopy tutorial: band structures

2022-04-07T21:52:02Z

Pierre: /* Plot atomic orbital projected Band structure */

This tutorial will show you how to visualise wave-function and band-structure related information following a DFT Quantum Espresso calculation using the "qepy" module of Yambopy.

'''Concerning the ICTP school: recall that for the yambopy tutorials you have to load the yambo, quantum-espresso and anaconda3 modules:'''
spack load yambo
spack load quantum-espressp
spack load anaconda3
'''You can copy the tutorial databases as explained in the virtual machine setup page: '''
cp -r /media/ictpuser/smr3694/ictptutor/YAMBOPY_TUTORIALS ~/

The full tutorial, including the Quantum espresso and Yambo databases that we will read, can be downloaded and extracted from the yambo website:
$wget www.yambo-code.org/educational/tutorials/files/databases_qepy.tar
$tar -xvf databases_qepy.tar
$cd databases_qepy

==Tutorial 1. BN (semiconductor). Band structure==

The qepy classes are useful both to execute Quantum Espresso and to analyze the results. The full qepy library is imported by simply doing:

from qepy import *

===Plot Band structure===

The qepy class '''PwXML''' reads the data file generated by Quantum Espresso and post-processes the data. The class is instanced by doing:

xml = PwXML(prefix='bn', path='bands')

The variable prefix corresponds to the same variable of the QE input. The folder location is indicated by variable path. In order to plot the bands, we also define the k-points path (in crystal coordinates) using the function Path:

npoints = 50
path_kpoints = Path([ [[0.0, 0.0, 0.0],'$\Gamma$'],
[[0.5, 0.0, 0.0],'M'],
[[1./3,1./3,0.0],'K'],
[[0.0, 0.0, 0.0],'$\Gamma$']], [int(npoints*2),int(npoints),int(sqrt(5)*npoints)])

It is worth to note that the path should coincide with the selected path for the QE band calculations.

In order to show the plot we call the '''plot_eigen''' method of the '''PwXML''' class:

xml.plot_eigen(path_kpoints)

This function will automatically plot the bands as shown below running the script:

python plot-qe-bands.py

[[File:Bands BN 1.png| 400 px | center |Band structure of BN calculated at the level of DFT-LDA]]

Alternatively, we can use the function '''plot_eigen_ax'''. This functions requires as input a matplotlib '''figure''' object with given axes, as you will see in the next example.

===Plot atomic orbital projected Band structure===

In addition to the band structure, useful information regarding the atomic orbital nature of the electronic wave functions can be displayed using the class '''ProjwfcXML'''.
In order to make quantum espresso calculate the relevant data, we need to use the QE executable '''projwfc.x''', which will create the file '''atomic_proj.xml'''. This executable projects the Kohn-Sham wavefunctions onto orthogonalized atomic orbitals, among others functionalities. The orbital indexing and ordering are explained in the input documentation of the projwfc.x function which you are invited to check (https://www.quantum-espresso.org/Doc/INPUT_PROJWFC.html#idm94). We can run '''projwf.x''' directly from the python script using the qepy class '''ProjwfcIn''':

proj = ProjwfcIn(prefix='bn')
proj.run(folder='bands')

Be aware that this can take a while in a large system with many k-points. As in the class '''PwXML''', we then define the path_kpoints and also the selected atomic orbitals to project our bands. We have chosen to represent the projection weight onto the nitrogen (1) and boron (2) orbitals, which according to the projection orderings is given by

atom_1 = list(range(16))
atom_2 = list(range(16,32))

The full list of orbitals is written in the file '''bands/prowfc.log'''. We have also defined the figure box

import matplotlib.pyplot as plt
fig = plt.figure(figsize=(5,7))
ax = fig.add_axes( [ 0.12, 0.10, 0.70, 0.80 ])

The class '''ProjwfcXML''' then runs as in this example:

band = ProjwfcXML(prefix='bn',path='bands',qe_version='7.0')
band.plot_eigen(ax,path_kpoints=path_kpoints,cmap='viridis',selected_orbitals=atom_1,selected_orbitals_2=atom_2)

We can run now the file:

python plot-qe-orbitals.py

[[File:Bands AOP BN 1.png| 400 px | center | Atomic orbital projected band structure of monolayer BN]]

We have chosen the colormap 'viridis' (variable cmap). You see that the colormap goes from maximum '''selected_orbitals''' content (in this case nitrogen) to the maximum '''selected_orbitals_2''' content (in this case boron).
The colormap can be represented in many ways and qepy does not include a particular function for this. An example is:

import matplotlib as mpl
cmap=plt.get_cmap('viridis')
bx = fig.add_axes( [ 0.88, 0.10, 0.05, 0.80 ])
norm = mpl.colors.Normalize(vmin=0.,vmax=1.)
cb1 = mpl.colorbar.ColorbarBase(bx, cmap=cmap, norm=norm,orientation='vertical',ticks=[0,1])
cb1.set_ticklabels(['B', 'N'])

We can also plot the band orbital charater with size variation instead of a color scale. In this case we have to pass only the variable selected_orbitals (see the next tutorial).

==Tutorial 2. Iron. Ferromagnetic metalic material==

In the case of spin-polarized calculations we can plot the spin up and down band structures. We have included in the tutorial a small workflow example to run quantum espresso calculations from scratch. This is done using the classes Tasks and Flows developed in yambopy. You can check the file flow-iron.py for an example.
Yambopy includes basic predefined workflows to run the common QE and Yambo calculations. In this case we are using the class '''PwBandsTasks'''.

python flow-iron.py

In order to plot the spin-polarized bands. After doing all the calculations from scratch with the flows (flow-iron.py), we can make the band plot by running the script:

python plot-qe-bands.py

The class PwXML automatically detects the spin polarized case (nspin=2 in the QE input file). The spin up channel is displayed with red and the spin down channel with blue. In the case of iron we have selected this k-point path:

npoints = 50
path_kpoints = Path([ [[0.0, 0.0, 0.0 ],'G'],
[[0.0, 0.0, 1.0 ],'H'],
[[1./2,0.0,1./2.],'N'],
[[0.0, 0.0, 0.0 ],'G'],
[[1./2, 1./2, 1./2 ],'P'],
[[1./2,0.0,1./2. ],'N']], [npoints,npoints,npoints,npoints,npoints])

xml = PwXML(prefix='pw',path='bands/t0')
xml.plot_eigen(path_kpoints)

[[File:Figure 3-iron-bands.png| 600 px | center |Spin polarized band structure of iron calculated by DFT]]

The analysis of the projected atomic orbitals is also implemented. In this case the results are more cumbersome because the projection is separated in spin up and down channels. We can use the dot size as a function of the weight of the orbitals, as we have done in this example by selecting exclusively the ''d'' orbitals:

atom_d = [3,4,5,6,7]
band = ProjwfcXML(prefix='pw',path='bands/t0',qe_version='7.0')
band.plot_eigen(ax,path_kpoints=path_kpoints,selected_orbitals=atom_d,color='red',color_2='blue')

If we run the file:

plot-qe-orbitals-size.py

[[File:Figure 4-iron-bands-size-d-orbitals.png|400px|center|Iron band structure. Size is proportional to the weights of the projection on atomic d-orbitals. Red (blue) is up (down) spin polarization.]]

From the plot is clear that ''d'' orbitals are mainly localized around the Fermi energy.

Another option is to plot the orbital composition as a colormap running the file:

plot-qe-orbitals-colormap.py

[[File:Figure 5-iron-bands-colormap.png|400px|center]]

Here we have added the p and d orbitals in the analysis:

atom_p = [0,1,2]
atom_d = [3,4,5,6,7]
band = ProjwfcXML(prefix='pw',path='bands/t0',qe_version='7.0')
band.plot_eigen(ax,path_kpoints=path_kpoints,cmap='viridis',cmap2='rainbow',selected_orbitals=atom_p,selected_orbitals_2=atom_d)

The colormap bar is added in the same way as in Tutorial 1 (see script), but this time we have a different colormap for each spin polarisation.

==Tutorial 3==

Yambopy can be used either to run Yambo and QE calculations, or to analyse the results of QE and Yambo by dealing with their generated databases. This is done with a variety of classes included in the qepy (for QE) or yambopy (for Yambo) modules.
In the case of Yambo GW quasi-particle calculations, we can use the yambopy class '''YamboQPDB''' to read the database produced by the simulation. We can use this to find the scissor operator, plot the GW bands and to interpolate the GW bands on a smoother k-path. The example runs by typing:

python plot-qp.py

See below for an explanation of the tutorial. As usual, we can import the qepy and yambopy libraries:

from qepy import *
from yambopy import *
import matplotlib.pyplot as plt

We define the k-points path.
npoints = 10
path = Path([ [[ 0.0, 0.0, 0.0],'$\Gamma$'],
[[ 0.5, 0.0, 0.0],'M'],
[[1./3.,1./3., 0.0],'K'],
[[ 0.0, 0.0, 0.0],'$\Gamma$']], [int(npoints*2),int(npoints),int(sqrt(5)*npoints)] )

Importantly, the number of points is a free choice. We can increase the variable npoints as we wish, it just means that the interpolation step will take more time. In order to analyse GW results we need to have the file related to the basic data of our Yambo calculation ('''SAVE/ns.db1''') and the netcdf file with the quasi-particle results ('''ndb.QP'''). We load the data calling the respective classes:

# Read Lattice information from SAVE
lat = YamboSaveDB.from_db_file(folder='SAVE',filename='ns.db1')
# Read QP database
ydb = YamboQPDB.from_db(filename='ndb.QP',folder='qp-gw')
(in the yambopy module, each class is specialised to read a specific Yambo database)

The first option is to plot the energies and calculate the ideal Kohn-Sham to GW scissor operator. We need to select the index of the top valence band:

n_top_vb = 4
ydb.plot_scissor_ax(ax,n_top_vb)

Yambopy displays the fitting and also the data of the slope of each fitting. Notice that this is also a test if the GW calculations are running well. '''If the dependence is not linear you should double-check your results!'''

[[File:Figure 6-slope-scissor.png|400px|center]]

In this case the slope is:

valence bands:
slope: 1.0515886598785766
conduction bands:
slope: 1.026524081134514
scissor list (shift,c,v) [eV,adim,adim]: [1.8985204833551723, 1.026524081134514, 1.0515886598785766]

In addition to the scissor operator, we can plot the GW (and DFT) band structure along the path. The first choice would be to plot the actual GW calculations, without interpolation, to check that the results are meaningful (or not). This plot is independent of the number of k-points ('''npoints''') that we want to put in the interpolation. The class YamboQPDB finds the calculated points that belong to the path and plots them. Be aware that if we use coarse grids the class would not find any point and the function will not work.

ks_bs_0, qp_bs_0 = ydb.get_bs_path(lat,path)
ks_bs_0.plot_ax(ax,legend=True,color_bands='r',label='KS')
qp_bs_0.plot_ax(ax,legend=True,color_bands='b',label='QP-GW')

[[File:Figure 7-GW-band-structure-non-interpolated.png|400px|center]]

The interpolation of the DFT and GW band structures looks similar:

ks_bs, qp_bs = ydb.interpolate(lat,path,what='QP+KS',lpratio=20)
ks_bs.plot_ax(ax,legend=True,color_bands='r',label='KS')
qp_bs.plot_ax(ax,legend=True,color_bands='b',label='QP-GW')

The '''lpratio''' can be increased if the interpolation does not work as well as intended. The interpolation is the same one implemented in abipy and currently requires abipy installed in order to work.

[[File:Figure 8-GW-band-structure-interpolated.png|400px|center]]

Finally, we can compare the calculated GW eigenvalues with the interpolation.

[[File:Figure 8-GW-band-structure-comparison.png|400px|center]]

==Links==
* Back to [[ICTP 2022#Tutorials]]

Optical properties at finite temperature

2022-04-07T11:07:03Z

Pierre: /* Result analysis */

[[File:High low.png|thumb | 200px|Yambo tutorial image]]

In this tutorial, we will show you how to calculate optical properties including thermal effects due to the electron-phonon coupling. 
This tutorial assumes that you have completed all the steps from the previous tutorial on [[Electron Phonon Coupling]]. 
The tutorial is dived into different steps first we will calculate absorption at independent particle approximation, then we will include excitonic effects,
and finally, we will show how to analyze the data. 

'''Notice that there are two ways to include finite temperature effects in the BSE''':
* The first is to use a full Bethe-Salpeter equation with all the coupling terms in such a way to treat complex quasi-particles that cannot be handled in standard (Hermitian) BSE, this is explained in section [[https://www.yambo-code.org/wiki/index.php?title=Optical_properties_at_finite_temperature#Bethe-Salpeter_at_finite_temperature_.28diagonalization_solver.29 Bethe-Salpeter at finite temperature (diagonalization solver)]]. This approach has the advantage that one can obtain excitonic eigenvalues and eigenfunctions at finite temperature and not only the optical spectrum.

* The second method is use of a standard Hermitinan (resonant) BSE plus the inversion solver <ref name='dbgrid'>D. Kammerlander et al. [https://arxiv.org/abs/1209.1509 Phys. Rev. B '''86''', 125203] </ref>. In this case, a Lo two-particle Green's function that contains the complex quasi-particles is used in the solver to get a finite temperature spectrum. This approach has the advantage that can be used with the double-grid but it does not give access to eigenvalues and eigenfunctions of the BSE. This approach is described in the section: [[https://www.yambo-code.org/wiki/index.php?title=Optical_properties_at_finite_temperature#Bethe-Salpeter_at_finite_temperature_.28inversion_solver.29 Bethe-Salpeter at finite temperature (inversion solver)]]

== Absorption at finite temperature ==

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 <code>yambo_ph -o c -V qp -F yambo.in_300K</code>:

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= 1500 # [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. Run your calculations with the command: <code>yambo_ph -F yambo.in_300K -J T300</code> and do the same for the 0K case. 
 You can plot the results with gnuplot using the command: <code> plot 'o-T300.eps_q1_ip' w l t '300 K', 'o-T0.eps_q1_ip' w l t '0 K'</code>. Hereafter the result without and with electron-phonon coupling at two different temperatures:
[[File:Si absorption finite t.png|800px|center |Absorption of bulk silicon at finite temperature]]

The temperature effect is clearly visible in the figure.

== Bethe-Salpeter at finite temperature (diagonalization solver) ==

In this section we will calculate the Bethe-Salpeter at finite temperature using the full-BSE.

You can generate the input using the command <code>yambo_ph -X s -o b -k sex -y d -V qp -F yambo.in_BSE</code>

em1s # [R][Xs] Statically Screened Interaction
optics # [R] Linear Response optical properties
bss # [R] BSE solver
bse # [R][BSE] Bethe Salpeter Equation.
dipoles # [R] Oscillator strenghts (or dipoles)
DIP_Threads=0 # [OPENMP/X] Number of threads for dipoles
X_Threads=0 # [OPENMP/X] Number of threads for response functions
K_Threads=0 # [OPENMP/BSK] Number of threads for response functions
Chimod= "HARTREE" # [X] IP/Hartree/ALDA/LRC/PF/BSfxc
BSEmod= "coupling" # [BSE] resonant/retarded/coupling
BSKmod= "SEX" # [BSE] IP/Hartree/HF/ALDA/SEX/BSfxc
BSSmod= "d" # [BSS] (h)aydock/(d)iagonalization/(s)lepc/(i)nversion/(t)ddft`
BSENGexx= 9257 RL # [BSK] Exchange components
BSENGBlk=-1 RL # [BSK] Screened interaction block size [if -1 uses all the G-vectors of W(q,G,Gp)]
#WehCpl # [BSK] eh interaction included also in coupling
KfnQPdb= "E W < T0/ndb.QP" # [EXTQP BSK BSS] Database action
% KfnQP_E
0.500000 | 1.000000 | 1.000000 | # [EXTQP BSK BSS] E parameters (c/v) eV|adim|adim
%
BSEprop= "abs" # [BSS] Can be any among abs/jdos/kerr/magn/dich/photolum/esrt
BSEdips= "none" # [BSS] Can be "trace/none" or "xy/xz/yz" to define off-diagonal rotation plane
% BSEQptR
1 | 1 | # [BSK] Transferred momenta range
%
% BSEBands
 2 | 7 | # [BSK] Bands range
%
% BEnRange
0.000000 | 5.000000 | eV # [BSS] Energy range
%
% BDmRange
0.010000 | 0.010000 | eV # [BSS] Damping range
%
BEnSteps= 1500 # [BSS] Energy steps
% BLongDir
1.000000 | 0.000000 | 0.000000 | # [BSS] [cc] Electric Field
%
#WRbsWF # [BSS] Write to disk excitonic the WFs
% BndsRnXs
1 | 12 | # [Xs] Polarization function bands
%
NGsBlkXs= 113 RL # [Xs] Response block size
% LongDrXs
1.000000 | 0.000000 | 0.000000 | # [Xs] [cc] Electric Field
%
XTermKind= "none" # [X] X terminator ("none","BG" Bruneval-Gonze)

In this input file we asked Yambo to include finite temperature quasi-particle in the BSE with the line <code>KfnQPdb= "E W < T0/ndb.QP" </code>, notice that we also introduce an additional rigid shift of 0.5 eV to mimic the GW correction with the line <code> 0.500000 | 1.000000 | 1.000000 | \# [EXTQP BSK BSS] E parameters (c/v) eV|adim|adim</code>. You can now run the calculation with the command <code> yambo_ph -F yambo.in_BSE -J T0_BSE</code> 
In principle you can calculate the GW correction following the [https://www.yambo-code.org/wiki/index.php?title=How_to_obtain_the_quasi-particle_band_structure_of_a_bulk_material:_h-BN this tutorial] and then merge the corresponding ndb.QP database with the one of the electron-phonon coupling using the command <code> ypp -qpdb m</code>.
Notice that we changed the BSE type to "coupling" because you need the full Bethe-Salpeter to deal with complex quasi-particles. 
Now we do the same calculation at T = 300 K. Copy the previous input file <code> cp yambo.in_BSE yambo.in_BSE_300</code>. Then modify where the code reads the QP database :
KfnQPdb= "E W < T300/ndb.QP" # [EXTQP BSK BSS] Database action
Run the command with <code>yambo_ph -F yambo.in_BSE_300 -J T300_BSE</code> 
Finally, you should be able to plot the absorption spectra from the outputs in gnuplot :
<code> plot 'o-T0_BSE.eps_q1_diago_bse' w l t '0 K', 'o-T300_BSE.eps_q1_diago_bse' w l t '300 K'</code>
[[File:Si bse optics.png|800px|center |Bethe-Salpeter at finite temperature for bulk silicom]]

'''Approximations:''' 
In this calculation we have made different approximations: 
1) we did not include the renormalization of excitons due to the change of the screening potential W with temperature. Including electron-phonon coupling
in the dielectric constant by changing the line <code> XfnQPdb= "none" </code> is unfortunately not enough, for a discussion see ref. <ref>L. Adamska and P. Umari, [https://journals.aps.org/prb/abstract/10.1103/PhysRevB.103.075201 Phys. Rev. B '''103''', 075201 (2021)] </ref> 
2) The electron-phonon coupling should enter in the BSE through the exciton-phonon matrix elements and not from the single-particle self-energy, as we have done in this tutorial. The approximation used in this tutorial is valid for not too strong bound excitons, and in general it generated a finite life-time also for the lowest exciton in direct materials that is not correct. For a discussion see refs.<ref>G. Antonius, S. G. Louie [https://arxiv.org/abs/1705.04245 Phys. Rev. B '''105''', 085111 (2022)]</ref><ref name='claudio'> see Supp. Mat. of [https://journals.aps.org/prb/abstract/10.1103/PhysRevB.99.081109 Phys. Rev. B '''99''', 081109(R) (2019)]</ref><ref name='fulvio'>F. Paleari, [https://orbilu.uni.lu/handle/10993/41058 Phd Thesis (2019)]</ref> and <ref>H. Chen, D. Sangalli, and M. Bernardi [https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.125.107401 Phys. Rev. Lett. '''125''', 107401 (2020)]</ref>. 
3) we did not include scattering between exciton and phonons, therefore our results do not contain phonon-assisted absorption peaks, you can include these terms by finite difference displacements see refs. <ref name='fulvio'></ref>.

== Result analysis ==
Results of BSE at finite temperature can be analized in the same way of standard BSE. You can plot exciton binging energy as a function of the temperature, exciton wave-function and so, see [https://www.yambo-code.org/wiki/index.php?title=How_to_analyse_excitons Exciton Plot tutorial].
For example if you sort exciton with the command <code> ypp_ph -e s -J T0_BSE</code> for the case T=0K, in addition to the exciton energy now you will find its line-width in the files .
<code>o-T0_BSE.exc_qpt1_E_sorted</code> and <code>o-T0_BSE.exc_qpt1_I_sorted</code>.
#
# Maximum Residual Value = .77136E+05
#
# E [ev] Strength Index W [meV]
#
...............
2.82974172 0.388171151E-1 335.000000 11.9707289
2.82976580 0.110338017 334.000000 11.9676380
2.82979798 0.992550552 333.000000 11.9646358
2.83366466 0.177820283E-2 336.000000 6.12744808
2.83368802 0.183327938E-3 338.000000 6.12689734
2.83369708 0.679280609E-3 337.000000 6.12836981
2.87943292 0.146849517E-3 332.000000 2.51892495
2.87945318 0.295515812E-3 331.000000 2.51745558
2.89804840 0.539837289E-3 330.000000 1.42265606
2.93557525 0.532299697 329.000000 12.5491867
2.93559408 0.214365616 391.000000 12.5504198
2.93561029 0.175688639 390.000000 12.5498590
2.97320676 0.158064649E-3 388.000000 9.35168839
2.97321057 0.300355150E-4 408.000000 9.35530090
2.97321224 0.931323011E-4 389.000000 9.35327435
2.97682190 0.551215744E-4 405.000000 8.87609196

Notice that in the files you will find many exciton with negative energy due the full-BSE nature, finally consider line-width numbers with care due to the approximation described above.

== Excitonic Eliashberg Functions ==
If you run the BSE and save the excitonic wave-functions (uncomment the flag <code>WRbsWF</code> in the BSE input file), tt is possible to plot the excitonic Eliashberg Functions with the command <code> ypp_ph -e e</code>.
For an interpretation of these functions see the discussion in reference <ref>Andrea Marini
[https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.101.106405 Phys. Rev. Lett. '''101''', 106405 (2008)]</ref>

== Bethe-Salpeter at finite temperature (inversion solver) ==

It is possible to take into account complex quasi-particles in the BSE solution by means of the inversion procedure.
The idea is to construct an Lo (Eq. 5 of Ref.<ref name='dbgrid'></ref>) that includes the electron-phonon corrections including the imaginary part.
In order to generate the input you can type: <code>yambo_ph -X s -o b -k sex -y i -V qp</code>

em1s # [R][Xs] Statically Screened Interaction
optics # [R] Linear Response optical properties
bss # [R] BSE solver
bse # [R][BSE] Bethe Salpeter Equation.
dipoles # [R] Oscillator strenghts (or dipoles)
DIP_Threads=0 # [OPENMP/X] Number of threads for dipoles
X_Threads=0 # [OPENMP/X] Number of threads for response functions
K_Threads=0 # [OPENMP/BSK] Number of threads for response functions
Chimod= "HARTREE" # [X] IP/Hartree/ALDA/LRC/PF/BSfxc
BSEmod= "resonant" # [BSE] resonant/retarded/coupling
BSKmod= "SEX" # [BSE] IP/Hartree/HF/ALDA/SEX/BSfxc
BSSmod= "i" # [BSS] (h)aydock/(d)iagonalization/(s)lepc/(i)nversion/(t)ddft`
BSENGexx= 9257 RL # [BSK] Exchange components
BSENGBlk=-1 RL # [BSK] Screened interaction block size [if -1 uses all the G-vectors of W(q,G,Gp)]
#WehCpl # [BSK] eh interaction included also in coupling
KfnQPdb= "E W < T0/ndb.QP" # [EXTQP BSK BSS] Database action
KfnQP_INTERP_NN= 1 # [EXTQP BSK BSS] Interpolation neighbours (NN mode)
KfnQP_INTERP_shells= 20.00000 # [EXTQP BSK BSS] Interpolation shells (BOLTZ mode)
KfnQP_DbGd_INTERP_mode= "NN" # [EXTQP BSK BSS] Interpolation DbGd mode
% KfnQP_E
 0.500000 | 1.000000 | 1.000000 | # [EXTQP BSK BSS] E parameters (c/v) eV|adim|adim
%
BSEprop= "abs" # [BSS] abs/kerr/magn/dichr trace
% BSEQptR
1 | 1 | # [BSK] Transferred momenta range
%
% BSEBands
 2 | 7 | # [BSK] Bands range
%
% BEnRange
 0.00000 | 5.00000 | eV # [BSS] Energy range
%
% BDmRange
 0.010000 | 0.010000 | eV # [BSS] Damping range
%
BEnSteps= 1500 # [BSS] Energy steps
% BLongDir
1.000000 | 0.000000 | 0.000000 | # [BSS] [cc] Electric Field
%
#WRbsWF # [BSS] Write to disk excitonic the WFs
XfnQPdb= "none" # [EXTQP Xd] Database action
% BndsRnXs
1 | 12 | # [Xs] Polarization function bands
%
NGsBlkXs= 113 RL # [Xs] Response block size
% LongDrXs
1.000000 | 0.000000 | 0.000000 | # [Xs] [cc] Electric Field
%
XTermKind= "none" # [X] X terminator ("none","BG" Bruneval-Gonze)

then if you run the job you will get the spectrum at finite temperature. Notice that with the inversion solver you can choose if you want the BSE with or without the coupling by changing the <code>BSEmod</code> variable.
This approach can be combined with double-grid as explained here: [http://www.attaccalite.com/speed-up-dielectric-constant-calculations-using-the-double-grid-method-with-yambo/ Speed up dielectric constant calculations using the double-grid method with Yambo ]. Yambo automatically will interpolate quasi-particle corrections on the fine grid. Notice that BSE inversion can become unfeasible for large BSE matrices.

Nota bene: if you use the inversion solver you just get the optical spectra, but you do not have access to the excitonic eigenvalues and eigenvectors and by consequence you cannot plot the excitonic wave-function or print the excitonic Eliashberg functions.

== Phonon-assisted density of states ==
Even if exciton-phonon scattering is not included in the BSE at finite temperature, it is possible to plot the phonon-assisted density of states (DOS) for light emission, defined as:

[[File:Formula exc ph.png | 600px | center | Phonon-assisted JDOS ]]

where '''Eql''' is the energy of the '''l'''-exciton at '''q'''-momentum, '''ωqλ''' is the phonon energy, '''nB''' is the Bose function, exciton are weighted with a Boltzamn factor and '''Emin''' is the lowest exciton energy in the all BZ.
In order to calculate these DOS you need to solve BSE for all the q-points for the lowest excitons, you can use SLEPC library to speed up calculations, then you need a converged phonon calculations.
Then you interpolate phonon on a dense phonon grid using <code>matdyn.x</code> utility in QE, and Yambo will interpolate exicton on the same grid.
Here we present an example for hBN that is an indirect semiconductor.
Running <code> ypp_ph -e p</code> you will get:

BoseTemp= 50 K # Bosonic Temperature
excitons # [R] Excitonic properties
ph_ass_dos # [R] Phonon-assisted DOS
States= "1 - 4" # Index of the BS state(s)
PHfreqF= "./bn.freq_54" # PWscf format file containing the phonon frequencies
% DOSERange
5.000000 | 5.500000 | eV # Energy range
%
DOSESteps= 1000 # Energy steps
DOS_broad= 0.004 eV # Broadening of the DOS

where "bn.freq_54" is the file produced by matdyn.x with the input:

&input
asr='simple',
flfrc='bn.fc',
flfrq='bn.freq_54',
dos=.true.,
fldos='bn.dos',
ndos=2
nk1=54, nk2=54, nk3=18
/

and the BSE was calculated on a 18x18x6 grid. The final spectra will be:

[[File:Test dos pl bulk.png| 700px | center| Phonon-assisted JDOS]]

where the DOS has been shifted to match the experimental peaks. Notice that the peak intensities are complitely off, because of the lack of the exciton-phonon matrix elements.

==Links==
* Back to [[ICTP 2022#Tutorials]]

==References==

Optical properties at finite temperature

2022-04-07T11:05:05Z

Pierre: /* Result analysis */

[[File:High low.png|thumb | 200px|Yambo tutorial image]]

In this tutorial, we will show you how to calculate optical properties including thermal effects due to the electron-phonon coupling. 
This tutorial assumes that you have completed all the steps from the previous tutorial on [[Electron Phonon Coupling]]. 
The tutorial is dived into different steps first we will calculate absorption at independent particle approximation, then we will include excitonic effects,
and finally, we will show how to analyze the data. 

'''Notice that there are two ways to include finite temperature effects in the BSE''':
* The first is to use a full Bethe-Salpeter equation with all the coupling terms in such a way to treat complex quasi-particles that cannot be handled in standard (Hermitian) BSE, this is explained in section [[https://www.yambo-code.org/wiki/index.php?title=Optical_properties_at_finite_temperature#Bethe-Salpeter_at_finite_temperature_.28diagonalization_solver.29 Bethe-Salpeter at finite temperature (diagonalization solver)]]. This approach has the advantage that one can obtain excitonic eigenvalues and eigenfunctions at finite temperature and not only the optical spectrum.

* The second method is use of a standard Hermitinan (resonant) BSE plus the inversion solver <ref name='dbgrid'>D. Kammerlander et al. [https://arxiv.org/abs/1209.1509 Phys. Rev. B '''86''', 125203] </ref>. In this case, a Lo two-particle Green's function that contains the complex quasi-particles is used in the solver to get a finite temperature spectrum. This approach has the advantage that can be used with the double-grid but it does not give access to eigenvalues and eigenfunctions of the BSE. This approach is described in the section: [[https://www.yambo-code.org/wiki/index.php?title=Optical_properties_at_finite_temperature#Bethe-Salpeter_at_finite_temperature_.28inversion_solver.29 Bethe-Salpeter at finite temperature (inversion solver)]]

== Absorption at finite temperature ==

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 <code>yambo_ph -o c -V qp -F yambo.in_300K</code>:

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= 1500 # [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. Run your calculations with the command: <code>yambo_ph -F yambo.in_300K -J T300</code> and do the same for the 0K case. 
 You can plot the results with gnuplot using the command: <code> plot 'o-T300.eps_q1_ip' w l t '300 K', 'o-T0.eps_q1_ip' w l t '0 K'</code>. Hereafter the result without and with electron-phonon coupling at two different temperatures:
[[File:Si absorption finite t.png|800px|center |Absorption of bulk silicon at finite temperature]]

The temperature effect is clearly visible in the figure.

== Bethe-Salpeter at finite temperature (diagonalization solver) ==

In this section we will calculate the Bethe-Salpeter at finite temperature using the full-BSE.

You can generate the input using the command <code>yambo_ph -X s -o b -k sex -y d -V qp -F yambo.in_BSE</code>

em1s # [R][Xs] Statically Screened Interaction
optics # [R] Linear Response optical properties
bss # [R] BSE solver
bse # [R][BSE] Bethe Salpeter Equation.
dipoles # [R] Oscillator strenghts (or dipoles)
DIP_Threads=0 # [OPENMP/X] Number of threads for dipoles
X_Threads=0 # [OPENMP/X] Number of threads for response functions
K_Threads=0 # [OPENMP/BSK] Number of threads for response functions
Chimod= "HARTREE" # [X] IP/Hartree/ALDA/LRC/PF/BSfxc
BSEmod= "coupling" # [BSE] resonant/retarded/coupling
BSKmod= "SEX" # [BSE] IP/Hartree/HF/ALDA/SEX/BSfxc
BSSmod= "d" # [BSS] (h)aydock/(d)iagonalization/(s)lepc/(i)nversion/(t)ddft`
BSENGexx= 9257 RL # [BSK] Exchange components
BSENGBlk=-1 RL # [BSK] Screened interaction block size [if -1 uses all the G-vectors of W(q,G,Gp)]
#WehCpl # [BSK] eh interaction included also in coupling
KfnQPdb= "E W < T0/ndb.QP" # [EXTQP BSK BSS] Database action
% KfnQP_E
0.500000 | 1.000000 | 1.000000 | # [EXTQP BSK BSS] E parameters (c/v) eV|adim|adim
%
BSEprop= "abs" # [BSS] Can be any among abs/jdos/kerr/magn/dich/photolum/esrt
BSEdips= "none" # [BSS] Can be "trace/none" or "xy/xz/yz" to define off-diagonal rotation plane
% BSEQptR
1 | 1 | # [BSK] Transferred momenta range
%
% BSEBands
 2 | 7 | # [BSK] Bands range
%
% BEnRange
0.000000 | 5.000000 | eV # [BSS] Energy range
%
% BDmRange
0.010000 | 0.010000 | eV # [BSS] Damping range
%
BEnSteps= 1500 # [BSS] Energy steps
% BLongDir
1.000000 | 0.000000 | 0.000000 | # [BSS] [cc] Electric Field
%
#WRbsWF # [BSS] Write to disk excitonic the WFs
% BndsRnXs
1 | 12 | # [Xs] Polarization function bands
%
NGsBlkXs= 113 RL # [Xs] Response block size
% LongDrXs
1.000000 | 0.000000 | 0.000000 | # [Xs] [cc] Electric Field
%
XTermKind= "none" # [X] X terminator ("none","BG" Bruneval-Gonze)

In this input file we asked Yambo to include finite temperature quasi-particle in the BSE with the line <code>KfnQPdb= "E W < T0/ndb.QP" </code>, notice that we also introduce an additional rigid shift of 0.5 eV to mimic the GW correction with the line <code> 0.500000 | 1.000000 | 1.000000 | \# [EXTQP BSK BSS] E parameters (c/v) eV|adim|adim</code>. You can now run the calculation with the command <code> yambo_ph -F yambo.in_BSE -J T0_BSE</code> 
In principle you can calculate the GW correction following the [https://www.yambo-code.org/wiki/index.php?title=How_to_obtain_the_quasi-particle_band_structure_of_a_bulk_material:_h-BN this tutorial] and then merge the corresponding ndb.QP database with the one of the electron-phonon coupling using the command <code> ypp -qpdb m</code>.
Notice that we changed the BSE type to "coupling" because you need the full Bethe-Salpeter to deal with complex quasi-particles. 
Now we do the same calculation at T = 300 K. Copy the previous input file <code> cp yambo.in_BSE yambo.in_BSE_300</code>. Then modify where the code reads the QP database :
KfnQPdb= "E W < T300/ndb.QP" # [EXTQP BSK BSS] Database action
Run the command with <code>yambo_ph -F yambo.in_BSE_300 -J T300_BSE</code> 
Finally, you should be able to plot the absorption spectra from the outputs in gnuplot :
<code> plot 'o-T0_BSE.eps_q1_diago_bse' w l t '0 K', 'o-T300_BSE.eps_q1_diago_bse' w l t '300 K'</code>
[[File:Si bse optics.png|800px|center |Bethe-Salpeter at finite temperature for bulk silicom]]

'''Approximations:''' 
In this calculation we have made different approximations: 
1) we did not include the renormalization of excitons due to the change of the screening potential W with temperature. Including electron-phonon coupling
in the dielectric constant by changing the line <code> XfnQPdb= "none" </code> is unfortunately not enough, for a discussion see ref. <ref>L. Adamska and P. Umari, [https://journals.aps.org/prb/abstract/10.1103/PhysRevB.103.075201 Phys. Rev. B '''103''', 075201 (2021)] </ref> 
2) The electron-phonon coupling should enter in the BSE through the exciton-phonon matrix elements and not from the single-particle self-energy, as we have done in this tutorial. The approximation used in this tutorial is valid for not too strong bound excitons, and in general it generated a finite life-time also for the lowest exciton in direct materials that is not correct. For a discussion see refs.<ref>G. Antonius, S. G. Louie [https://arxiv.org/abs/1705.04245 Phys. Rev. B '''105''', 085111 (2022)]</ref><ref name='claudio'> see Supp. Mat. of [https://journals.aps.org/prb/abstract/10.1103/PhysRevB.99.081109 Phys. Rev. B '''99''', 081109(R) (2019)]</ref><ref name='fulvio'>F. Paleari, [https://orbilu.uni.lu/handle/10993/41058 Phd Thesis (2019)]</ref> and <ref>H. Chen, D. Sangalli, and M. Bernardi [https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.125.107401 Phys. Rev. Lett. '''125''', 107401 (2020)]</ref>. 
3) we did not include scattering between exciton and phonons, therefore our results do not contain phonon-assisted absorption peaks, you can include these terms by finite difference displacements see refs. <ref name='fulvio'></ref>.

== Result analysis ==
Results of BSE at finite temperature can be analized in the same way of standard BSE. You can plot exciton binging energy as a function of the temperature, exciton wave-function and so, see [https://www.yambo-code.org/wiki/index.php?title=How_to_analyse_excitons Exciton Plot tutorial].
For example if you sort exciton with the command <code> ypp_ph -e s -J T0_BSE</code> for the case T=0K, in addition to the exciton energy now you will find its line-width in the files .
<code>o.exc_qpt1_E_sorted</code> and <code>o.exc_qpt1_I_sorted</code>.
#
# Maximum Residual Value = .77136E+05
#
# E [ev] Strength Index W [meV]
#
...............
2.82974172 0.388171151E-1 335.000000 11.9707289
2.82976580 0.110338017 334.000000 11.9676380
2.82979798 0.992550552 333.000000 11.9646358
2.83366466 0.177820283E-2 336.000000 6.12744808
2.83368802 0.183327938E-3 338.000000 6.12689734
2.83369708 0.679280609E-3 337.000000 6.12836981
2.87943292 0.146849517E-3 332.000000 2.51892495
2.87945318 0.295515812E-3 331.000000 2.51745558
2.89804840 0.539837289E-3 330.000000 1.42265606
2.93557525 0.532299697 329.000000 12.5491867
2.93559408 0.214365616 391.000000 12.5504198
2.93561029 0.175688639 390.000000 12.5498590
2.97320676 0.158064649E-3 388.000000 9.35168839
2.97321057 0.300355150E-4 408.000000 9.35530090
2.97321224 0.931323011E-4 389.000000 9.35327435
2.97682190 0.551215744E-4 405.000000 8.87609196

Notice that in the files you will find many exciton with negative energy due the full-BSE nature, finally consider line-width numbers with care due to the approximation described above.

== Excitonic Eliashberg Functions ==
If you run the BSE and save the excitonic wave-functions (uncomment the flag <code>WRbsWF</code> in the BSE input file), tt is possible to plot the excitonic Eliashberg Functions with the command <code> ypp_ph -e e</code>.
For an interpretation of these functions see the discussion in reference <ref>Andrea Marini
[https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.101.106405 Phys. Rev. Lett. '''101''', 106405 (2008)]</ref>

== Bethe-Salpeter at finite temperature (inversion solver) ==

It is possible to take into account complex quasi-particles in the BSE solution by means of the inversion procedure.
The idea is to construct an Lo (Eq. 5 of Ref.<ref name='dbgrid'></ref>) that includes the electron-phonon corrections including the imaginary part.
In order to generate the input you can type: <code>yambo_ph -X s -o b -k sex -y i -V qp</code>

em1s # [R][Xs] Statically Screened Interaction
optics # [R] Linear Response optical properties
bss # [R] BSE solver
bse # [R][BSE] Bethe Salpeter Equation.
dipoles # [R] Oscillator strenghts (or dipoles)
DIP_Threads=0 # [OPENMP/X] Number of threads for dipoles
X_Threads=0 # [OPENMP/X] Number of threads for response functions
K_Threads=0 # [OPENMP/BSK] Number of threads for response functions
Chimod= "HARTREE" # [X] IP/Hartree/ALDA/LRC/PF/BSfxc
BSEmod= "resonant" # [BSE] resonant/retarded/coupling
BSKmod= "SEX" # [BSE] IP/Hartree/HF/ALDA/SEX/BSfxc
BSSmod= "i" # [BSS] (h)aydock/(d)iagonalization/(s)lepc/(i)nversion/(t)ddft`
BSENGexx= 9257 RL # [BSK] Exchange components
BSENGBlk=-1 RL # [BSK] Screened interaction block size [if -1 uses all the G-vectors of W(q,G,Gp)]
#WehCpl # [BSK] eh interaction included also in coupling
KfnQPdb= "E W < T0/ndb.QP" # [EXTQP BSK BSS] Database action
KfnQP_INTERP_NN= 1 # [EXTQP BSK BSS] Interpolation neighbours (NN mode)
KfnQP_INTERP_shells= 20.00000 # [EXTQP BSK BSS] Interpolation shells (BOLTZ mode)
KfnQP_DbGd_INTERP_mode= "NN" # [EXTQP BSK BSS] Interpolation DbGd mode
% KfnQP_E
 0.500000 | 1.000000 | 1.000000 | # [EXTQP BSK BSS] E parameters (c/v) eV|adim|adim
%
BSEprop= "abs" # [BSS] abs/kerr/magn/dichr trace
% BSEQptR
1 | 1 | # [BSK] Transferred momenta range
%
% BSEBands
 2 | 7 | # [BSK] Bands range
%
% BEnRange
 0.00000 | 5.00000 | eV # [BSS] Energy range
%
% BDmRange
 0.010000 | 0.010000 | eV # [BSS] Damping range
%
BEnSteps= 1500 # [BSS] Energy steps
% BLongDir
1.000000 | 0.000000 | 0.000000 | # [BSS] [cc] Electric Field
%
#WRbsWF # [BSS] Write to disk excitonic the WFs
XfnQPdb= "none" # [EXTQP Xd] Database action
% BndsRnXs
1 | 12 | # [Xs] Polarization function bands
%
NGsBlkXs= 113 RL # [Xs] Response block size
% LongDrXs
1.000000 | 0.000000 | 0.000000 | # [Xs] [cc] Electric Field
%
XTermKind= "none" # [X] X terminator ("none","BG" Bruneval-Gonze)

then if you run the job you will get the spectrum at finite temperature. Notice that with the inversion solver you can choose if you want the BSE with or without the coupling by changing the <code>BSEmod</code> variable.
This approach can be combined with double-grid as explained here: [http://www.attaccalite.com/speed-up-dielectric-constant-calculations-using-the-double-grid-method-with-yambo/ Speed up dielectric constant calculations using the double-grid method with Yambo ]. Yambo automatically will interpolate quasi-particle corrections on the fine grid. Notice that BSE inversion can become unfeasible for large BSE matrices.

Nota bene: if you use the inversion solver you just get the optical spectra, but you do not have access to the excitonic eigenvalues and eigenvectors and by consequence you cannot plot the excitonic wave-function or print the excitonic Eliashberg functions.

== Phonon-assisted density of states ==
Even if exciton-phonon scattering is not included in the BSE at finite temperature, it is possible to plot the phonon-assisted density of states (DOS) for light emission, defined as:

[[File:Formula exc ph.png | 600px | center | Phonon-assisted JDOS ]]

where '''Eql''' is the energy of the '''l'''-exciton at '''q'''-momentum, '''ωqλ''' is the phonon energy, '''nB''' is the Bose function, exciton are weighted with a Boltzamn factor and '''Emin''' is the lowest exciton energy in the all BZ.
In order to calculate these DOS you need to solve BSE for all the q-points for the lowest excitons, you can use SLEPC library to speed up calculations, then you need a converged phonon calculations.
Then you interpolate phonon on a dense phonon grid using <code>matdyn.x</code> utility in QE, and Yambo will interpolate exicton on the same grid.
Here we present an example for hBN that is an indirect semiconductor.
Running <code> ypp_ph -e p</code> you will get:

BoseTemp= 50 K # Bosonic Temperature
excitons # [R] Excitonic properties
ph_ass_dos # [R] Phonon-assisted DOS
States= "1 - 4" # Index of the BS state(s)
PHfreqF= "./bn.freq_54" # PWscf format file containing the phonon frequencies
% DOSERange
5.000000 | 5.500000 | eV # Energy range
%
DOSESteps= 1000 # Energy steps
DOS_broad= 0.004 eV # Broadening of the DOS

where "bn.freq_54" is the file produced by matdyn.x with the input:

&input
asr='simple',
flfrc='bn.fc',
flfrq='bn.freq_54',
dos=.true.,
fldos='bn.dos',
ndos=2
nk1=54, nk2=54, nk3=18
/

and the BSE was calculated on a 18x18x6 grid. The final spectra will be:

[[File:Test dos pl bulk.png| 700px | center| Phonon-assisted JDOS]]

where the DOS has been shifted to match the experimental peaks. Notice that the peak intensities are complitely off, because of the lack of the exciton-phonon matrix elements.

==Links==
* Back to [[ICTP 2022#Tutorials]]

==References==

Optical properties at finite temperature

2022-04-07T11:03:04Z

Pierre: /* Bethe-Salpeter at finite temperature (diagonalization solver) */

[[File:High low.png|thumb | 200px|Yambo tutorial image]]

In this tutorial, we will show you how to calculate optical properties including thermal effects due to the electron-phonon coupling. 
This tutorial assumes that you have completed all the steps from the previous tutorial on [[Electron Phonon Coupling]]. 
The tutorial is dived into different steps first we will calculate absorption at independent particle approximation, then we will include excitonic effects,
and finally, we will show how to analyze the data. 

'''Notice that there are two ways to include finite temperature effects in the BSE''':
* The first is to use a full Bethe-Salpeter equation with all the coupling terms in such a way to treat complex quasi-particles that cannot be handled in standard (Hermitian) BSE, this is explained in section [[https://www.yambo-code.org/wiki/index.php?title=Optical_properties_at_finite_temperature#Bethe-Salpeter_at_finite_temperature_.28diagonalization_solver.29 Bethe-Salpeter at finite temperature (diagonalization solver)]]. This approach has the advantage that one can obtain excitonic eigenvalues and eigenfunctions at finite temperature and not only the optical spectrum.

* The second method is use of a standard Hermitinan (resonant) BSE plus the inversion solver <ref name='dbgrid'>D. Kammerlander et al. [https://arxiv.org/abs/1209.1509 Phys. Rev. B '''86''', 125203] </ref>. In this case, a Lo two-particle Green's function that contains the complex quasi-particles is used in the solver to get a finite temperature spectrum. This approach has the advantage that can be used with the double-grid but it does not give access to eigenvalues and eigenfunctions of the BSE. This approach is described in the section: [[https://www.yambo-code.org/wiki/index.php?title=Optical_properties_at_finite_temperature#Bethe-Salpeter_at_finite_temperature_.28inversion_solver.29 Bethe-Salpeter at finite temperature (inversion solver)]]

== Absorption at finite temperature ==

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 <code>yambo_ph -o c -V qp -F yambo.in_300K</code>:

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= 1500 # [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. Run your calculations with the command: <code>yambo_ph -F yambo.in_300K -J T300</code> and do the same for the 0K case. 
 You can plot the results with gnuplot using the command: <code> plot 'o-T300.eps_q1_ip' w l t '300 K', 'o-T0.eps_q1_ip' w l t '0 K'</code>. Hereafter the result without and with electron-phonon coupling at two different temperatures:
[[File:Si absorption finite t.png|800px|center |Absorption of bulk silicon at finite temperature]]

The temperature effect is clearly visible in the figure.

== Bethe-Salpeter at finite temperature (diagonalization solver) ==

In this section we will calculate the Bethe-Salpeter at finite temperature using the full-BSE.

You can generate the input using the command <code>yambo_ph -X s -o b -k sex -y d -V qp -F yambo.in_BSE</code>

em1s # [R][Xs] Statically Screened Interaction
optics # [R] Linear Response optical properties
bss # [R] BSE solver
bse # [R][BSE] Bethe Salpeter Equation.
dipoles # [R] Oscillator strenghts (or dipoles)
DIP_Threads=0 # [OPENMP/X] Number of threads for dipoles
X_Threads=0 # [OPENMP/X] Number of threads for response functions
K_Threads=0 # [OPENMP/BSK] Number of threads for response functions
Chimod= "HARTREE" # [X] IP/Hartree/ALDA/LRC/PF/BSfxc
BSEmod= "coupling" # [BSE] resonant/retarded/coupling
BSKmod= "SEX" # [BSE] IP/Hartree/HF/ALDA/SEX/BSfxc
BSSmod= "d" # [BSS] (h)aydock/(d)iagonalization/(s)lepc/(i)nversion/(t)ddft`
BSENGexx= 9257 RL # [BSK] Exchange components
BSENGBlk=-1 RL # [BSK] Screened interaction block size [if -1 uses all the G-vectors of W(q,G,Gp)]
#WehCpl # [BSK] eh interaction included also in coupling
KfnQPdb= "E W < T0/ndb.QP" # [EXTQP BSK BSS] Database action
% KfnQP_E
0.500000 | 1.000000 | 1.000000 | # [EXTQP BSK BSS] E parameters (c/v) eV|adim|adim
%
BSEprop= "abs" # [BSS] Can be any among abs/jdos/kerr/magn/dich/photolum/esrt
BSEdips= "none" # [BSS] Can be "trace/none" or "xy/xz/yz" to define off-diagonal rotation plane
% BSEQptR
1 | 1 | # [BSK] Transferred momenta range
%
% BSEBands
 2 | 7 | # [BSK] Bands range
%
% BEnRange
0.000000 | 5.000000 | eV # [BSS] Energy range
%
% BDmRange
0.010000 | 0.010000 | eV # [BSS] Damping range
%
BEnSteps= 1500 # [BSS] Energy steps
% BLongDir
1.000000 | 0.000000 | 0.000000 | # [BSS] [cc] Electric Field
%
#WRbsWF # [BSS] Write to disk excitonic the WFs
% BndsRnXs
1 | 12 | # [Xs] Polarization function bands
%
NGsBlkXs= 113 RL # [Xs] Response block size
% LongDrXs
1.000000 | 0.000000 | 0.000000 | # [Xs] [cc] Electric Field
%
XTermKind= "none" # [X] X terminator ("none","BG" Bruneval-Gonze)

In this input file we asked Yambo to include finite temperature quasi-particle in the BSE with the line <code>KfnQPdb= "E W < T0/ndb.QP" </code>, notice that we also introduce an additional rigid shift of 0.5 eV to mimic the GW correction with the line <code> 0.500000 | 1.000000 | 1.000000 | \# [EXTQP BSK BSS] E parameters (c/v) eV|adim|adim</code>. You can now run the calculation with the command <code> yambo_ph -F yambo.in_BSE -J T0_BSE</code> 
In principle you can calculate the GW correction following the [https://www.yambo-code.org/wiki/index.php?title=How_to_obtain_the_quasi-particle_band_structure_of_a_bulk_material:_h-BN this tutorial] and then merge the corresponding ndb.QP database with the one of the electron-phonon coupling using the command <code> ypp -qpdb m</code>.
Notice that we changed the BSE type to "coupling" because you need the full Bethe-Salpeter to deal with complex quasi-particles. 
Now we do the same calculation at T = 300 K. Copy the previous input file <code> cp yambo.in_BSE yambo.in_BSE_300</code>. Then modify where the code reads the QP database :
KfnQPdb= "E W < T300/ndb.QP" # [EXTQP BSK BSS] Database action
Run the command with <code>yambo_ph -F yambo.in_BSE_300 -J T300_BSE</code> 
Finally, you should be able to plot the absorption spectra from the outputs in gnuplot :
<code> plot 'o-T0_BSE.eps_q1_diago_bse' w l t '0 K', 'o-T300_BSE.eps_q1_diago_bse' w l t '300 K'</code>
[[File:Si bse optics.png|800px|center |Bethe-Salpeter at finite temperature for bulk silicom]]

'''Approximations:''' 
In this calculation we have made different approximations: 
1) we did not include the renormalization of excitons due to the change of the screening potential W with temperature. Including electron-phonon coupling
in the dielectric constant by changing the line <code> XfnQPdb= "none" </code> is unfortunately not enough, for a discussion see ref. <ref>L. Adamska and P. Umari, [https://journals.aps.org/prb/abstract/10.1103/PhysRevB.103.075201 Phys. Rev. B '''103''', 075201 (2021)] </ref> 
2) The electron-phonon coupling should enter in the BSE through the exciton-phonon matrix elements and not from the single-particle self-energy, as we have done in this tutorial. The approximation used in this tutorial is valid for not too strong bound excitons, and in general it generated a finite life-time also for the lowest exciton in direct materials that is not correct. For a discussion see refs.<ref>G. Antonius, S. G. Louie [https://arxiv.org/abs/1705.04245 Phys. Rev. B '''105''', 085111 (2022)]</ref><ref name='claudio'> see Supp. Mat. of [https://journals.aps.org/prb/abstract/10.1103/PhysRevB.99.081109 Phys. Rev. B '''99''', 081109(R) (2019)]</ref><ref name='fulvio'>F. Paleari, [https://orbilu.uni.lu/handle/10993/41058 Phd Thesis (2019)]</ref> and <ref>H. Chen, D. Sangalli, and M. Bernardi [https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.125.107401 Phys. Rev. Lett. '''125''', 107401 (2020)]</ref>. 
3) we did not include scattering between exciton and phonons, therefore our results do not contain phonon-assisted absorption peaks, you can include these terms by finite difference displacements see refs. <ref name='fulvio'></ref>.

== Result analysis ==
Results of BSE at finite temperature can be analized in the same way of standard BSE. You can plot exciton binging energy as a function of the temperature, exciton wave-function and so, see [https://www.yambo-code.org/wiki/index.php?title=How_to_analyse_excitons Exciton Plot tutorial].
For example if you sort exciton with the command <code> ypp_ph -e s</code> for the case T=0K, in addition to the exciton energy now you will find its line-width in the files .
<code>o.exc_qpt1_E_sorted</code> and <code>o.exc_qpt1_I_sorted</code>.
#
# Maximum Residual Value = .77136E+05
#
# E [ev] Strength Index W [meV]
#
...............
2.82974172 0.388171151E-1 335.000000 11.9707289
2.82976580 0.110338017 334.000000 11.9676380
2.82979798 0.992550552 333.000000 11.9646358
2.83366466 0.177820283E-2 336.000000 6.12744808
2.83368802 0.183327938E-3 338.000000 6.12689734
2.83369708 0.679280609E-3 337.000000 6.12836981
2.87943292 0.146849517E-3 332.000000 2.51892495
2.87945318 0.295515812E-3 331.000000 2.51745558
2.89804840 0.539837289E-3 330.000000 1.42265606
2.93557525 0.532299697 329.000000 12.5491867
2.93559408 0.214365616 391.000000 12.5504198
2.93561029 0.175688639 390.000000 12.5498590
2.97320676 0.158064649E-3 388.000000 9.35168839
2.97321057 0.300355150E-4 408.000000 9.35530090
2.97321224 0.931323011E-4 389.000000 9.35327435
2.97682190 0.551215744E-4 405.000000 8.87609196

Notice that in the files you will find many exciton with negative energy due the full-BSE nature, finally consider line-width numbers with care due to the approximation described above.

== Excitonic Eliashberg Functions ==
If you run the BSE and save the excitonic wave-functions (uncomment the flag <code>WRbsWF</code> in the BSE input file), tt is possible to plot the excitonic Eliashberg Functions with the command <code> ypp_ph -e e</code>.
For an interpretation of these functions see the discussion in reference <ref>Andrea Marini
[https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.101.106405 Phys. Rev. Lett. '''101''', 106405 (2008)]</ref>

== Bethe-Salpeter at finite temperature (inversion solver) ==

It is possible to take into account complex quasi-particles in the BSE solution by means of the inversion procedure.
The idea is to construct an Lo (Eq. 5 of Ref.<ref name='dbgrid'></ref>) that includes the electron-phonon corrections including the imaginary part.
In order to generate the input you can type: <code>yambo_ph -X s -o b -k sex -y i -V qp</code>

em1s # [R][Xs] Statically Screened Interaction
optics # [R] Linear Response optical properties
bss # [R] BSE solver
bse # [R][BSE] Bethe Salpeter Equation.
dipoles # [R] Oscillator strenghts (or dipoles)
DIP_Threads=0 # [OPENMP/X] Number of threads for dipoles
X_Threads=0 # [OPENMP/X] Number of threads for response functions
K_Threads=0 # [OPENMP/BSK] Number of threads for response functions
Chimod= "HARTREE" # [X] IP/Hartree/ALDA/LRC/PF/BSfxc
BSEmod= "resonant" # [BSE] resonant/retarded/coupling
BSKmod= "SEX" # [BSE] IP/Hartree/HF/ALDA/SEX/BSfxc
BSSmod= "i" # [BSS] (h)aydock/(d)iagonalization/(s)lepc/(i)nversion/(t)ddft`
BSENGexx= 9257 RL # [BSK] Exchange components
BSENGBlk=-1 RL # [BSK] Screened interaction block size [if -1 uses all the G-vectors of W(q,G,Gp)]
#WehCpl # [BSK] eh interaction included also in coupling
KfnQPdb= "E W < T0/ndb.QP" # [EXTQP BSK BSS] Database action
KfnQP_INTERP_NN= 1 # [EXTQP BSK BSS] Interpolation neighbours (NN mode)
KfnQP_INTERP_shells= 20.00000 # [EXTQP BSK BSS] Interpolation shells (BOLTZ mode)
KfnQP_DbGd_INTERP_mode= "NN" # [EXTQP BSK BSS] Interpolation DbGd mode
% KfnQP_E
 0.500000 | 1.000000 | 1.000000 | # [EXTQP BSK BSS] E parameters (c/v) eV|adim|adim
%
BSEprop= "abs" # [BSS] abs/kerr/magn/dichr trace
% BSEQptR
1 | 1 | # [BSK] Transferred momenta range
%
% BSEBands
 2 | 7 | # [BSK] Bands range
%
% BEnRange
 0.00000 | 5.00000 | eV # [BSS] Energy range
%
% BDmRange
 0.010000 | 0.010000 | eV # [BSS] Damping range
%
BEnSteps= 1500 # [BSS] Energy steps
% BLongDir
1.000000 | 0.000000 | 0.000000 | # [BSS] [cc] Electric Field
%
#WRbsWF # [BSS] Write to disk excitonic the WFs
XfnQPdb= "none" # [EXTQP Xd] Database action
% BndsRnXs
1 | 12 | # [Xs] Polarization function bands
%
NGsBlkXs= 113 RL # [Xs] Response block size
% LongDrXs
1.000000 | 0.000000 | 0.000000 | # [Xs] [cc] Electric Field
%
XTermKind= "none" # [X] X terminator ("none","BG" Bruneval-Gonze)

then if you run the job you will get the spectrum at finite temperature. Notice that with the inversion solver you can choose if you want the BSE with or without the coupling by changing the <code>BSEmod</code> variable.
This approach can be combined with double-grid as explained here: [http://www.attaccalite.com/speed-up-dielectric-constant-calculations-using-the-double-grid-method-with-yambo/ Speed up dielectric constant calculations using the double-grid method with Yambo ]. Yambo automatically will interpolate quasi-particle corrections on the fine grid. Notice that BSE inversion can become unfeasible for large BSE matrices.

Nota bene: if you use the inversion solver you just get the optical spectra, but you do not have access to the excitonic eigenvalues and eigenvectors and by consequence you cannot plot the excitonic wave-function or print the excitonic Eliashberg functions.

== Phonon-assisted density of states ==
Even if exciton-phonon scattering is not included in the BSE at finite temperature, it is possible to plot the phonon-assisted density of states (DOS) for light emission, defined as:

[[File:Formula exc ph.png | 600px | center | Phonon-assisted JDOS ]]

where '''Eql''' is the energy of the '''l'''-exciton at '''q'''-momentum, '''ωqλ''' is the phonon energy, '''nB''' is the Bose function, exciton are weighted with a Boltzamn factor and '''Emin''' is the lowest exciton energy in the all BZ.
In order to calculate these DOS you need to solve BSE for all the q-points for the lowest excitons, you can use SLEPC library to speed up calculations, then you need a converged phonon calculations.
Then you interpolate phonon on a dense phonon grid using <code>matdyn.x</code> utility in QE, and Yambo will interpolate exicton on the same grid.
Here we present an example for hBN that is an indirect semiconductor.
Running <code> ypp_ph -e p</code> you will get:

BoseTemp= 50 K # Bosonic Temperature
excitons # [R] Excitonic properties
ph_ass_dos # [R] Phonon-assisted DOS
States= "1 - 4" # Index of the BS state(s)
PHfreqF= "./bn.freq_54" # PWscf format file containing the phonon frequencies
% DOSERange
5.000000 | 5.500000 | eV # Energy range
%
DOSESteps= 1000 # Energy steps
DOS_broad= 0.004 eV # Broadening of the DOS

where "bn.freq_54" is the file produced by matdyn.x with the input:

&input
asr='simple',
flfrc='bn.fc',
flfrq='bn.freq_54',
dos=.true.,
fldos='bn.dos',
ndos=2
nk1=54, nk2=54, nk3=18
/

and the BSE was calculated on a 18x18x6 grid. The final spectra will be:

[[File:Test dos pl bulk.png| 700px | center| Phonon-assisted JDOS]]

where the DOS has been shifted to match the experimental peaks. Notice that the peak intensities are complitely off, because of the lack of the exciton-phonon matrix elements.

==Links==
* Back to [[ICTP 2022#Tutorials]]

==References==