Nonequilibrium absorption in bulk silicon: Difference between revisions

From The Yambo Project
Jump to navigation Jump to search
No edit summary
 
(34 intermediate revisions by 2 users not shown)
Line 5: Line 5:
= The material: Silicon =
= The material: Silicon =


We will study nonequilibrium absorption in bulk silicon. The same material used for this [Silicon|GW tutorial]  
We will study nonequilibrium absorption in bulk silicon. The same material used for this [[Silicon|GW tutorial]]  


* [http://cst-www.nrl.navy.mil/lattice/struk/a1.html FCC] lattice
* [http://cst-www.nrl.navy.mil/lattice/struk/a1.html FCC] lattice
Line 45: Line 45:


= Equilibrium optical properties =
= Equilibrium optical properties =
We start from a density-functional theory (DFT) calculation. Then we apply MBPT to obtain the QP energies within the GW approximation and the absorption spectra from the solution of the BSE. The QP band structure will show the effect of the screened electron-electron interaction, while the absorption spectrum within the BSE accounts for the effect of the electron-hole interaction. For more information about solving the Equilibrium BSE equation refer to the following tutorial: [https://wiki.yambo-code.eu/wiki/index.php/Calculating%20optical%20spectra%20including%20excitonic%20effects:%20a%20step-by-step%20guide Calculating optical spectra including excitonic effects: a step-by-step guide]
We start from a density-functional theory (DFT) calculation. Then we apply MBPT to obtain the QP energies within the GW approximation and the absorption spectra from the solution of the BSE. The QP band structure will show the effect of the screened electron-electron interaction, while the absorption spectrum within the BSE accounts for the effect of the electron-hole interaction. When QP corrections are not available, they can be included using the scissor operator <code>KfnQP_E</code>, and the numerical values for "corrections" can be deduced from a GW calculation. For more information about solving the Equilibrium BSE equation refer to the following tutorial: [https://wiki.yambo-code.eu/wiki/index.php/Calculating%20optical%20spectra%20including%20excitonic%20effects:%20a%20step-by-step%20guide Calculating optical spectra including excitonic effects: a step-by-step guide]


(to delete open) In this case, do we actually compute the qp eng withing gw approximation? or we directy solve the bse eq (to delete close)  
(to delete open) In this case, do we actually compute the qp eng withing gw approximation? or we directy solve the bse eq (to delete close)  
Line 58: Line 58:
First, to solve the BSE equation we need to calculate the statically screened electron-electron interaction (W), which is related to static diaelectric screening <math> \epsilon^{-1} </math>.  
First, to solve the BSE equation we need to calculate the statically screened electron-electron interaction (W), which is related to static diaelectric screening <math> \epsilon^{-1} </math>.  


<math display="block">\epsilon^{-1}_{\mathbf{G},\mathbf{G'}} (\mathbf{q}, \omega = 0) = \delta_{\mathbf{G},\mathbf{G'}} + v(\mathbf{v} + \mathbf{G})_{\chi,\mathbf{G},\mathbf{G'}}(\mathbf{q}, \omega = 0)</math>
<math> \Large
\epsilon^{-1}_{\mathbf{G},\mathbf{G'}} (\mathbf{q}, \omega = 0) = \delta_{\mathbf{G},\mathbf{G'}} + v(\mathbf{v} + \mathbf{G}) \chi_{\mathbf{G},\mathbf{G'}(\mathbf{q}}, \omega = 0)</math>


To do this create the following input file:  
To do this create the following input file:  
  $ yambo -X s -F 01_EQ_BSE_screening.in
  $ yambo -X s -F 01_EQ_BSE_screening.in


Line 83: Line 82:
</div>
</div>


Where we changed the response block size <code>NGsBlkXs = 3</code> and number of polarization function bands <code> BndsRnXs </code> to be in the range of 1 to 30.  
Where we changed the response block size <code>NGsBlkXs = 3 Ry</code> and number of polarization function bands <code>BndsRnXs</code> to be in the range of 1 to 30. Now run the file:
 
$ yambo -F 01_EQ_BSE_screening.in -J EQ_BSE


== Solving the Bethe-Salpeter equation ==  
== Solving the Bethe-Salpeter equation ==  
To solve the Bethe-Salpeter equation, the latter is usually rewritten in the space of transitions between valence and conduction states as the (pseudo)eigenvalue problem for a two-particle Hamiltonian. For a non spin-polarized system and in the <math> q = 0 </math> limit, the two-particle Hamiltonian matrix elements are given by
<math>
\Large
H_{vc\mathbf{k},v'c'\mathbf{k'}} = (\varepsilon_{c\mathbf{k}} - \varepsilon_{v\mathbf{k}})\delta_{vv'}\delta_{cc'}\delta_{\mathbf{k},\mathbf{k'}} + (f_{c\mathbf{k}} - f_{v\mathbf{k}})\left[ 2V_{vc\mathbf{k},v'c'\mathbf{k'}} - W_{vc\mathbf{k},v'c'\mathbf{k'}}\right]
</math>
where <math>vc\mathbf{k}</math> indicates the pair of quasiparticle states <math>v\mathbf{k}</math> and <math>c\mathbf{k}</math>. The first term on the RHS is the quasiparticle energy differences (diagonal only). The second term is the kernel which is the sum of the electron-hole exchange part V (which stems from the Hartree potential) and the electron-hole attraction part W (which stems from the screened exchange potential). The electron-hole attraction part W, depends on <math> \epsilon^{-1} </math> which was calculated in the previous section.
After calculating the equilibrium Bethe-Salpeter kernel, the resulting two-particle Hamiltonian can be diagonalized to compute the macroscopic dielectric function, which can be used to get absorption and EELS spectra. Calculation of BS kernel and diagonalization of Hamiltonian can be done in one step.
To generate the input file for the calculation of BSE kernel, run the following command:


$ yambo -o b -k sex -F 02_EQ_BSE_kernel_diag.in -y d -J EQ_BSE
with the following input
<div class="mw-collapsible mw-collapsed" style="width:1100px; overflow:auto;">
<div style="font-weight:bold;line-height:1.6;"> 02_EQ_BSE_kernel_diag.in file (Expandable) </div>
<div class="mw-collapsible-content">
<pre>
optics                          # [R] Linear Response optical properties
bss                              # [R] BSE solver
bse                              # [R][BSE] Bethe Salpeter Equation.
K_Threads=0                      # [OPENMP/BSK] Number of threads for response functions
BSKmod= "SEX"                    # [BSE] IP/Hartree/HF/ALDA/SEX/BSfxc
BSEmod= "resonant"              # [BSE] resonant/retarded/coupling
BSSmod= "d"                      # [BSS] (h)aydock/(d)iagonalization/(s)lepc/(i)nversion/(t)ddft`
BSENGexx= 2085            RL    # [BSK] Exchange components
BSENGBlk=-1                RL    # [BSK] Screened interaction block size [if -1 uses all the G-vectors of W(q,G,Gp)]
% BSEQptR
1 | 1 |                            # [BSK] Transferred momenta range
%
% BSEBands
  2 |  7 |                        # [BSK] Bands range
%
% BEnRange
  2.20000 | 5.00000 |        eV    # [BSS] Energy range
%
% KfnQP_E
0.70000 | 1.000000 | 1.000000 |      # [EXTQP BSK BSS] E parameters (c/v)
%
% BDmRange
0.100000 | 0.100000 |        eV    # [BSS] Damping range
%
BEnSteps= 200                    # [BSS] Energy steps
% BLongDir
1.000000 | 0.000000 | 0.000000 |        # [BSS] [cc] Electric Field versor
%
BSEprop= "abs"                  # [BSS] Can be any among abs/jdos/kerr/asymm/anHAll/magn/dich/photolum/esrt
BSEdips= "none"                  # [BSS] Can be "trace/none" or "xy/xz/yz" to define off-diagonal rotation plane
</pre>
</div>
</div>
Where <code>BSEBands</code> was changed to include bands from 3 to 8, which defines the basis of quasiparticle pairs that are used to describe excitons. 2 valence and 4 conduction bands where included. <code>BEnRange</code> was increased to from 0-10 eV to 0-16 eV, and <code>BEnSteps</code> was increased from 100 to 200. The scissors operator <code>KfnQP_E</code> was used to match the positon calculated spectra with experimental values.
To run the file:
$ yambo -F 02_EQ_BSE_kernel_diag.in -J EQ_BSE
After the calculation has finished, the absorption spectra can be plotted with the following:
<div class="mw-collapsible mw-collapsed" style="width:1100px; overflow:auto;">
<div style="font-weight:bold;line-height:1.6;"> gnuplot commands (Expandable) </div>
<div class="mw-collapsible-content">
<pre>
$ gnuplot
...
gnuplot> set x label "Energy (eV)"
gnuplot> set y label "Absorption"
gnuplot> plot 'o-EQ_BSE.eps_q1_diago_bse' u 1:2 w l t '-k sex'
</pre>
</div>
</div>
and the following plot will results (the "experimental" and "-k hartree" plots are added additionally for comparison)
<div class="mw-collapsible mw-collapsed" style="width:1200px; overflow:auto;">
<div style="font-weight:bold;line-height:1.6;"> Plot of the computed spectra (Expandable) </div>
<div class="mw-collapsible-content">
[[File:Plot of Si absorption spectra.png|1200px|frameless|Plots of computed absorption spectra for the equilibrium case.]]
</div>
</div>


= Generating non-equilibrium carriers =
= Generating non-equilibrium carriers =
To simulate non-equilibrium absorption, we need to define the number of non-equilibrium carriers, specifically electrons in the conduction band and holes in the valence band. These carriers follow the Fermi-Dirac distribution given by:
<math>\Large f(E) = 1 / (exp((E - μ) / k_B T) + 1)</math>
Here:
* E is the energy level,
* <math>\mu</math> is the chemical potential (specific to electrons or holes),
* <math>T</math> is the effective temperature,
* <math>k_B</math> is the Boltzmann constant.
For the simulation, both electrons and holes are characterized by their own <math>\mu_i</math> and <math>T_i</math> .
The distribution of non-equilibrium carriers can be generated using Yambo’s post processing for real-time module (<code>ypp_rt</code>), which allows modeling of these conditions.
To generate to non-equilibrium carriers databases create <code>input_ypp1 input_ypp2</code> files, with the following content:
<div class="mw-collapsible mw-collapsed" style="width:800px; overflow:auto;">
<div style="font-weight:bold;line-height:1.6;"> input_ypp1 (Expandable) </div>
<div class="mw-collapsible-content">
<pre>
RTDBs                      # [R] Real-Time databases
Select_Fermi              # [R] NEQ DBs according to Fermi distribution
% RTBands
  1 | 8 |              # [RT] Bands range
%
% RTmuEh
-0.80 | 0.80 | eV          # [RT] Chemical potentials hole | electron
%
% RTtempEh
0.04 | 0.04 | eV      # [RT] Effective temperature hole | electron
%
RTautotuneThr= 0.000100    # [RT] Threshold to match no. of pumped holes and electrons.
</pre>
</div>
</div>
<div class="mw-collapsible mw-collapsed" style="width:800px; overflow:auto;">
<div style="font-weight:bold;line-height:1.6;"> input_ypp2 (Expandable) </div>
<div class="mw-collapsible-content">
<pre>
RealTime                  # [R] Real-Time Post-Processing
RToccupations              # [R] Analize time-dependent occupations
RTenergy                  # [R] Post-Processing kind: function of energy
TimeStep= 0.000000 fs  # Time step
% TimeRange
-1.000000 |-1.000000 | fs  # Time-window where processing is done
%
#IncludeEQocc              # [NEGF] Include also equilibrium occupations
#SkipFermiFIT              # [NEGF] Do a Fermi Fit of occupations (and lifetimes ratio)
%QPkrange                #  generalized Kpoint/Band indices
  1| 3|  1|8|
%
</pre>
</div>
</div>


And then run the files:
$ ypp_rt
$ ypp_rt -F input_ypp1 -J carriers -C carriers
$ ypp_rt -F input_ypp2 -J carriers -C carriers
The file <code>input_ypp1</code> will generate a <code>ndb.RT_carriers</code> database which will contain the Fermi distribution of the electrons and the hole according to the specified chemical potential, <code>RTmuEh</code>, and for a specific temperature given by <code>RTtempEh</code>.  The file <code>input_ypp2</code> will generate a series of plots from the carriers database.


= Energy shift in the band structure =
= Energy shift in the band structure =
Line 95: Line 241:
The GW method is the standard approach to compute quasi-particle corrections. However, in presence of non-equilibrium carriers, or even at finite temperature the formulation of the GW self-energy needs to be refined. This is due to its frequency dependence, which results in extra terms when performing the analytic continuation from the Keldish contour (non-equilibrium case) or from the Matsubara axis (finite temperature case) to the real time/frequency axis. See for example this reference <ref name="faleev2006"/> for the finite temperature case.
The GW method is the standard approach to compute quasi-particle corrections. However, in presence of non-equilibrium carriers, or even at finite temperature the formulation of the GW self-energy needs to be refined. This is due to its frequency dependence, which results in extra terms when performing the analytic continuation from the Keldish contour (non-equilibrium case) or from the Matsubara axis (finite temperature case) to the real time/frequency axis. See for example this reference <ref name="faleev2006"/> for the finite temperature case.
On the other hand the COHSEX self-energy, being static, avoids this complication. This is why we will compute changes in the QP corrections within the COHSEX approximation.
On the other hand the COHSEX self-energy, being static, avoids this complication. This is why we will compute changes in the QP corrections within the COHSEX approximation.
To compute the Non-Equilibrium Absorption Spectra, we need to diagonalize a previously defined two-particle Hamiltonian. To account for the Non-Equilibrium effect, the Kohn-Sham energies, <math>\varepsilon_{c\mathbf{k}}, \varepsilon_{v\mathbf{k}}</math>, calculated at the DFT level, need to be corrected. This requires calculating the following: COHSEX self-energy at equilibrium, COHSEX self-energy at non-equilibrium, and Equilibrium GW self-energy at Plasmon-Pole Approximation level.
Using these, the energy correction, <math>\varepsilon^{corr,neq}</math>, can be computed as:
<math>\Large \varepsilon^{corr,neq}_{c\mathbf{k}} = \varepsilon^{GW,eq}_{c\mathbf{k}} + (\varepsilon^{cohsex,neq}_{c\mathbf{k}} - \varepsilon^{cohsex,eq}_{c\mathbf{k}})</math>
The KS energies will be corrected as follows:
<math>\Large \varepsilon^{KS,neq}_{c\mathbf{k}} = \varepsilon_{c\mathbf{k}} + \varepsilon^{corr,neq}_{c\mathbf{k}}</math>
These corrected KS energies will be used when diagonalizing the two-particle Hamiltonian.


== COHSEX corrections at equilibrium ==
== COHSEX corrections at equilibrium ==
To make a COHSEX calculation at equilibrium, we do the following:


== Screening in presence of non-equilibrium carriers ==
$ yambo -p c -F 03_eq_cohsex.in
 
which will create a file with the following content:
 
<div class="mw-collapsible mw-collapsed" style="width:1000px; overflow:auto;">
<div style="font-weight:bold;line-height:1.6;"> 03_eq_cohsex.in (Expandable) </div>
<div class="mw-collapsible-content">
<pre>
gw0                              # [R] GW approximation
cohsex                          # [R][Xp] COlumb Hole Screened EXchange
dyson                            # [R] Dyson Equation solver
el_el_corr                      # [R] Electron-Electron Correlation
HF_and_locXC                    # [R] Hartree-Fock
em1s                            # [R][Xs] Statically Screened Interaction
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
SE_Threads=0                    # [OPENMP/GW] Number of threads for self-energy
EXXRLvcs= 2085            RL    # [XX] Exchange    RL components
VXCRLvcs= 2085            RL    # [XC] XCpotential RL components
Chimod= "HARTREE"                # [X] IP/Hartree/ALDA/LRC/PF/BSfxc
% BndsRnXs
  1 |  10 |                        # [Xs] Polarization function bands
%
NGsBlkXs= 1                RL    # [Xs] Response block size
% LongDrXs
1.000000 | 0.000000 | 0.000000 |        # [Xs] [cc] Electric Field
%
%QPkrange                        # [GW] QP generalized Kpoint/Band indices
1|8|1|10|
%
</pre>
</div>
</div>
 
Where the number of Polarization function bands, <code>BndsRnXs</code> was changed from <code>1|50</code> to <code>1|10</code>, and the quasi-particles generalized Kpoint and Band indices,<code>QPkrange</code>, was changed from <code>1|29|1|50</code> to <code>1|8|1|10</code>. And after it we run the file with the following command:
 
$ yambo_rt -F 03_eq_cohsex.in -J cohsex_eq_emp -C cohsex_eq_emp


== COHSEX corrections in presence of non-equilibrium carriers ==
== COHSEX corrections in presence of non-equilibrium carriers ==
Before making a COHSEX calculation, we need to recalculate the static screening in the presence of non-equilibrium carriers. After it, we can make a COHSEX calculation in the presence of non-equilibrium carriers using the <code>ndb.RT_carriers</code> database generated in the previous step with the <code>ypp_rt</code> module. To do this, we do the following:
$ yambo -p c -F 03_neq_cohsex.in
where we need to add the following lines:
XfnRTdb= "f @ 0. ps < carriers/ndb.RT_carriers"
GfnRTdb= "f @ 0. ps < carriers/ndb.RT_carriers"


and the resulting file is:
<div class="mw-collapsible mw-collapsed" style="width:1000px; overflow:auto;">
<div style="font-weight:bold;line-height:1.6;"> 03_neq_cohsex.in (Expandable) </div>
<div class="mw-collapsible-content">
<pre>
gw0                              # [R] GW approximation
cohsex                          # [R][Xp] COlumb Hole Screened EXchange
dyson                            # [R] Dyson Equation solver
el_el_corr                      # [R] Electron-Electron Correlation
HF_and_locXC                    # [R] Hartree-Fock
em1s                            # [R][Xs] Statically Screened Interaction
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
SE_Threads=0                    # [OPENMP/GW] Number of threads for self-energy
EXXRLvcs= 2085            RL    # [XX] Exchange    RL components
VXCRLvcs= 2085            RL    # [XC] XCpotential RL components
Chimod= "HARTREE"                # [X] IP/Hartree/ALDA/LRC/PF/BSfxc
% BndsRnXs
  1 |  10 |                        # [Xs] Polarization function bands
%
NGsBlkXs= 1                RL    # [Xs] Response block size
% LongDrXs
1.000000 | 0.000000 | 0.000000 |        # [Xs] [cc] Electric Field
%
%QPkrange                        # [GW] QP generalized Kpoint/Band indices
1|8|1|10|
%
</pre>
</div>
</div>
where we added the respective two lines with <code>XfnRTdb</code>, and <code>GfnRTdb</code>, and made changes to <code>BndsRnXs</code> and <code>QPkrange</code> similar to the equilibrium case.
The line
> em1s                            # [R][Xs] Statically Screened Interaction
will tell Yambo to make the static screening calculation, and the lines:
> gw0                              # [R] GW approximation
> cohsex                          # [R][Xp] COlumb Hole Screened EXchange
will tell Yambo to do make a COHSEX calculation.
After it we can run the file with the command:
$ yambo_rt -F 03_neq_cohsex.in -J cohsex_neq -C cohsex_neq
== GW PPA corrections for the equilibrium case ==
This subchapter can be skipped, and replaced by using scissors operator when diagonalizing the two particle Hamiltonian.
To get accurate values for the band gap value we need to go beyond COHSEX approximation. Therefore we will obtain the Quasi-particle corrections in the plasmon pole approximation for equilibrium case, and in next steps, to this value we will add the correction calculated at COHSEX approximation to account for the non-equilibrium state of the system.
To make this calculation, we do the following:
$ yambo -p p -F 03_eq_qp_ppa.in
which will create a file with the following content:
<div class="mw-collapsible mw-collapsed" style="width:1000px; overflow:auto;">
<div style="font-weight:bold;line-height:1.6;"> 03_eq_qp_ppa.in (Expandable) </div>
<div class="mw-collapsible-content">
<pre>
gw0                              # [R] GW approximation
ppa                              # [R][Xp] Plasmon Pole Approximation for the Screened Interaction
el_el_corr                      # [R] Electron-Electron Correlation
HF_and_locXC                    # [R] Hartree-Fock
em1d                            # [R][X] Dynamically Screened Interaction
X_Threads=0                      # [OPENMP/X] Number of threads for response functions
DIP_Threads=0                    # [OPENMP/X] Number of threads for dipoles
SE_Threads=0                    # [OPENMP/GW] Number of threads for self-energy
EXXRLvcs= 2085            RL    # [XX] Exchange    RL components
VXCRLvcs= 2085            RL    # [XC] XCpotential RL components
Chimod= "HARTREE"                # [X] IP/Hartree/ALDA/LRC/PF/BSfxc
% BndsRnXp
  1 |  10 |                        # [Xp] Polarization function bands
%
NGsBlkXp= 1                RL    # [Xp] Response block size
% LongDrXp
1.000000 | 0.000000 | 0.000000 |        # [Xp] [cc] Electric Field
%
PPAPntXp= 27.21138        eV    # [Xp] PPA imaginary energy
% GbndRnge
  1 |  10 |                        # [GW] G[W] bands range
%
GTermKind= "none"                # [GW] GW terminator ("none","BG" Bruneval-Gonze,"BRS" Berger-Reining-Sottile)
DysSolver= "n"                  # [GW] Dyson Equation solver ("n","s","g","q")
%QPkrange                        # [GW] QP generalized Kpoint/Band indices
1|8|1|10|
%       
</pre>
</div>
</div>
where the following changes have been made: the number of polarization function bands, <code>BndsRnXp</code> was changed from <code>1|50</code> to <code>1|10</code>, the G[W] bands range, <code>GbndRnge</code> was changed from <code>1|50</code> to <code>1|10</code>, and quasi-particle generalized k points  and band indices, <code>QPkrange</code> was changed from <code>1|29|1|50</code> to <code>1|8|1|10</code>. And after it we run the file:
$ yambo_rt -F 03_eq_qp_ppa.in -J gw_eq -C gw_eq
== Data post-processing ==
Now after obtaining all the energy components that we need to calculate corrections to Kohn-Sham energies for the Non-Equilibrium case, we can manipulate them using the yambo's real time post processing module, <code>ypp_rt</code>. We want to obtain energies defined by the following equation:
<math>\Large \varepsilon^{corr,neq}_{c\mathbf{k}} = \varepsilon^{GW,eq}_{c\mathbf{k}} + (\varepsilon^{cohsex,neq}_{c\mathbf{k}} - \varepsilon^{cohsex,eq}_{c\mathbf{k}})</math>.
Therefore we create a yambo post processing file with the following contents:
<div class="mw-collapsible mw-collapsed" style="width:1000px; overflow:auto;">
<div style="font-weight:bold;line-height:1.6;"> input_ypp3 (Expandable) </div>
<div class="mw-collapsible-content">
<pre>
QPDBs                      # [R] Quasi-particle databases
QPDB_merge                # [R] Mergering
%Actions_and_names      # [QPDB] QP databases and Actions (format is "what"|"OP"|"prefactor"|"DB"|). OP can be ±/x.
"E" | "+" | "1" | "cohsex_neq/ndb.QP" |
"E" | "-" | "1" | "cohsex_eq_emp/ndb.QP" |
"Z" | "+" | "1" | "gw_eq/ndb.QP" |
"E" | "+" | "1" | "08_gw_eq/ndb.QP" |
%     
</pre>
</div>
</div>
where we tell ypp to take from the quasi particle correction databases the values of <math>\varepsilon^{cohsex,neq}_{c\mathbf{k}}</math> subtract from them <math>\varepsilon^{cohsex,eq}_{c\mathbf{k}}</math>, multiply the resulting value by <code>?????</code>  and add it to <math>\varepsilon^{GW,eq}_{c\mathbf{k}}</math> and save them in the the new database <code>ndb.QP_merged_1_gw_cohsex_gw_ppa</code>.
We run the file by the following command:
$ ypp_rt -F input_ypp3 -J QP_merge -C QP_merge
Or alternatively, if no calculation for the quasiparticle correction has been made (except the COHSEX calculation), and a scissors operator is going to be used, we can use the following yambo post processing file:
<div class="mw-collapsible mw-collapsed" style="width:1000px; overflow:auto;">
<div style="font-weight:bold;line-height:1.6;"> input_ypp3_scissors (Expandable) </div>
<div class="mw-collapsible-content">
<pre>
QPDBs                      # [R] Quasi-particle databases
QPDB_merge                # [R] Mergering
%Actions_and_names      # [QPDB] QP databases and Actions (format is "what"|"OP"|"prefactor"|"DB"|). OP can be ±/x.
"E" | "+" | "1" | "cohsex_neq/ndb.QP" |
"E" | "-" | "1" | "cohsex_eq_emp/ndb.QP" |
</pre>
</div>
</div>
which will just take from the quasi particle correction databases the values of <math>\varepsilon^{cohsex,neq}_{c\mathbf{k}}</math> subtract from them <math>\varepsilon^{cohsex,eq}_{c\mathbf{k}}</math>. And we run this file by:
$ ypp_rt -F input_ypp3_scissors -J QP_merge_shissor -C QP_merge_shissor


= Renormalization of the exciton binding energy =
= Renormalization of the exciton binding energy =


We finally perform a BSE calculation loading both the nonequilibrium carriers and the quasi-particle energy shifts.
We finally perform a BSE calculation loading both the nonequilibrium carriers and the quasi-particle energy shifts.
We create a file using the following command:
$ yambo -o b -y d -V qp -F 04_neq_bse.in
== With GW qp corrections ==
where we add the lines:
XfnRTdb= "f @ 0. ps < carriers/ndb.RT_carriers"
KfnRTdb= "f @ 0. ps < carriers/ndb.RT_carriers"
RTOccMode="KR"
KfnQPdb= "E < QP_merge/ndb.QP_merged_1_gw_cohsex_gw_ppa"
the resulting structure of the file is:
<div class="mw-collapsible mw-collapsed" style="width:1000px; overflow:auto;">
<div style="font-weight:bold;line-height:1.6;"> 04_neq_bse.in (Expandable) </div>
<div class="mw-collapsible-content">
<pre>
bss                              # [R] BSE solver
optics                          # [R] Linear Response optical properties
dipoles                          # [R] Oscillator strenghts (or dipoles)
bse                              # [R][BSE] Bethe Salpeter Equation.
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
BSKmod= "SEX"                    # [BSE] IP/Hartree/HF/ALDA/SEX/BSfxc
BSEmod= "resonant"              # [BSE] resonant/retarded/coupling
BSSmod= "d"                      # [BSS] (h)aydock/(d)iagonalization/(s)lepc/(i)nversion/(t)ddft`
BSENGexx= 2085            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= "none"                  # [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.000000 | 1.000000 | 1.000000 |        # [EXTQP BSK BSS] E parameters  (c/v) eV|adim|adim
%
KfnQP_Z= ( 1.000000 , 0.000000 )        # [EXTQP BSK BSS] Z factor  (c/v)
KfnQP_Wv_E= 0.000000      eV    # [EXTQP BSK BSS] W Energy reference  (valence)
% KfnQP_Wv
0.000000 | 0.000000 | 0.000000 |        # [EXTQP BSK BSS] W parameters  (valence) eV| 1|eV^-1
%
KfnQP_Wv_dos= 0.000000    eV    # [EXTQP BSK BSS] W dos pre-factor  (valence)
KfnQP_Wc_E= 0.000000      eV    # [EXTQP BSK BSS] W Energy reference  (conduction)
% KfnQP_Wc
0.000000 | 0.000000 | 0.000000 |        # [EXTQP BSK BSS] W parameters  (conduction) eV| 1 |eV^-1
%
XfnRTdb= "f @ 0. ps < carriers/ndb.RT_carriers"
KfnRTdb= "f @ 0. ps < carriers/ndb.RT_carriers"
KfnQPdb= "E < QP_merge/ndb.QP_merged_1_gw_cohsex_gw_ppa"
RTOccMode="KR"
KfnQP_Wc_dos= 0.000000    eV    # [EXTQP BSK BSS] W dos pre-factor  (conduction)
% BSEQptR
1 | 1 |                            # [BSK] Transferred momenta range
%
% BSEBands
  2 |  7 |                        # [BSK] Bands range
%
% BEnRange
  2.00000 | 6.00000 |        eV    # [BSS] Energy range
%
% BDmRange
0.100000 | 0.100000 |        eV    # [BSS] Damping range
%
BEnSteps= 200                    # [BSS] Energy steps
% BLongDir
1.000000 | 0.000000 | 0.000000 |        # [BSS] [cc] Electric Field versor
%
BSEprop= "abs"                  # [BSS] Can be any among abs/jdos/kerr/asymm/anHAll/magn/dich/photolum/esrt
BSEdips= "none"                  # [BSS] Can be "trace/none" or "xy/xz/yz" to define off-diagonal rotation plane
#WRbsWF                        # [BSS] Write to disk excitonic the WFs
</pre>
</div>
</div>
where we added the lines, and changed the bands range, <code>BSEBands</code> from <code>1|50</code> to <code>2|7</code>, energy range, <code>BEnRange</code> from <code>0.00000 | 10.00000 </code> to <code>2.00000 | 10.00000 |</code> and energy steps, <code>BEnSteps</code> from <code>100</code> to <code>200</code>. And we run the file:
$ yambo -F 04_neq_bse.in -J "neq_bse,gw_eq"
Alternatively, when scissors are used, the lines can be added, to account for the lack of a proper GW correction.
== With shissor corrections ==
where we add the lines:
XfnRTdb= "f @ 0. ps < carriers/ndb.RT_carriers"
KfnRTdb= "f @ 0. ps < carriers/ndb.RT_carriers"
RTOccMode="KR"
% KfnQP_E
  0.70000 | 1.000000 | 1.000000 |      # [EXTQP BSK BSS] E parameters (c/v)
%
KfnQPdb= "E < QP_merge_shissor/ndb.QP_merged_1_gw_cohsex"
the resulting structure of the file is:
<div class="mw-collapsible mw-collapsed" style="width:1000px; overflow:auto;">
<div style="font-weight:bold;line-height:1.6;"> 04_neq_bse.in (Expandable) </div>
<div class="mw-collapsible-content">
<pre>
bss                              # [R] BSE solver
optics                          # [R] Linear Response optical properties
dipoles                          # [R] Oscillator strenghts (or dipoles)
bse                              # [R][BSE] Bethe Salpeter Equation.
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
BSKmod= "SEX"                    # [BSE] IP/Hartree/HF/ALDA/SEX/BSfxc
BSEmod= "resonant"              # [BSE] resonant/retarded/coupling
BSSmod= "d"                      # [BSS] (h)aydock/(d)iagonalization/(s)lepc/(i)nversion/(t)ddft`
BSENGexx= 2085            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
% KfnQP_E
0.700000 | 1.000000 | 1.000000 |        # [EXTQP BSK BSS] E parameters  (c/v) eV|adim|adim
%
XfnRTdb= "f @ 0. ps < carriers/ndb.RT_carriers"
KfnRTdb= "f @ 0. ps < carriers/ndb.RT_carriers"
KfnQPdb= "E < QP_merge/ndb.QP_merged_1_gw_cohsex"
RTOccMode="KR"
% BSEQptR
1 | 1 |                            # [BSK] Transferred momenta range
%
% BSEBands
  2 |  7 |                        # [BSK] Bands range
%
% BEnRange
  2.00000 | 6.00000 |        eV    # [BSS] Energy range
%
% BDmRange
0.100000 | 0.100000 |        eV    # [BSS] Damping range
%
BEnSteps= 200                    # [BSS] Energy steps
% BLongDir
1.000000 | 0.000000 | 0.000000 |        # [BSS] [cc] Electric Field versor
%
BSEprop= "abs"                  # [BSS] Can be any among abs/jdos/kerr/asymm/anHAll/magn/dich/photolum/esrt
BSEdips= "none"                  # [BSS] Can be "trace/none" or "xy/xz/yz" to define off-diagonal rotation plane
#WRbsWF                        # [BSS] Write to disk excitonic the WFs
</pre>
</div>
</div>
where we added the lines, and changed the bands range, <code>BSEBands</code> from <code>1|50</code> to <code>2|7</code>, energy range, <code>BEnRange</code> from <code>0.00000 | 10.00000 </code> to <code>2.00000 | 10.00000 |</code> and energy steps, <code>BEnSteps</code> from <code>100</code> to <code>200</code>. And we run the file:


= References =
= References =

Latest revision as of 15:56, 21 November 2024

In this tutorial you will learn the basic concepts for computing changes in the optical properties of a semi-conductor in presence of a non-equilibrium electrons and holes distribution in conduction and valence band respectively. This tutorial is based on the results published in Phys. Rev. B[1]

Under construction

The material: Silicon

We will study nonequilibrium absorption in bulk silicon. The same material used for this GW tutorial

  • FCC lattice
  • Two atoms per cell (8 electrons)
  • Lattice constant 10.183 [a.u.]
  • Plane waves cutoff 15 Rydberg
  • Direct gap 3.4 eV at Gamma
  • Indirect gap 1.1 eV between Gamma= (0 0 0) and a point X', close to X=(0 1 0)
Silicon Band Structure

Tutorial files and Tutorial structure

Follow the instructions in Tutorials#Files and download/unpack the Silicon.tar.gz. Once the tutorial archive file is unzipped the following folder structure will appear

COPYING  README  Silicon/

with the Solid_Si folder containing

> ls Silicon/ 
PWSCF/  YAMBO/

In the Pwscf folder the student will find an input/output directory with input/output files for pw.x. The Silicon pseudopotential file is also provided.

> ls PWSCF/
convergence_scripts  input  output  psps

In the convergence_scripts you will find some useful shell scripts to run the ground state convergence runs for Silicon. The YAMBO folder contains the Yambo input/output files and core databases.

> ls YAMBO/
2x2x2/  4x4x4/  6x6x6/  8x8x8/  Convergence_Plots_and_Scripts/  GAMMA/

The core databases are provided for several k-points grids. In addition the folder Convergence_Plots_and_Scripts contains some scripts used for the [Silicon|GW tutorial] . Here we will just use the 8x8x8 (which is still very far from convergence) folder for computing (nonequilibrium) optical properties.

To run the tutorial you will need the standard executables

yambo
ypp 

plus the executables of the real time module of the Yambo code

yambo_rt
ypp_rt

Equilibrium optical properties

We start from a density-functional theory (DFT) calculation. Then we apply MBPT to obtain the QP energies within the GW approximation and the absorption spectra from the solution of the BSE. The QP band structure will show the effect of the screened electron-electron interaction, while the absorption spectrum within the BSE accounts for the effect of the electron-hole interaction. When QP corrections are not available, they can be included using the scissor operator KfnQP_E, and the numerical values for "corrections" can be deduced from a GW calculation. For more information about solving the Equilibrium BSE equation refer to the following tutorial: Calculating optical spectra including excitonic effects: a step-by-step guide

(to delete open) In this case, do we actually compute the qp eng withing gw approximation? or we directy solve the bse eq (to delete close)

Enter the folder 8x8x8, and initialize yambo in the respective folder

$ cd 8x8x8
$ yambo

For this step you can either compute static screening at equilibrium, or use the screening computed for the GW step in the Bethe-Salpeter.

Static screening at equilibrium

First, to solve the BSE equation we need to calculate the statically screened electron-electron interaction (W), which is related to static diaelectric screening [math]\displaystyle{ \epsilon^{-1} }[/math].

[math]\displaystyle{ \Large \epsilon^{-1}_{\mathbf{G},\mathbf{G'}} (\mathbf{q}, \omega = 0) = \delta_{\mathbf{G},\mathbf{G'}} + v(\mathbf{v} + \mathbf{G}) \chi_{\mathbf{G},\mathbf{G'}(\mathbf{q}}, \omega = 0) }[/math]

To do this create the following input file:

$ yambo -X s -F 01_EQ_BSE_screening.in

with the following contents:

01_EQ_BSE_screening.in file (Expandable)
em1s                                 # [R][Xs] Statically Screened Interaction
% BndsRnXs
   1 |  30 |                         # [Xs] Polarization function bands
%
NGsBlkXs= 3                    Ry    # [Xs] Response block size
% LongDrXs
 1.000000 | 0.000000 | 0.000000 |    # [Xs] [cc] Electric Field
%

Where we changed the response block size NGsBlkXs = 3 Ry and number of polarization function bands BndsRnXs to be in the range of 1 to 30. Now run the file:

$ yambo -F 01_EQ_BSE_screening.in -J EQ_BSE

Solving the Bethe-Salpeter equation

To solve the Bethe-Salpeter equation, the latter is usually rewritten in the space of transitions between valence and conduction states as the (pseudo)eigenvalue problem for a two-particle Hamiltonian. For a non spin-polarized system and in the [math]\displaystyle{ q = 0 }[/math] limit, the two-particle Hamiltonian matrix elements are given by

[math]\displaystyle{ \Large H_{vc\mathbf{k},v'c'\mathbf{k'}} = (\varepsilon_{c\mathbf{k}} - \varepsilon_{v\mathbf{k}})\delta_{vv'}\delta_{cc'}\delta_{\mathbf{k},\mathbf{k'}} + (f_{c\mathbf{k}} - f_{v\mathbf{k}})\left[ 2V_{vc\mathbf{k},v'c'\mathbf{k'}} - W_{vc\mathbf{k},v'c'\mathbf{k'}}\right] }[/math]

where [math]\displaystyle{ vc\mathbf{k} }[/math] indicates the pair of quasiparticle states [math]\displaystyle{ v\mathbf{k} }[/math] and [math]\displaystyle{ c\mathbf{k} }[/math]. The first term on the RHS is the quasiparticle energy differences (diagonal only). The second term is the kernel which is the sum of the electron-hole exchange part V (which stems from the Hartree potential) and the electron-hole attraction part W (which stems from the screened exchange potential). The electron-hole attraction part W, depends on [math]\displaystyle{ \epsilon^{-1} }[/math] which was calculated in the previous section.

After calculating the equilibrium Bethe-Salpeter kernel, the resulting two-particle Hamiltonian can be diagonalized to compute the macroscopic dielectric function, which can be used to get absorption and EELS spectra. Calculation of BS kernel and diagonalization of Hamiltonian can be done in one step.

To generate the input file for the calculation of BSE kernel, run the following command:

$ yambo -o b -k sex -F 02_EQ_BSE_kernel_diag.in -y d -J EQ_BSE

with the following input

02_EQ_BSE_kernel_diag.in file (Expandable)
optics                           # [R] Linear Response optical properties
bss                              # [R] BSE solver
bse                              # [R][BSE] Bethe Salpeter Equation.
K_Threads=0                      # [OPENMP/BSK] Number of threads for response functions
BSKmod= "SEX"                    # [BSE] IP/Hartree/HF/ALDA/SEX/BSfxc
BSEmod= "resonant"               # [BSE] resonant/retarded/coupling
BSSmod= "d"                      # [BSS] (h)aydock/(d)iagonalization/(s)lepc/(i)nversion/(t)ddft`
BSENGexx= 2085             RL    # [BSK] Exchange components
BSENGBlk=-1                RL    # [BSK] Screened interaction block size [if -1 uses all the G-vectors of W(q,G,Gp)]
% BSEQptR
 1 | 1 |                             # [BSK] Transferred momenta range
%
% BSEBands
   2 |  7 |                         # [BSK] Bands range
%
% BEnRange
  2.20000 | 5.00000 |         eV    # [BSS] Energy range
%
% KfnQP_E
 0.70000 | 1.000000 | 1.000000 |      # [EXTQP BSK BSS] E parameters (c/v)
%
% BDmRange
 0.100000 | 0.100000 |         eV    # [BSS] Damping range
%
BEnSteps= 200                    # [BSS] Energy steps
% BLongDir
 1.000000 | 0.000000 | 0.000000 |        # [BSS] [cc] Electric Field versor
%
BSEprop= "abs"                   # [BSS] Can be any among abs/jdos/kerr/asymm/anHAll/magn/dich/photolum/esrt
BSEdips= "none"                  # [BSS] Can be "trace/none" or "xy/xz/yz" to define off-diagonal rotation plane

Where BSEBands was changed to include bands from 3 to 8, which defines the basis of quasiparticle pairs that are used to describe excitons. 2 valence and 4 conduction bands where included. BEnRange was increased to from 0-10 eV to 0-16 eV, and BEnSteps was increased from 100 to 200. The scissors operator KfnQP_E was used to match the positon calculated spectra with experimental values.

To run the file:

$ yambo -F 02_EQ_BSE_kernel_diag.in -J EQ_BSE

After the calculation has finished, the absorption spectra can be plotted with the following:

gnuplot commands (Expandable)
$ gnuplot
...
gnuplot> set x label "Energy (eV)"
gnuplot> set y label "Absorption" 
gnuplot> plot 'o-EQ_BSE.eps_q1_diago_bse' u 1:2 w l t '-k sex'

and the following plot will results (the "experimental" and "-k hartree" plots are added additionally for comparison)

Plot of the computed spectra (Expandable)

Plots of computed absorption spectra for the equilibrium case.

Generating non-equilibrium carriers

To simulate non-equilibrium absorption, we need to define the number of non-equilibrium carriers, specifically electrons in the conduction band and holes in the valence band. These carriers follow the Fermi-Dirac distribution given by:

[math]\displaystyle{ \Large f(E) = 1 / (exp((E - μ) / k_B T) + 1) }[/math]

Here:

  • E is the energy level,
  • [math]\displaystyle{ \mu }[/math] is the chemical potential (specific to electrons or holes),
  • [math]\displaystyle{ T }[/math] is the effective temperature,
  • [math]\displaystyle{ k_B }[/math] is the Boltzmann constant.

For the simulation, both electrons and holes are characterized by their own [math]\displaystyle{ \mu_i }[/math] and [math]\displaystyle{ T_i }[/math] .

The distribution of non-equilibrium carriers can be generated using Yambo’s post processing for real-time module (ypp_rt), which allows modeling of these conditions.

To generate to non-equilibrium carriers databases create input_ypp1 input_ypp2 files, with the following content:

input_ypp1 (Expandable)
RTDBs                      	# [R] Real-Time databases
Select_Fermi               	# [R] NEQ DBs according to Fermi distribution
% RTBands
   1 | 8 |               	# [RT] Bands range
%
% RTmuEh
-0.80 | 0.80 | eV  	        # [RT] Chemical potentials hole | electron
%
% RTtempEh
 0.04 	| 0.04 	| eV  	    # [RT] Effective temperature hole | electron
%
RTautotuneThr= 0.000100    	# [RT] Threshold to match no. of pumped holes and electrons.
input_ypp2 (Expandable)
RealTime                   	# [R] Real-Time Post-Processing
RToccupations              	# [R] Analize time-dependent occupations
RTenergy                   	# [R] Post-Processing kind: function of energy
TimeStep= 0.000000 	fs  	# Time step
% TimeRange
-1.000000 |-1.000000 | fs  	# Time-window where processing is done
%
#IncludeEQocc              	# [NEGF] Include also equilibrium occupations
#SkipFermiFIT              	# [NEGF] Do a Fermi Fit of occupations (and lifetimes ratio)
%QPkrange                	#  generalized Kpoint/Band indices
  1| 3|  1|8|
%

And then run the files:

$ ypp_rt 
$ ypp_rt -F input_ypp1 -J carriers -C carriers 
$ ypp_rt -F input_ypp2 -J carriers -C carriers 

The file input_ypp1 will generate a ndb.RT_carriers database which will contain the Fermi distribution of the electrons and the hole according to the specified chemical potential, RTmuEh, and for a specific temperature given by RTtempEh. The file input_ypp2 will generate a series of plots from the carriers database.

Energy shift in the band structure

The GW method is the standard approach to compute quasi-particle corrections. However, in presence of non-equilibrium carriers, or even at finite temperature the formulation of the GW self-energy needs to be refined. This is due to its frequency dependence, which results in extra terms when performing the analytic continuation from the Keldish contour (non-equilibrium case) or from the Matsubara axis (finite temperature case) to the real time/frequency axis. See for example this reference [2] for the finite temperature case. On the other hand the COHSEX self-energy, being static, avoids this complication. This is why we will compute changes in the QP corrections within the COHSEX approximation.

To compute the Non-Equilibrium Absorption Spectra, we need to diagonalize a previously defined two-particle Hamiltonian. To account for the Non-Equilibrium effect, the Kohn-Sham energies, [math]\displaystyle{ \varepsilon_{c\mathbf{k}}, \varepsilon_{v\mathbf{k}} }[/math], calculated at the DFT level, need to be corrected. This requires calculating the following: COHSEX self-energy at equilibrium, COHSEX self-energy at non-equilibrium, and Equilibrium GW self-energy at Plasmon-Pole Approximation level.

Using these, the energy correction, [math]\displaystyle{ \varepsilon^{corr,neq} }[/math], can be computed as:

[math]\displaystyle{ \Large \varepsilon^{corr,neq}_{c\mathbf{k}} = \varepsilon^{GW,eq}_{c\mathbf{k}} + (\varepsilon^{cohsex,neq}_{c\mathbf{k}} - \varepsilon^{cohsex,eq}_{c\mathbf{k}}) }[/math]

The KS energies will be corrected as follows:

[math]\displaystyle{ \Large \varepsilon^{KS,neq}_{c\mathbf{k}} = \varepsilon_{c\mathbf{k}} + \varepsilon^{corr,neq}_{c\mathbf{k}} }[/math]

These corrected KS energies will be used when diagonalizing the two-particle Hamiltonian.

COHSEX corrections at equilibrium

To make a COHSEX calculation at equilibrium, we do the following:

$ yambo -p c -F 03_eq_cohsex.in 

which will create a file with the following content:

03_eq_cohsex.in (Expandable)
gw0                              # [R] GW approximation
cohsex                           # [R][Xp] COlumb Hole Screened EXchange
dyson                            # [R] Dyson Equation solver
el_el_corr                       # [R] Electron-Electron Correlation
HF_and_locXC                     # [R] Hartree-Fock
em1s                             # [R][Xs] Statically Screened Interaction
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
SE_Threads=0                     # [OPENMP/GW] Number of threads for self-energy
EXXRLvcs= 2085             RL    # [XX] Exchange    RL components
VXCRLvcs= 2085             RL    # [XC] XCpotential RL components
Chimod= "HARTREE"                # [X] IP/Hartree/ALDA/LRC/PF/BSfxc
% BndsRnXs
   1 |  10 |                         # [Xs] Polarization function bands
%
NGsBlkXs= 1                RL    # [Xs] Response block size
% LongDrXs
 1.000000 | 0.000000 | 0.000000 |        # [Xs] [cc] Electric Field
%
%QPkrange                        # [GW] QP generalized Kpoint/Band indices
1|8|1|10|
%

Where the number of Polarization function bands, BndsRnXs was changed from 1|50 to 1|10, and the quasi-particles generalized Kpoint and Band indices,QPkrange, was changed from 1|29|1|50 to 1|8|1|10. And after it we run the file with the following command:

$ yambo_rt -F 03_eq_cohsex.in -J cohsex_eq_emp -C cohsex_eq_emp

COHSEX corrections in presence of non-equilibrium carriers

Before making a COHSEX calculation, we need to recalculate the static screening in the presence of non-equilibrium carriers. After it, we can make a COHSEX calculation in the presence of non-equilibrium carriers using the ndb.RT_carriers database generated in the previous step with the ypp_rt module. To do this, we do the following:

$ yambo -p c -F 03_neq_cohsex.in 

where we need to add the following lines:

XfnRTdb= "f @ 0. ps < carriers/ndb.RT_carriers"
GfnRTdb= "f @ 0. ps < carriers/ndb.RT_carriers"

and the resulting file is:

03_neq_cohsex.in (Expandable)
gw0                              # [R] GW approximation
cohsex                           # [R][Xp] COlumb Hole Screened EXchange
dyson                            # [R] Dyson Equation solver
el_el_corr                       # [R] Electron-Electron Correlation
HF_and_locXC                     # [R] Hartree-Fock
em1s                             # [R][Xs] Statically Screened Interaction
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
SE_Threads=0                     # [OPENMP/GW] Number of threads for self-energy
EXXRLvcs= 2085             RL    # [XX] Exchange    RL components
VXCRLvcs= 2085             RL    # [XC] XCpotential RL components
Chimod= "HARTREE"                # [X] IP/Hartree/ALDA/LRC/PF/BSfxc
% BndsRnXs
   1 |  10 |                         # [Xs] Polarization function bands
%
NGsBlkXs= 1                RL    # [Xs] Response block size
% LongDrXs
 1.000000 | 0.000000 | 0.000000 |        # [Xs] [cc] Electric Field
%
%QPkrange                        # [GW] QP generalized Kpoint/Band indices
1|8|1|10|
%

where we added the respective two lines with XfnRTdb, and GfnRTdb, and made changes to BndsRnXs and QPkrange similar to the equilibrium case.

The line

> em1s                             # [R][Xs] Statically Screened Interaction

will tell Yambo to make the static screening calculation, and the lines:

> gw0                              # [R] GW approximation
> cohsex                           # [R][Xp] COlumb Hole Screened EXchange

will tell Yambo to do make a COHSEX calculation.

After it we can run the file with the command:

$ yambo_rt -F 03_neq_cohsex.in -J cohsex_neq -C cohsex_neq

GW PPA corrections for the equilibrium case

This subchapter can be skipped, and replaced by using scissors operator when diagonalizing the two particle Hamiltonian.

To get accurate values for the band gap value we need to go beyond COHSEX approximation. Therefore we will obtain the Quasi-particle corrections in the plasmon pole approximation for equilibrium case, and in next steps, to this value we will add the correction calculated at COHSEX approximation to account for the non-equilibrium state of the system.

To make this calculation, we do the following:

$ yambo -p p -F 03_eq_qp_ppa.in 

which will create a file with the following content:

03_eq_qp_ppa.in (Expandable)
gw0                              # [R] GW approximation
ppa                              # [R][Xp] Plasmon Pole Approximation for the Screened Interaction
el_el_corr                       # [R] Electron-Electron Correlation
HF_and_locXC                     # [R] Hartree-Fock
em1d                             # [R][X] Dynamically Screened Interaction
X_Threads=0                      # [OPENMP/X] Number of threads for response functions
DIP_Threads=0                    # [OPENMP/X] Number of threads for dipoles
SE_Threads=0                     # [OPENMP/GW] Number of threads for self-energy
EXXRLvcs= 2085             RL    # [XX] Exchange    RL components
VXCRLvcs= 2085             RL    # [XC] XCpotential RL components
Chimod= "HARTREE"                # [X] IP/Hartree/ALDA/LRC/PF/BSfxc
% BndsRnXp
   1 |  10 |                         # [Xp] Polarization function bands
%
NGsBlkXp= 1                RL    # [Xp] Response block size
% LongDrXp
 1.000000 | 0.000000 | 0.000000 |        # [Xp] [cc] Electric Field
%
PPAPntXp= 27.21138         eV    # [Xp] PPA imaginary energy
% GbndRnge
   1 |  10 |                         # [GW] G[W] bands range
%
GTermKind= "none"                # [GW] GW terminator ("none","BG" Bruneval-Gonze,"BRS" Berger-Reining-Sottile)
DysSolver= "n"                   # [GW] Dyson Equation solver ("n","s","g","q")
%QPkrange                        # [GW] QP generalized Kpoint/Band indices
1|8|1|10|
%        

where the following changes have been made: the number of polarization function bands, BndsRnXp was changed from 1|50 to 1|10, the G[W] bands range, GbndRnge was changed from 1|50 to 1|10, and quasi-particle generalized k points and band indices, QPkrange was changed from 1|29|1|50 to 1|8|1|10. And after it we run the file:

$ yambo_rt -F 03_eq_qp_ppa.in -J gw_eq -C gw_eq

Data post-processing

Now after obtaining all the energy components that we need to calculate corrections to Kohn-Sham energies for the Non-Equilibrium case, we can manipulate them using the yambo's real time post processing module, ypp_rt. We want to obtain energies defined by the following equation:

[math]\displaystyle{ \Large \varepsilon^{corr,neq}_{c\mathbf{k}} = \varepsilon^{GW,eq}_{c\mathbf{k}} + (\varepsilon^{cohsex,neq}_{c\mathbf{k}} - \varepsilon^{cohsex,eq}_{c\mathbf{k}}) }[/math].

Therefore we create a yambo post processing file with the following contents:

input_ypp3 (Expandable)
QPDBs                      	# [R] Quasi-particle databases
QPDB_merge                 	# [R] Mergering
%Actions_and_names       	# [QPDB] QP databases and Actions (format is "what"|"OP"|"prefactor"|"DB"|). OP can be ±/x.
 "E" | "+" | "1" | "cohsex_neq/ndb.QP" |
 "E" | "-" | "1" | "cohsex_eq_emp/ndb.QP" |
 "Z" | "+" | "1" | "gw_eq/ndb.QP" |
 "E" | "+" | "1" | "08_gw_eq/ndb.QP" |
%       

where we tell ypp to take from the quasi particle correction databases the values of [math]\displaystyle{ \varepsilon^{cohsex,neq}_{c\mathbf{k}} }[/math] subtract from them [math]\displaystyle{ \varepsilon^{cohsex,eq}_{c\mathbf{k}} }[/math], multiply the resulting value by ????? and add it to [math]\displaystyle{ \varepsilon^{GW,eq}_{c\mathbf{k}} }[/math] and save them in the the new database ndb.QP_merged_1_gw_cohsex_gw_ppa.

We run the file by the following command:

$ ypp_rt -F input_ypp3 -J QP_merge -C QP_merge

Or alternatively, if no calculation for the quasiparticle correction has been made (except the COHSEX calculation), and a scissors operator is going to be used, we can use the following yambo post processing file:

input_ypp3_scissors (Expandable)
QPDBs                      	# [R] Quasi-particle databases
QPDB_merge                 	# [R] Mergering
%Actions_and_names       	# [QPDB] QP databases and Actions (format is "what"|"OP"|"prefactor"|"DB"|). OP can be ±/x.
 "E" | "+" | "1" | "cohsex_neq/ndb.QP" |
 "E" | "-" | "1" | "cohsex_eq_emp/ndb.QP" |
%  

which will just take from the quasi particle correction databases the values of [math]\displaystyle{ \varepsilon^{cohsex,neq}_{c\mathbf{k}} }[/math] subtract from them [math]\displaystyle{ \varepsilon^{cohsex,eq}_{c\mathbf{k}} }[/math]. And we run this file by:

$ ypp_rt -F input_ypp3_scissors -J QP_merge_shissor -C QP_merge_shissor

Renormalization of the exciton binding energy

We finally perform a BSE calculation loading both the nonequilibrium carriers and the quasi-particle energy shifts.

We create a file using the following command:

$ yambo -o b -y d -V qp -F 04_neq_bse.in

With GW qp corrections

where we add the lines:

XfnRTdb= "f @ 0. ps < carriers/ndb.RT_carriers"
KfnRTdb= "f @ 0. ps < carriers/ndb.RT_carriers"
RTOccMode="KR"
KfnQPdb= "E < QP_merge/ndb.QP_merged_1_gw_cohsex_gw_ppa"

the resulting structure of the file is:

04_neq_bse.in (Expandable)
bss                              # [R] BSE solver
optics                           # [R] Linear Response optical properties
dipoles                          # [R] Oscillator strenghts (or dipoles)
bse                              # [R][BSE] Bethe Salpeter Equation.
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
BSKmod= "SEX"                    # [BSE] IP/Hartree/HF/ALDA/SEX/BSfxc
BSEmod= "resonant"               # [BSE] resonant/retarded/coupling
BSSmod= "d"                      # [BSS] (h)aydock/(d)iagonalization/(s)lepc/(i)nversion/(t)ddft`
BSENGexx= 2085             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= "none"                  # [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.000000 | 1.000000 | 1.000000 |        # [EXTQP BSK BSS] E parameters  (c/v) eV|adim|adim
%
KfnQP_Z= ( 1.000000 , 0.000000 )         # [EXTQP BSK BSS] Z factor  (c/v)
KfnQP_Wv_E= 0.000000       eV    # [EXTQP BSK BSS] W Energy reference  (valence)
% KfnQP_Wv
 0.000000 | 0.000000 | 0.000000 |        # [EXTQP BSK BSS] W parameters  (valence) eV| 1|eV^-1
%
KfnQP_Wv_dos= 0.000000     eV    # [EXTQP BSK BSS] W dos pre-factor  (valence)
KfnQP_Wc_E= 0.000000       eV    # [EXTQP BSK BSS] W Energy reference  (conduction)
% KfnQP_Wc
 0.000000 | 0.000000 | 0.000000 |        # [EXTQP BSK BSS] W parameters  (conduction) eV| 1 |eV^-1
%
XfnRTdb= "f @ 0. ps < carriers/ndb.RT_carriers"
KfnRTdb= "f @ 0. ps < carriers/ndb.RT_carriers"
KfnQPdb= "E < QP_merge/ndb.QP_merged_1_gw_cohsex_gw_ppa"
RTOccMode="KR"
KfnQP_Wc_dos= 0.000000     eV    # [EXTQP BSK BSS] W dos pre-factor  (conduction)
% BSEQptR
 1 | 1 |                             # [BSK] Transferred momenta range
%
% BSEBands
   2 |  7 |                         # [BSK] Bands range
%
% BEnRange
  2.00000 | 6.00000 |         eV    # [BSS] Energy range
%
% BDmRange
 0.100000 | 0.100000 |         eV    # [BSS] Damping range
%
BEnSteps= 200                    # [BSS] Energy steps
% BLongDir
 1.000000 | 0.000000 | 0.000000 |        # [BSS] [cc] Electric Field versor
%
BSEprop= "abs"                   # [BSS] Can be any among abs/jdos/kerr/asymm/anHAll/magn/dich/photolum/esrt
BSEdips= "none"                  # [BSS] Can be "trace/none" or "xy/xz/yz" to define off-diagonal rotation plane
#WRbsWF                        # [BSS] Write to disk excitonic the WFs

where we added the lines, and changed the bands range, BSEBands from 1|50 to 2|7, energy range, BEnRange from 0.00000 | 10.00000 to 2.00000 | 10.00000 | and energy steps, BEnSteps from 100 to 200. And we run the file:

$ yambo -F 04_neq_bse.in -J "neq_bse,gw_eq"

Alternatively, when scissors are used, the lines can be added, to account for the lack of a proper GW correction.

With shissor corrections

where we add the lines:

XfnRTdb= "f @ 0. ps < carriers/ndb.RT_carriers"
KfnRTdb= "f @ 0. ps < carriers/ndb.RT_carriers"
RTOccMode="KR"
% KfnQP_E
 0.70000 | 1.000000 | 1.000000 |      # [EXTQP BSK BSS] E parameters (c/v)
%
KfnQPdb= "E < QP_merge_shissor/ndb.QP_merged_1_gw_cohsex"

the resulting structure of the file is:

04_neq_bse.in (Expandable)
bss                              # [R] BSE solver
optics                           # [R] Linear Response optical properties
dipoles                          # [R] Oscillator strenghts (or dipoles)
bse                              # [R][BSE] Bethe Salpeter Equation.
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
BSKmod= "SEX"                    # [BSE] IP/Hartree/HF/ALDA/SEX/BSfxc
BSEmod= "resonant"               # [BSE] resonant/retarded/coupling
BSSmod= "d"                      # [BSS] (h)aydock/(d)iagonalization/(s)lepc/(i)nversion/(t)ddft`
BSENGexx= 2085             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
% KfnQP_E
 0.700000 | 1.000000 | 1.000000 |        # [EXTQP BSK BSS] E parameters  (c/v) eV|adim|adim
%
XfnRTdb= "f @ 0. ps < carriers/ndb.RT_carriers"
KfnRTdb= "f @ 0. ps < carriers/ndb.RT_carriers"
KfnQPdb= "E < QP_merge/ndb.QP_merged_1_gw_cohsex"
RTOccMode="KR"
% BSEQptR
 1 | 1 |                             # [BSK] Transferred momenta range
%
% BSEBands
   2 |  7 |                         # [BSK] Bands range
%
% BEnRange
  2.00000 | 6.00000 |         eV    # [BSS] Energy range
%
% BDmRange
 0.100000 | 0.100000 |         eV    # [BSS] Damping range
%
BEnSteps= 200                    # [BSS] Energy steps
% BLongDir
 1.000000 | 0.000000 | 0.000000 |        # [BSS] [cc] Electric Field versor
%
BSEprop= "abs"                   # [BSS] Can be any among abs/jdos/kerr/asymm/anHAll/magn/dich/photolum/esrt
BSEdips= "none"                  # [BSS] Can be "trace/none" or "xy/xz/yz" to define off-diagonal rotation plane
#WRbsWF                        # [BSS] Write to disk excitonic the WFs

where we added the lines, and changed the bands range, BSEBands from 1|50 to 2|7, energy range, BEnRange from 0.00000 | 10.00000 to 2.00000 | 10.00000 | and energy steps, BEnSteps from 100 to 200. And we run the file:

References