The Yambo Project - User contributions [en]

Tutorials all

2025-10-13T19:49:16Z

Yambowikiadmin: /* Non-Linear response */

This page contains the full list of available tutorials. These cover a variety of mixed topics, both physical and methodological. Some of these are om "Standalone mode" and are designed to be followed from start to finish in one page. Others have been ported in "Modular form". These are also available under the link [[Tutorials modular]]. The tutorials do not require previous knowledge of yambo, unless explicitly specified. Each tutorial requires download of a specific core database, and typically they cover a specific physical system (like bulk GaSb or a hydrogen chain). Ground state input files and pseudopotentials are provided. Output files are also provided for reference.

These tutorials can be accessed directly from this page of from the side bar. They include different kind of subjects:

'''Warning''': These tutorials were prepared using previous version of the Yambo code: some command lines, variables, reports and outputs can be slightly different from the last version of the code. Scripts for parsing output cannot work anymore and should be edited to work with the new outputs. New command lines can be accessed typing <code>yambo -h </code>

== Basic ==
* [[First steps: a walk through from DFT to optical properties]] (modular tutorial)
* [[Silicon|GW on bulk silicon]]
* [[GW on h-BN (standalone)|GW on h-BN]] (standalone version); [[How to obtain the quasi-particle band structure of a bulk material: h-BN|GW on h-BN]] (modular version)
* [[LiF|BSE in lithium-floride]]
* [[BSE tutorial on hBN]]
* [[How to analyse excitons]]

== More on GW and Quasi-particles ==
* [[Real_Axis_and_Lifetimes|Real Axis and Lifetimes]]
* [[Self-consistent GW on eigenvalues only]]
* [[GW parallel strategies|GW parallel strategies v1]]; [[GW parallel strategies CECAM|GW parallel strategies v2]]; [[Pushing convergence in parallel|GW convergence in parallel]] (modular tutorials)
* [[GW tutorial on HPC]] (external tutorial on HackMD)
* [[Quasi-particles and Self-energy within the Multipole Approximation (MPA)]]
* [http://www.attaccalite.com/gw-correction-on-an-arbitrary-point-of-the-brillouin-zone GW on arbitrary k point] (external tutorial on www.attaccalite.com)

== More on Linear response and BSE ==
* [[Hydrogen chain|TDDFT Failure and long range correlations]]
* [[SOC|Spin-orbit coupling in GaSb]]
* [[The magneto-optical Kerr effect|The magneto-optical Kerr effect in bulk iron, RPA]]
* [[The magneto-optical Kerr effect including excitonic effects|The magneto-optical Kerr effect in XXX, excitonic effects]] (todo)
* [[Dichroism in molecules]]
* [[Surface spectroscopy]] (to restore from older version)
* [[Si_Surface|Linear Response of a silicon surface (2D)]]
* [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] (to be imported)
* [[Si_wire|Linear Response from a silicon wire (1D)]] (to restore from old version)
* [[H2|H2 molecule: Linear Response & TDDFT (0D)]] (to restore from old version)
* [[SiH4|SiH4: isolated molecule TDDFT & BSE (0D)]] (restored from old version, still to be finalized)

== Post Processing with ypp ==
* [[Yambo Post Processing (ypp)]]

== Electron phonon coupling ==
* [[Electron Phonon Coupling|Electron Phonon Coupling]]
* [[Optical properties at finite temperature]]
* [[Phonon-assisted luminescence by finite atomic displacements]]
* [[Exciton-phonon coupling and luminescence]]
* [[Thermal lines and special displacements with YamboPy]]

== Real time & Non linear response ==
=== Linear response ===
* [[Linear response from real time simulations (density matrix only)|Linear response using the density matrix (IP)]]
* [[Linear response using Dynamical Berry Phase|Linear response using wave-functions and Dynamical Berry Phase (IP)]]
* [[Linear response in velocity gauge|Linear response using wave-functions in velocity gauge (IP)]]
* [[Real time Bethe-Salpeter Equation (density_matrix_only)|Linear response using the density matrix (TD-HSEX)]]
* [[Linear response from real time simulations|Linear response from real time simulations (IP & TD-HSEX)]] (modular tutorial)

=== Non-Linear response ===
* [[Real time approach to non-linear response (SHG)]]
* [[Correlation effects in the non-linear response]]
* [[SHG within the TD-HSEX level (also called TD-aGW or TD-BSE)]]
* [[A fast approach to excitonic effects in linear/non-linear response]]
* [[Angular dependence of non-linear response]]
* [[Third Harmonic Generation (THG)]]
* [http://www.attaccalite.com/lumen/spin_orbit.html Spin-orbit coupling and non-linear response] (external tutorial on www.attaccalite.com)
* [[Two-photon absorption]]
* [[Sum frequency generation]]
* [[Shift current]]
* [[Parallelization for non-linear response calculations]]

=== Pump and Probe experiments ===
* [[Nonequilibrium absorption in bulk silicon|Pump and Probe: Transient Absorption from BSE with NEQ occupations]]
* [[Pump and Probe|Pump and Probe: Transient Absorption from real time propagation with pump and probe]]
* [[Pump and Probe: Time resolved ARPES|Pump and Probe: Time resolved ARPES from real time propagation with pump]]
* [[Pump and Probe: reading electronic occuations from Perturbo code]]

== Developing Yambo ==
* [[How to create a new project in Yambo]]
* [[How to create a new ypp interface]]
* [[Some hints on github]]

SHG within the TD-aGW level (also called TD-HSEX, TD-BSE)

2025-10-13T19:48:16Z

Yambowikiadmin: Yambowikiadmin moved page SHG within the TD-aGW level (also called TD-HSEX, TD-BSE) to SHG within the TD-HSEX level (also called TD-aGW or TD-BSE)

#REDIRECT [[SHG within the TD-HSEX level (also called TD-aGW or TD-BSE)]]

SHG within the TD-HSEX level (also called TD-aGW or TD-BSE)

2025-10-13T19:48:16Z

Yambowikiadmin: Yambowikiadmin moved page SHG within the TD-aGW level (also called TD-HSEX, TD-BSE) to SHG within the TD-HSEX level (also called TD-aGW or TD-BSE)

[[File:TD-aGW.png|right|300px |Sum frequency generation]]

In this tutorial we show step by step how to calculate SHG within the TD-aGW approximation.<ref name="Attaccalite2013"/>
This approximation was also called TD-BSE or TD-HSEX in previous papers.<ref name="Attaccalite2011"/>

This calculation is quite complicated, here we present a simple strategy to do it.
We will consider hBN monolayer as example because calculations are very fast.

== Code ==
Download the code from github

[https://github.com/yambo-code/yambo https://github.com/yambo-code/yambo]

and compiler all yambo projects by doing

./configure
make all

because this tutorial will require different yambo executables

== Calculations for the quasi-particle band structure ==

We start from a self-consistent(scf) and a non-self-consistent (nscf) of hBN with QuantumEspresso.
All inputs file are here: [https://www.attaccalite.com/tutorials_yambo/hBNsymm.tgz hBNsymm.tgz]

We included 40 bands in the DFT calculation that are used to converge the dielectric constant the enters in the GW and in the COLLISIONS calculations. 
We run all the calculation in the folder hBNsymm.

* Run the SCF and NSCF calculations with QuantumEspresso
* Go in the <code>bn.save</code> folder and run p2y
* Run the setup, just execute yambo
* Generate the input file for the GW calculation: yambo -g n -r -X p -gw0 p -F gw.in. Notice that if you want to use a scissor operator in the non-linear response you can skip this calculation. Notice also that we added the flag <code>-r</code> to introduce a cut-off to the Coulomb interaction for low-dimensional systems.

dyson # [R] Dyson Equation solver
rim_cut # [R] Coulomb potential
ppa # [R][Xp] Plasmon Pole Approximation for the Screened Interaction
gw0 # [R] GW approximation
el_el_corr # [R] Electron-Electron Correlation
HF_and_locXC # [R] Hartree-Fock
em1d # [R][X] Dynamically Screened Interaction
RandQpts= 5000000 # [RIM] Number of random q-points in the BZ
RandGvec= 29 RL # [RIM] Coulomb interaction RS components
CUTGeo= "slab z" # [CUT] Coulomb Cutoff geometry: box/cylinder/sphere/ws/slab X/Y/Z/XY..
% CUTBox
0.000000 | 0.000000 | 0.000000 | # [CUT] [au] Box sides
%
Chimod= "HARTREE" # [X] IP/Hartree/ALDA/LRC/PF/BSfxc
% BndsRnXp
1 | 40 | # [Xp] Polarization function bands
%
NGsBlkXp= 3000 mHa # [Xp] Response block size
% LongDrXp
0.000000 | 1.000000 | 0.000000 | # [Xp] [cc] Electric Field
%
PPAPntXp= 27.21138 eV # [Xp] PPA imaginary energy
% GbndRnge
1 | 40 | # [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|14|3|6|
%
Notice that we calculate the GW corrections only for the bands '''3-6''' we want to use in the non-linear response.

* Run the calculation with the command: <code>yambo -F gw.in</code>
* Optional step: you can calculate the BSE response that can be usefull to interprate the non-linear response with the command: yambo -o b -k sex -r -X p -V qp -y d -F bse.in (do not forget to read the Quasi-Particle database by setting KfnQPdb= "E < SAVE/ndb.QP" in the bse.in file.

== Calculations of the non-linear response ==
* Generate input file to remove symmatries: ypp_nl -y
fixsyms
# [R] Remove symmetries not consistent with an external perturbation
% Efield1
0.000000 | 1.000000 | 0.000000 | # First external Electric Field
%
% Efield2
0.000000 | 0.000000 | 0.000000 | # Additional external Electric Field
%
BField= 0.000000 T # [MAG] Magnetic field modulus
Bpsi= 0.000000 deg # [MAG] Magnetic field psi angle [degree]
Btheta= 0.000000 deg # [MAG] Magnetic field theta angle [degree]
RRmTimeRev # Remove Time Reversal

here we remove all symmetries but you it is possible to remove only the ones not compatible with the external field, 
by putting the field version in the ypp input and removing the time reversal symmetry (<code>RmTimeRev</code>) that is not compatible with the real-time dynamics, 
see for instance here: [[Prerequisites for Real Time propagation with Yambo]]

* Go in the <code>FixSymm</code> folder and run setup again

* Calculate the static dielectric constant and the collisions integrals: yambo_rt -X s -collisions -potential h+sex -r -F collisions.in

em1s # [R][Xs] Statically Screened Interaction
collisions # [R] Collisions
dipoles # [R] Oscillator strenghts (or dipoles)
rim_cut # [R] Coulomb potential
RandQpts=5000000 # [RIM] Number of random q-points in the BZ
RandGvec= 29 RL # [RIM] Coulomb interaction RS components
CUTGeo= "slab z" # [CUT] Coulomb Cutoff geometry: box/cylinder/sphere/ws/slab X/Y/Z/XY..
% CUTBox
0.000000 | 0.000000 | 0.000000 | # [CUT] [au] Box sides
%
Chimod= "HARTREE" # [X] IP/Hartree/ALDA/LRC/PF/BSfxc
% BndsRnXs
 1 | 40 | # [Xs] Polarization function bands
%
NGsBlkXs= 3000 mHa # [Xs] Response block size
% LongDrXs
0.000000 | 1.00000 | 0.00000 # [Xs] [cc] Electric Field
%
% COLLBands
 4 | 5 | # [COLL] Bands for the collisions
%
HXC_Potential= "SEX+HARTREE" # [SC] SC HXC Potential
HARRLvcs= 3000 mHa # [HA] Hartree RL components
EXXRLvcs= 3000 mHa # [XX] Exchange RL components
CORRLvcs= 3000 mHa # [GW] Correlation RL components

For the Collisions it is important to use <code>yambo_rt</code> because it is in single-precision and much faster than <code>yambo_nl</code>. 
Notice that we also use the cutoff on the Coulomb interaction <code>-r</code> because this is a 2D system.

* Run the calculations <code>yambo_rt -F collisions.in</code>.
* Generate input file for the non-linear response: yambo_nl -u n -V qp -F nonlinear.in 
* Run the SHG calculation including collisions, with the command: yambo_nl -nl n -V qp

nloptics # [R] Non-linear spectroscopy
DIP_Threads=0 # [OPENMP/X] Number of threads for dipoles
NL_Threads=0 # [OPENMP/NL] Number of threads for nl-optics
% NLBands
 4 | 5 | # [NL] Bands range
%
NLverbosity= "high" # [NL] Verbosity level (low | high)
NLtime=42.000000 fs # [NL] Simulation Time
NLintegrator= "CRANKNIC" # [NL] Integrator ("EULEREXP/RK2/RK4/RK2EXP/HEUN/INVINT/CRANKNIC")
NLCorrelation= "SEX" # [NL] Correlation ("IPA/HARTREE/TDDFT/LRC/LRW/JGM/SEX")
NLLrcAlpha= 0.000000 # [NL] Long Range Correction
% NLEnRange
 2.000000 | 8.000000 | eV # [NL] Energy range
%
NLEnSteps= 80 # [NL] Energy steps
NLDamping= 0.200000 eV # [NL] Damping (or dephasing)
RADLifeTime=-1.000000 fs # [RT] Radiative life-time (if negative Yambo sets it equal to Phase_LifeTime in NL)
#EvalCurrent # [NL] Evaluate the current
#FrPolPerdic # [DIP] Force periodicity of polarization respect to the external field
HARRLvcs= 3000 mHa # [HA] Hartree RL components
EXXRLvcs= 3000 mHa # [XX] Exchange RL components
GfnQPdb= " E < SAVE/ndb.QP" # [EXTQP G] Database action
.............
% Field1_Freq
0.100000 | 0.100000 | eV # [RT Field1] Frequency
%
Field1_NFreqs= 1 # [RT Field1] Frequency
Field1_Int= 1000.00 kWLm2 # [RT Field1] Intensity
Field1_Width= 0.000000 fs # [RT Field1] Width
Field1_kind= "SOFTSIN" # [RT Field1] Kind(SIN|COS|RES|ANTIRES|GAUSS|DELTA|QSSIN)
Field1_pol= "linear" # [RT Field1] Pol(linear|circular)
% Field1_Dir
 0.000000 | 1.000000 | 0.000000 | # [RT Field1] Versor
%

== Analysis of the non-linear response ==
* Use the command ypp_nl -u to analyze the result
nonlinear # [R] Non-linear response analysis
Xorder= 4 # Max order of the response/exc functions
% TimeRange
40.000 | -1.00000 | fs # Time-window where processing is done
%
ETStpsRt= 200 # Total Energy steps
% EnRngeRt
0.00000 | 20.00000 | eV # Energy range
%
DampMode= "NONE" # Damping type ( NONE | LORENTZIAN | GAUSSIAN )
DampFactor= 0.000000 eV # Damping parameter
PumpPATH= "none" # Path of the simulation with the Pump only

now you can plot columns 4 and 5 of the o.YPP-X_probe_order_2 file that are the real and imaginary part of Xhi2 along yyy direction.
Using gnuplot, just plot

p 'o.YPP-X_probe_order_2' u 1:(sqrt($4**2+$5**2)) t 'SHG in hBN monolayer' w lp

and you will get

[[File:Xhi2 hbn.png|center]]

Alternative to ypp_nl you can also use YamboPy to analyze and plot the data see here: [https://wiki.yambo-code.eu/wiki/index.php?title=Real_time_approach_to_non-linear_response_(SHG)#Analysis_of_the_results_using_YamboPy Analysis of the results using_YamboPy]

Here after we resume with a simple scheme all the steps to calculate SHG at the GoWo+BSE level (courtesy of Sohan, Central University of Punjab, India):

[[File:SHG including GoWo BSE.png|center|SHG GoWo+BSE]]

If you arrived up here, you have the right to listen to the [https://www.youtube.com/watch?v=gz2WMvHCJOU Yambo song].

== Tips and tricks ==

* When you calculate the dielectric constant add the flag -V par to set the parallelization and put as much as possible processors on conduction <code>c</code> and <code>v</code> bands
to decrease the memory
* Regenerate inputs with flag -V RL and decrease the number of total plane-wave to half, this will speed up calculation and save memory, results are the same
* In the non-linear calculation put all processors on the frequencies <code>w</code> this is a super-efficent parallelization but does not distributed the memory

== References ==
<references>
<ref name="Attaccalite2011">C. Attaccalite, M. Grüning, and A. Marini, ''Real-time approach to the optical properties of solids and nanostructures: Time-dependent Bethe-Salpeter equation'', [https://doi.org/10.1103/PhysRevB.84.245110 Phys. Rev. B '''84''', 245110 (2011)]</ref>
<ref name="Attaccalite2013">C. Attaccalite and M. Grüning, Nonlinear optics from an ab initio approach by means of the dynamical Berry phase: Application to second- and third-harmonic generation in semiconductors, [ https://doi.org/10.1103/PhysRevB.88.235113 Phys. Rev. B '''88''', 235113 (2013)]</ref>
</references>

2025-05-26T10:35:33Z

Yambowikiadmin: Page upgraded to contain the full list of tutorials

This page contains the full list of available tutorials. These cover a variety of mixed topics, both physical and methodological. Some of these are om "Standalone mode" and are designed to be followed from start to finish in one page. Others have been ported in "Modular form". These are also available under the link [[Tutorials modular]]. The tutorials do not require previous knowledge of yambo, unless explicitly specified. Each tutorial requires download of a specific core database, and typically they cover a specific physical system (like bulk GaSb or a hydrogen chain). Ground state input files and pseudopotentials are provided. Output files are also provided for reference.

These tutorials can be accessed directly from this page of from the side bar. They include different kind of subjects:

'''Warning''': These tutorials were prepared using previous version of the Yambo code: some command lines, variables, reports and outputs can be slightly different from the last version of the code. Scripts for parsing output cannot work anymore and should be edited to work with the new outputs. New command lines can be accessed typing <code>yambo -h </code>

== Basic ==
* [[Silicon|GW on bulk silicon]]
* [[GW on h-BN (standalone)| GW in practice: how to obtain the quasi-particle band structure of a bulk material]]
* [[LiF|BSE in lithium-floride]]
* [[How to analyse excitons]]

== More on GW and Quasi-particles ==
* [[Real_Axis_and_Lifetimes|Real Axis and Lifetimes]]
* [[Self-consistent GW on eigenvalues only]]
* [[GW tutorial on HPC]]
* [[Quasi-particles and Self-energy within the Multipole Approximation (MPA)]]
* [http://www.attaccalite.com/gw-correction-on-an-arbitrary-point-of-the-brillouin-zone GW on arbitrary k point] (tutorial to be imported into the wiki)
* [[SOC|Spin-orbit coupling in GaSb]]

== More on Linear response and BSE ==
* [[Hydrogen chain|TDDFT Failure and long range correlations]]
* [[SOC|Spin-orbit coupling in GaSb]]
* [[The magneto-optical Kerr effect|The magneto-optical Kerr effect in bulk iron, RPA]]
* [[The magneto-optical Kerr effect including excitonic effects|The magneto-optical Kerr effect in XXX, excitonic effects]] (todo)
* [[Dichroism in molecules]]
* [[Surface spectroscopy]] (to restore from older version)
* [[Si_Surface|Linear Response of a silicon surface (2D)]]
* [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] (to be imported)
* [[Si_wire|Linear Response from a silicon wire (1D)]] (to restore from old version)
* [[H2|H2 molecule: Linear Response & TDDFT (0D)]] (to restore from old version)
* [[SiH4|SiH4: isolated molecule TDDFT & BSE (0D)]] (restored from old version, still to be finalized)

== Post Processing with ypp ==
* [[Yambo Post Processing (ypp)]]

== Electron phonon coupling ==
* [[Electron Phonon Coupling|Electron Phonon Coupling]]
* [[Optical properties at finite temperature]]
* [[Phonon-assisted luminescence by finite atomic displacements]]
* [[Exciton-phonon coupling and luminescence]]

== Real time & Non linear response ==
* [[Linear response using Dynamical Berry Phase]]
* [[Linear response in velocity gauge]]
* [[Linear response from real time simulations]] (modular tutorial)
* [[Nonequilibrium absorption in bulk silicon]] (under construction)
* [[Real time approach to non-linear response (SHG)]]
* [[Correlation effects in the non-linear response]]
* [[SHG within the TD-aGW level (also called TD-HSEX, TD-BSE)]]
* [[A fast approach to excitonic effects in linear/non-linear response]]
* [[Angular dependence of non-linear response]]
* [[Third Harmonic Generation (THG)]]
* [http://www.attaccalite.com/lumen/spin_orbit.html Spin-orbit coupling and non-linear response] (tutorial to be imported into the wiki)
* [[Two-photon absorption]]
* [[Sum frequency generation]]
* [[Shift current]]
* [[Pump and Probe|Pump and Probe: Transient Absorption]]
* [[Pump and Probe: Time resolved ARPES]]
* [[Parallelization for non-linear response calculations]]

== Developing Yambo ==
* [[How to create a new project in Yambo]]
* [[How to create a new ypp interface]]
* [[Some hints on github]]

2025-05-19T11:53:11Z

Yambowikiadmin: Uploaded own work with UploadWizard

Second-harmonic generation of 2D-hBN (5.3)

2025-05-19T11:17:51Z

Yambowikiadmin:

= Step 0: Theoretical framework =

In this tutorial, we compute the Second-harmonic generation (SHG) from the dynamics of the Bloch-states. The equation-of-motion of the Bloch-states is explained in the tutorial on [[Dielectric function from Bloch-states dynamics (5.3)|Dielectric function from Bloch-states dynamics]]. Independently of the 'experiment' we simulate, the part of the integration of motion stays the same. This means that any experiment, including non-linear optics, can be performed at the level of the theory listed for the computation of the dielectric function:
* independent (quasi)particle
* time-dependent Hartree (RPA level)
* time-dependent DFT (and DPFT)
* time-dependent Hartree+Screened exchange (BSE level)

What changes when we want to calculate the SHG (or higher harmonics) is:
* the time-dependence of the input electric field and
* the post-processing of the signal

'''Time dependence of the electric Field''': we choose a sinusoidal electric field, thus a monochromatic laser field. This allows one to expand the polarization in the form <math>\bf{P}(t) = \sum_{n=-\infty}^{+\infty} \bf{p}_n e^{-i\omega_n t} </math> where the coefficient <math>\bf{p}_1,...,\bf{p}_n </math> are related to <math>\chi^{(1)},...,\chi^{(n)} </math> through a coefficient depending on a power of the strength of the field (with the power depending on the order of the response).

At difference with a delta-like perturbation, a real-time simulation gives the response at the laser-field only. Then, to obtain the spectrum for the desired range of frequency, we have to perform so many simulations as the frequencies in the desired range.

'''Post-processing of the signal''': The switch-on of the electric field corresponds to a weak delta-like kick. Thus, even if it may not be noticeable from the polarization plot, the polarization results from both the sinusoidal field and the delta-like kick. The latter excites all eigenfrequencies of the system. Though weak, the resulting signal is comparable with the second-harmonic signal. To eliminate the signal from the eigenfrequencies, we apply a dephasing (<math>\gamma_{deph}</math>). To have a signal useful to sample the second-harmonic signal, we need to wait a time <math>\bar t</math> much larger than the dephasing-time <math>1/\gamma_{deph} </math>. Note that we cannot choose a too short dephasing time (to shorten the simulation time), as this would result in a too large broadening of the spectrum. After <math>\bar t</math>, we sample 2N+1 times in a period and use discrete Fourier transform to extract <math>{p}_n</math> for <math> n = 0,...,N</math>, where n=2 gives the SHG. This is schematically illustrated in figure:

[[File:Pt analysis.png|600px|center|Non-linear response analysis]]

The scheme from Ref. <ref name="Attaccalite2013">C. Attaccalite and M. Gruning [https://arxiv.org/abs/1309.4012v2 Phys. Rev. B, '''88''', 235113 (2013)]</ref> summarised the workflow for computing the SHG (or higher-harmonics) in a energy range <math>[\Omega_1,\Omega_2]</math>.

[[File:Scheme nl.png|350px|center|Schematic representation of real-time calculations]]

= Step 1: Prerequisites =
You should be in the folder <code>YAMBO_TUTORIALS/hBN-2D/YAMBO/</code>.
In this example, we will consider a single layer of hexagonal boron nitride (hBN).
Execute the commands (if running on 4 cores interactively you may need to add -nompi to the command line)
yambo_nl
ypp_nl -fixsym -F 00_removesym.in
to create the input, to remove the symmetries which are not compatible with the field:
fixsyms # [R] Remove symmetries not consistent with an external perturbation
% Efield1
1.000000 | 1.000000 | 0.000000 | # First external Electric Field
%
% Efield2
0.000000 | 0.000000 | 0.000000 | # Additional external Electric Field
%
BField= 0.000000 T # [MAG] Magnetic field modulus
Bpsi= 0.000000 deg # [MAG] Magnetic field psi angle [degree]
Btheta= 0.000000 deg # [MAG] Magnetic field theta angle [degree]
#RmAllSymm # Remove all symmetries
RmTimeRev # Remove Time Reversal
#RmSpaceInv # Remove Spatial Inversion
#KeepKGrid # Do not expand the k-grid
where we chose the field along the xy direction and we remove time-reversal symmetry.

Then run the pre-processing job:
ypp_nl -F 00_removesym.in
This creates the directory <code>FixSymm</code>. Change directory to <code>FixSymm</code>. This contains the <code>SAVE</code> directory with reduced symmetries, compatible with the direction of the field. Now run the setup again (<code>yambo_nl</code>).

= Step 2: Independent-particle approximation =

Use the command:
cd FixSymm
yambo_nl
yambo_nl -nl n -F 01_SHG_ip.in
to generate the input:

nloptics # [R] Non-linear spectroscopy
% NLBands
3 | 6 | # [NL] Bands range
%
NLverbosity= "high" # [NL] Verbosity level (low | high)
NLtime=-1.000000 fs # [NL] Simulation Time
NLintegrator= "INVINT" # [NL] Integrator ("EULEREXP/RK2/RK4/RK2EXP/HEUN/INVINT/CRANKNIC")
NLCorrelation= "IPA" # [NL] Correlation ("IPA/HARTREE/TDDFT/LRC/LRW/JGM/SEX")
NLLrcAlpha= 0.000000 # [NL] Long Range Correction
% NLEnRange
0.5000000 |8.000000 | eV # [NL] Energy range
%
NLEnSteps= 12 # [NL] Energy steps
NLDamping= 0.200000 eV # [NL] Damping (or dephasing)
RADLifeTime=-1.000000 fs # [RT] Radiative life-time (if negative Yambo sets it equal to Phase_LifeTime in NL)
#EvalCurrent # [NL] Evaluate the current
#FrPolPerdic # [DIP] Force periodicity of polarization respect to the external field
% Field1_Freq
0.100000 | 0.100000 | eV # [RT Field1] Frequency
%
Field1_NFreqs= 1 # [RT Field1] Frequency
Field1_Int= 1000.00 kWLm2 # [RT Field1] Intensity
Field1_Width= 0.000000 fs # [RT Field1] Width
Field1_kind= "SIN" # [RT Field1] Kind(SIN|COS|RES|ANTIRES|GAUSS|DELTA|QSSIN)
Field1_pol= "linear" # [RT Field1] Pol(linear|circular)
% Field1_Dir
1.000000 | 1.000000 | 0.000000 | # [RT Field1] Versor
%
Field1_Tstart= 0.010000 fs # [RT Field1] Initial Time

To describe the Bloch-dynamics we will use as a basis the KS-states 3-6 for the k-grid used in the non-scf DFT calculations, that is a 6x6x1 (<code>NLBands</code>). We choose 12 equally spaced frequencies (<code>NLEnSteps</code>) for the applied field in the range (<code>NLEnRange</code>) between 0.5 and 8 eV. We choose a sinusoidal time-dependence (<code>Field1_kind</code>) and set the field direction (<code>Field1_kind</code>) along xy (in plane). The negative simulation time (<code>NLtime=-1.000000</code>) means that the code will choose it based on the value of the <code>NLDamping</code>. While this is a good lower bound for the simulation time, one should check it is sufficient to get accurate results.

Run the simulation, possibly in parallel e. g. in a interactive session on 4 cores.

'''[In order to get 4 cores while in an interacting node, exit your current interactive run and rerun <code>srun</code> asking for the proper number of tasks:'''
srun -A tra25_yambo -p dcgp_usr_prod --reservation=s_tra_yambo -N 1 -n 4 -c 4 -t 04:00:00 --gres=tmpfs:10g --pty /bin/bash
'''Otherwise, you can submit the job with a slurm script as explained in the introductory page]'''
Note that when executing commands that don't need mpi (such as generating an input file) you may need to add -nompi to the command line

Now execute the command:
mpirun -n 4 yambo_nl -F 01_SHG_ip.in -J 01_SHG_ip -C 01_SHG_ip_files

This run will take several minutes. While it is running, you can check the progress of one specific frequency:
tail -f 01_SHG_ip_files/o-01_SHG_ip.polarization_F2
If you want to know how many fs the simulation will last (since this was decided by the code), you can get this information in the report file:
$ grep 'Total simulation' 01_SHG_ip_files/r-01_SHG_ip_nloptics
Total simulation time : 47.81405 [fs]

If you instead want to track the progress of the full simulation, you should check the log file(s). If you are doing a parallel run, you can type
tail -f 01_SHG_ip_files/LOG/l-01_SHG_ip_nloptics_CPU_1

Once the run is finished, check the report <code>01_SHG_ip_files/r-01_SHG_ip_nloptics</code>. In particular, in the appended input files, you can see the default parallelization applied by the code for calculating the dipoles and run the time-propagation. For the latter, the main parallelization is on the frequencies and then on the k-points (<code>NL_ROLEs= "w.k"</code>).

| NL_CPU= "2.2" # [PARALLEL] CPUs for each role
| NL_ROLEs= "w.k" # [PARALLEL] CPUs roles (w,k)
| DIP_CPU= "2.2.1" # [PARALLEL] CPUs for each role
| DIP_ROLEs= "v.c.k" # [PARALLEL] CPUs roles (k,c,v)

Also, look at the simulation time. This is approximately 45-50 fs.

Plot the polarization for a couple of frequencies, e.g. the first and last of the chosen range (the corresponding frequency can be read in the file, search for "Frequency value"):
set ylabel "x-component polarization" offset -0.7
set xlabel "time (fs)"
plot '01_SHG_ip_files/o-01_SHG_ip.polarization_F1' u 1:2 w l lw 2 title "omega = 0.5 eV"
replot '01_SHG_ip_files/o-01_SHG_ip.polarization_F12' u 1:2 w l lw 2 title "omega = 7.375 eV"

[[File:Berry_phase_polarization_of_2D_h-BN_obtained_from_the_previous_run_for_two_laser_frequencies.png|600px|center|Berry_phase_polarization_of_2D_h-BN_obtained from the run for two laser frequencies]]

After the applied field is 'switched on', the polarization is oscillating at the frequency of the applied field and at the eigenfrequencies of the system as described above in the cartoon. This is visible particularly for the higher frequency. Due to the dephasing, after about 20 fs the polarization is oscillating (at least to the eye) at the frequency of the applied field. Nonlinear components are not visible since they are several orders of magnitude smaller than the first order.

== Output post-processing: the dielectric function ==

Once we obtained the polarization, we Fourier-transform the polarization to obtain the dielectric function.

Use the command:

ypp_nl -nl -F 02_SHG_ip_pp.in -C 01_SHG_ip_files

to generate the input file:

nonlinear # [R] Non-linear response analysis
Xorder= 4 # Max order of the response/exc functions
% TimeRange
-1.000000 |-1.000000 | fs # Time-window where processing is done
%
ETStpsRt= 200 # Total Energy steps
% EnRngeRt
0.00000 | 20.00000 | eV # Energy range
%
DampMode= "NONE" # Damping type ( NONE | LORENTZIAN | GAUSSIAN )
DampFactor= 0.000000 eV # Damping parameter
PumpPATH= "none" # Path of the simulation with the Pump only

Where we changed the <code>Xorder</code> to 4. This variable controls the terms in the Fourier expansion of the polarization. We choose up to expand up to the fourth order. This should ensure the second order to be accurate enough. As an exercise, you can change this value (e.g. choose 2,3,5,6) and see how the result below is changing. The time-window where processing is done is chosen automatically, since values are negative. The code chooses by default the last period to carry out the analysis.

To run the post-processing, use the command:
ypp_nl -F 02_SHG_ip_pp.in -J 01_SHG_ip -C 01_SHG_ip_files

In <code>01_SHG_ip_files</code> the following files were generated
o-01_SHG_ip.YPP-X_probe_order_1
o-01_SHG_ip.YPP-X_probe_order_2
o-01_SHG_ip.YPP-X_probe_order_3
o-01_SHG_ip.YPP-X_probe_order_4
r-01_SHG_ip_nonlinear

The output files contain respectively the first, second, third and fourth order response(corresponding to <code>Xorder=4</code>) at the 12 frequencies in the chosen range.

Inspect <code>o-01_SHG_ip.YPP-X_probe_order_2</code>
# E [eV] X/Im[cm/stV](x) X/Re[cm/stV](x) X/Im[cm/stV](y) X/Re[cm/stV](y) X/Im[cm/stV](z) X/Re[cm/stV](z)
#
0.50000000 -0.63057187E-09 -0.19596550E-07 -0.36238688E-10 -0.72052719E-09 0.0000000 0.0000000
1.1250000 -0.19604352E-08 -0.23352876E-07 -0.18347724E-09 -0.86959709E-09 0.0000000 0.0000000
1.7500000 -0.64748829E-08 -0.36274951E-07 -0.27614124E-09 -0.17103902E-08 0.0000000 0.0000000
2.37500000 -0.241985981E-7 -0.864285561E-8 -0.160616725E-8 -0.190064969E-8 0.00000000 0.00000000

The first column are the frequencies of the applied field, the second and third is the polarization along x (imaginary and real part), the fourth and fifth is the polarization along y (imaginary and real part)
and the sixth and seventh is the polarization along z (imaginary and real part).

Plot the absolute value of the SHG for the "xxy" component (that is, the applied field is in the x and y directions and one look at the polarization along x):

set ylabel "{/Symbol=26 \174}{/Symbol=26 c}_{xyx}^{(2)}{/Symbol=26 \174} (pm/V)" offset -0.7
set xlabel "Energy (eV)"
pmVm1toCGS=2.38721e-09
plot '01_SHG_ip_files/o-01_SHG_ip.YPP-X_probe_order_2' u 1:(sqrt($2**2+$3**2)/pmVm1toCGS) w p title "calc"

[[File:02 SHG IP.png|600px|center|SHG_intensity_of_2D_h-BN obtained from the current run (12 frequencies) compared with a run with 112 frequencies]]

In the figure above, the results of the current run are compared with those for a run with 112 frequencies (one can obtain these results by repeating the calculations, changing the number of frequencies in <code>NLEnSteps</code> though this take about 10 times more computational time than the run with 12).

As an exercise, you can repeat the simulations for e,g, a time of <code>35 fs</code> nand <code>60 fs</code>, that is setting <code>NLtime</code> to the desired time. The comparison of the susceptibilities at different simulation times (as in the plot below) shows the importance of running long enough simulations. In too short simulations, the response signal coming from the switch-on of the electric field (delta-kick) is still of comparable intensities as the relevant second-order signal.

[[File:02 SHG IP l.png|600px|center|SHG of 2D-hBN with IPA. Effect of different simulation time]]

As a final note, the SHG depends on the size of the vacuum added in the supercell. Instead, one should consider the surface SHG which obtained by considering an effective thickness for the layer. This is usually chosen to be similar to the layer separation in the bulk counterpart.

= Step 3: Hartree+Screened exchange approximation (BSE level) =

== Hartree+Screened exchange collisions ==
Similarly to the linear part, one needs to generate the Screened exchange collisions. The procedure is exactly the same.

Create the input:
yambo_nl -X s -e -v hsex -F 03_coll_hsex.in
modify the input ad follows
em1s # [R][Xs] Statically Screened Interaction
collisions # [R] Collisions
dipoles # [R] Oscillator strenghts (or dipoles)
Chimod= "HARTREE" # [X] IP/Hartree/ALDA/LRC/PF/BSfxc
% BndsRnXs
1 | 20 | # [Xs] Polarization function bands
%
NGsBlkXs= 1000 mHa # [Xs] Response block size
% LongDrXs
1.000000 | 1.000000 | 0.000000 | # [Xs] [cc] Electric Field
%
% COLLBands
4 | 5 | # [COLL] Bands for the collisions
%
HXC_Potential= "SEX+HARTREE" # [SC] SC HXC Potential
HARRLvcs= 1000 mHa # [HA] Hartree RL components
EXXRLvcs= 1000 mHa # [XX] Exchange RL components
CORRLvcs= 1000 mHa # [GW] Correlation RL components

Calculate the collisions (you may want to run in parallel):
yambo_nl -F 03_coll_hsex.in -J 03_coll -C 04_hsex_files

== Run the simulation ==

Create the input:
yambo_nl -nl n -V qp -F 04_SHG_hsex.in

modify the following parts of the input (by this time all the innput variables should be known from the previous exercises):

% NLBands
4 | 5 | # [NL] Bands range
%

NLintegrator= "CRANKNIC" # [NL] Integrator ("EULEREXP/RK2/RK4/RK2EXP/HEUN/INVINT/CRANKNIC")
NLCorrelation= "SEX" # [NL] Correlation ("IPA/HARTREE/TDDFT/LRC/LRW/JGM/SEX")

% NLEnRange
0.500000 |8.000000 | eV # [NL] Energy range
%
NLEnSteps= 12 # [NL] Energy steps

You can choose 4 Energy steps instead to shorten the calculations (otherwise it takes about 8 minutes on 4 cores)

% GfnQP_E
3.000000 | 1.000000 | 1.000000 | # [EXTQP G] E parameters (c/v) eV|adim|adim
%

Field1_kind= "SIN" # [RT Field1] Kind(SIN|COS|RES|ANTIRES|GAUSS|DELTA|QSSIN)

% Field1_Dir
1.000000 | 1.000000 | 0.000000 | # [RT Field1] Versor
%
HARRLvcs= 1000 mHa # [HA] Hartree RL components
EXXRLvcs= 1000 mHa # [XX] Exchange RL components
CORRLvcs= 1000 mHa # [XX] Exchange RL components

Then run the simulation using e.g. mpi and 4 cores:
mpirun -n 4 yambo_nl -F 04_SHG_hsex.in -J 04_SHG_hsex,03_coll -C 04_SHG_hsex_files

== Postprocessing: obtain the nonlinear response ==

Similarly to what we did for the independent particle case, we create the input:
ypp_nl -nl -F 05_SHG_hsex_pp.in

we change <code>Xorder = 4 </code> and then run:
ypp_nl -F 05_SHG_hsex_pp.in -J 04_SHG_hsex -C 04_hsex_files

Then we plot (as before we compare with a simulation ran for 112 frequencies).

Plot the results:

set ylabel "{/Symbol=26 \174}{/Symbol=26 c}_{xyx}^{(2)}{/Symbol=26 \174} (pm/V)" offset -0.7
set xlabel "Energy (eV)"
pmVm1toCGS=2.38721e-09
plot '04_hsex_files/o-04_SHG_hsex.YPP-X_probe_order_2' u 1:(sqrt($2**2+$3**2)/pmVm1toCGS) w p title "calc"

[[File:02 SHG SEX.png|600px|center|SHG at the Hartree + Screened exchange from the current run with 12 frequencies (crosses) and a run with 112 frequencies (keeping the other parameters the same)]]

We note some differences between the two simulations, this is due to the too short simulation time. Repeating the calculations with a longer time gives more accurate results.
In addition, note that in the interest of time, all parameters are '''under converged'''. In principle, to calculate the collisions one should use the same parameters of a converged BSE calculation. As well, to eliminate spurious interaction with the periodic images, one should use the Coulomb cut-off as it was used for this system in the [[Quasi-particles of a 2D system| tutorial on quasiparticle calculations]]. As a reference, a converged calculation for this system is in <ref name="Gruning2014">M. Grüning and C. Attaccalite, [https://journals.aps.org/prb/abstract/10.1103/PhysRevB.89.081102 Phys. Rev. B '''89''', 081102]</ref>.

= References =
<references/>

=== Links for 2025 school ===
* Back to [[Modena 2025#Tutorials]]

Dielectric function from Bloch-states dynamics (5.3)

2025-05-19T11:17:25Z

Yambowikiadmin:

= Step 0: Theoretical background =

The dynamics of Bloch-states is found by integrating a set of Time Dependent Effective Schrödinger Equations (TD-ESEs). Effective means that an effective Hamiltonian is considered. We consider (see also below) the many-body effective Hamiltonians including electron-electron correlation at mean-field (Hartree) and at the screened-exchange level. In addition, one can also consider effective Hamiltonians based on density-functional theory. The latter choice gives time-dependent density-functional theory.

The coupling between electrons and the external field is described by means of the [http://www.theochem.unito.it/didattica/tec-comp_sdm/note1.pdf Modern Theory of Polarization]:

[[File:Tdse new.png|center|400px|Time-Dependent Schrödinger Equation]]

Where <math> | v_{i \mathbf{k}} \rangle </math> are the time-dependent Bloch states (corresponding to the filled Bloch-states at zero-field), <math>\mathcal E(t)</math> is the external field and the term <math> i e \partial k</math> is the dipole operator consistent with periodic boundary condition. For more details on this last term, see Ref. <ref>I. Souza, J. Íñiguez, and D. Vanderbilt, [https://cfm.ehu.es/ivo/publications/souza_prb04.pdf PRB 69, 085106 (2004)]</ref> and <ref name="Attaccalite2013">C. Attaccalite and M. Gruning [https://arxiv.org/abs/1309.4012v2 Rev. B, '''88''', 235113 (2013)]</ref>.

<math>h_{rt}</math> is the Hamiltonian,
[[File:H mb.png|400px|center|Many-Body Hamiltonian]].
This contains the equilibrium Hamiltonian, that is the zero-field eigenvalues evaluated through DFT; the quasiparticle corrections, evaluated from GW or estimated from experiment, and the variation of the self-energy. The latter is a functional of the density matrix <math>\rho</math>, which is obtained from the Bloch states.

From the solution of the equation of motion for the Bloch states, the polarization is calculated by means of Berry's phase

[[File:Berry polarization.png|center|350px|Berry's polarization]].

One can decide the level of theory by selecting the terms to include in the Hamiltonian:

* by including only the first term, one selects the ''independent-particle approximation''
* by including the first and second terms, one selects the ''independent-quasiparticle approximation''
* by including the first, second terms and the Hartree part of the self-energy, one selects the ''time-dependent Hartree or random-phase approximation''
* by including all terms, one selects the ''time-dependent Screened-exchange or Bethe-Salpeter equation'' level of theory.

The type of "experiment" to perform, or of the response one would like to extract, depend on the time-dependent field in input and the post-processing of the Berry's phase polarization. The induced polarization is the sum of the responses of the system to all orders in the external field:
<math>\bf{P}(t) = \sum_{n=-\infty}^{+\infty} \bf{P}_n (t) </math>

In the case we are interested in linear response only, we use a delta-like electric field that excites all frequencies of the system. We choose the intensity of the field <math>E</math> to be weak enough so that the response is dominated by the first order term. The response at all frequencies is obtained by taking the Fourier-transform of the polarization as

<math> \chi(\omega) = \frac{P(\omega)}{E} </math>.

The response <math> \chi</math> is related to the dielectric function through the relation <math>\epsilon(\omega) = 1 + 4 \pi \chi(\omega) </math>.
See also Ref. <ref name="Attaccalite2013"></ref>.

== Notes: ==
* Length vs Velocity gauge: ''the equations of motion of Yambo are in the length gauge. Other codes choose instead to work in the velocity gauge. If you want to know more about the advantages and disadvantages of the two gauges read section 2.7 of Ref. <ref>C. Attaccalite, [https://arxiv.org/abs/1609.09639 arXiv 1609.09639 (2017)]</ref>.''

* Advantages and disadvantages with respect to the density-matrix based formalism: ''Yambo implements two distinct formalisms for studying the real-time response, the one based on density-matrix and the one based on Boch states. Which formalism is best to use depends on what one wants to simulate. The dynamics formulation in terms of density matrix or non-equilibrium Green's functions allows a systematic treatment of correlation effects and electron-phonon coupling<ref name="DavideEPL">[https://arxiv.org/abs/1409.1706 D. Sangalli, and A. Marini EPL '''110''', 47004 (2015)]</ref>, but the price to pay is that the coupling with the external field is correct only to the linear order. This formalism is important in all those phenomena where electronic relaxation plays a major role. On the other hand, the formulation in terms of the Bloch states treats coupling with the external field in an exact way even beyond the linear regime, and for this reason, it allows the description of phenomena coherent with the external fields and to access, for example, non-linear responses.''

= Step 1: Prerequisites =
In this example, we will consider a single layer of hexagonal boron nitride (hBN). Before running real time simulations in yambo you need to:
* download input files and Yambo databases for this tutorial here: [https://media.yambo-code.eu/educational/tutorials/files/hBN-2D-RT.tar.gz hBN-2D-RT.tar.gz].
* perform the steps described in [[Prerequisites for Real Time propagation with Yambo (5.3)]].

= Step 2: Independent-particle approximation =

Let's start from the lowest level of theory. Here it is assumed you are familiar with the optical spectrum of 2D h-BN at this level of theory.

== Run the simulations ==
You should now be inside the folder <code>YAMBO_TUTORIALS/hBN-2D-RT/YAMBO/FixSymm</code>.
First of all, create a directory for the inputs to be run by <code>yambo_nl</code>:
mkdir Inputs_nl

Then, use the command:
yambo_nl -nl n -F Inputs_nl/01_td_ip.in

to generate the input (OMP, angular and extra fields keywords are omitted for brevity):

nloptics # [R] Non-linear spectroscopy
% NLBands
3 | 6 | # [NL] Bands range
%
NLverbosity= "high" # [NL] Verbosity level (low | high)
NLtime= 55.000000 fs # [NL] Simulation Time
NLintegrator= "INVINT" # [NL] Integrator ("EULEREXP/RK2/RK4/RK2EXP/HEUN/INVINT/CRANKNIC")
NLCorrelation= "IPA" # [NL] Correlation ("IPA/HARTREE/TDDFT/LRC/LRW/JGM/SEX")
NLLrcAlpha= 0.000000 # [NL] Long Range Correction
% NLEnRange
-1.000000 |-1.000000 | eV # [NL] Energy range
%
NLEnSteps= 1 # [NL] Energy steps
NLDamping= 0.000000 eV # [NL] Damping (or dephasing)
RADLifeTime=-1.000000 fs # [RT] Radiative life-time (if negative Yambo sets it equal to Phase_LifeTime in NL)
#EvalCurrent # [NL] Evaluate the current
#FrPolPerdic # [DIP] Force periodicity of polarization respect to the external field
HARRLvcs= 18475 RL # [HA] Hartree RL components
EXXRLvcs= 18475 RL # [XX] Exchange RL components
% Field1_Freq
0.100000 | 0.100000 | eV # [RT Field1] Frequency
%
Field1_NFreqs= 1 # [RT Field1] Frequency
Field1_Int= 1000.00 kWLm2 # [RT Field1] Intensity
Field1_Width= 0.000000 fs # [RT Field1] Width
Field1_kind= "DELTA" # [RT Field1] Kind(SIN|COS|RES|ANTIRES|GAUSS|DELTA|QSSIN)
Field1_pol= "linear" # [RT Field1] Pol(linear|circular)
% Field1_Dir
0.000000 | 1.000000 | 0.000000 | # [RT Field1] Versor
%
Field1_Tstart= 0.010000 fs # [RT Field1] Initial Time

Note that the standard input of <code>yambo_nl</code> has default settings for nonlinear response (e.g. the <code>SOFTSIN</code> for the <code>Field1_kind</code> that set the time-dependence 'kind' of the field).
Set the field direction <code>Field1_Dir</code> along y, the field kind to <code>DELTA</code>, the length of the simulation <code>NLtime</code> to 55 fs, select the bands from 3 to 6 (<code>NLBands</code>) and set the <code>NLDamping</code> to zero. The fields you need to change with respect to the default settings are shown above in red.

Now run
yambo_nl -F Inputs_nl/01_td_ip.in -J TD-IP_nl -C TD-IP_nl &
and then monitor the progress with
tail -f TD-IP_nl/o-TD-IP_nl.polarization_F1
ctrl+C

The code produces different files (in the <code>TD-IP_nl</code> directory, all with the suffix <code>TD-IP_nl</code>: <code>o-TD-IP_nl.polarization_F1</code> that contains the polarization, <code>o-TD-IP_nl.external_potential_F1</code> the external field we used, and finally <code>r-TD-IP_nl_nloptics</code> a report with all information about the simulation.

Plot the third column of <code>o-TD-IP_nl.polarization_F1</code> versus the first one (time-variable) to get the time-dependent polarization along the y-direction:
gnuplot> p 'TD-IP_nl/o-TD-IP_nl.polarization_F1' u 1:3 w l

[[File:01 TDpolarization delta NL.png|center|600px|TD polarization from a delta-kick in the IPA in 2D-hBN]]

From the plot, it is apparent that the polarization oscillates at multiple frequencies.
== Output postprocessing: the dielectric function ==

Once we obtained the polarization, we Fourier transform the polarization to obtain the dielectric function.

Use the command:

ypp_nl -nl -F Inputs_nl/02_ypp_abs.in

to generate the input file:

nonlinear # [R] NonLinear Optics Post-Processing
Xorder= 1 # Max order of the response functions
% TimeRange
-1.000000 |-1.000000 | fs # Time-window where processing is done
%
ETStpsRt= 1001 # Total Energy steps
% EnRngeRt
0.00000 | 10.00000 | eV # Energy range
%
DampMode= "LORENTZIAN" # Damping type ( NONE | LORENTZIAN | GAUSSIAN )
DampFactor= 0.10000 eV # Damping parameter

where we set a Lorentzian smearing corresponding to 0.1 eV. Note that due to the finite time of our simulation, we need to introduce a finite smearing to Fourier transform the result.

To run the post-processing, use the command:
ypp_nl -F Inputs_nl/02_ypp_abs.in -J TD-IP_nl -C TD-IP_nl
and obtain the files for the dielectric constant along with the field direction, the EELS along with the same direction, and the damped polarization.
o-TD-IP_nl.YPP-eps_along_E
o-TD-IP_nl.YPP-eels_along_E
o-TD-IP_nl.YPP-damped_polarization
inside the folder <code>TD-IP_nl</code>

Plot the dielectric constant and compare it with the linear response and the density-matrix formalism (available from the [[Real time approach to linear response]] tutorial):
gnuplot
plot "../CHI_IP/o-CHI_IP.eps_q1_ip" u 1:2 w l
rep "TD-IP_nl/o-TD-IP_nl.YPP-eps_along_E" u 1:2 w l

[[File:Bn optics.png|thumb|900px|center|Imaginary part of the dielectric constant]]

The differences observed with the spectra obtained from the linear-response formalisms are due to how the dipoles are evaluated. Within the linear-response formalism, dipoles are computed as <math>\frac{1}{E_{c \mathbf{k}}-E_{v \mathbf{k}}} \left\langle v \mathbf{k}\left|\mathbf{p}_{\alpha}+\mathrm i\left[V^{\mathrm{NL}}, \mathbf{r}_{\alpha}\right]\right| c \mathbf{k}\right\rangle\</math> (which is valid in the linear response regime). The Bloch-states formalism uses instead covariant dipoles evaluated numerically on the k-mesh (valid at all orders). The differences are due to numerical errors in evaluating the covariant dipoles. These differences disappear with a denser k-mesh.

'''It is a good practice to converge the dielectric function obtained from the Bloch-state dynamics with respect to the k-points using the linear-response spectrum as a reference. This gives the minimum number of k-points to be used to obtain accurate results.
'''

= Step 3: Hartree+Screened exchange approximation =

== Evaluating the collisions==

[[File:Self-Energy Rho Linear Expansion.png|center|200px|Expression for the many body self-energy doing a linear expansion in terms of the density matrix]].
K is the functional derivative of the self--energy with respect to the density matrix and corresponds to the kernel of the Bethe--Salpeter equation. In the language of the yambo code we will refer to it with the name "real time collisions" or simply "collisions". Choosing the HSEX approximation and using the GW energies for the quasi--particle corrections, the time propagation corresponds to a real-time version of the Bethe-Salpeter Equation<ref>[http://bcsbec.df.unicam.it/files/LaRNC11_N12%281988%29.pdf G. Strinati, La Rivista del Nuovo Cimento '''11''' (12), 1-86 (1988)]</ref>
More details can be found in Ref. <ref name="Attaccalite2011">[https://arxiv.org/abs/1109.2424 C. Attaccalite, M. Grüning, and A. Marini PRB '''84''', 245110 (2011)]</ref>
(This part is common in between the <code>yambo_nl</code> and the <code>yambo_rt</code>).

Generate the input file to calculate the collisions use the command :
yambo_nl -X s -e -v hsex -F Inputs_nl/03_coll_hsex.in

The flag -X will tell the code to calculate the dielectric constant that is required for the screened interaction.
It could be even computed in an independent calculation and then loaded using the <code>-J</code> flags

em1s # [R Xs] Static Inverse Dielectric Matrix
dipoles # [R ] Compute the dipoles
Chimod= "HARTREE" # [X] IP/Hartree/ALDA/LRC/PF/BSfxc
% BndsRnXs
1 | 20 | # [Xs] Polarization function bands
%
NGsBlkXs= 1000 mHa # [Xs] Response block size
% DmRngeXs
0.10000 | 0.10000 | eV # [Xs] Damping range
%
% LongDrXs
1.000000 | 0.000000 | 0.000000 | # [Xs] [cc] Electric Field
%

While this is the part of the input file specific for the evaluation of the collisions
collisions # [R] Eval the extended Collisions
% [[Variables#COLLBands|COLLBands]]
4 | 5 | # [COLL] Bands for the collisions
%
HXC_Potential= "SEX+HARTREE" # [SC] SC HXC Potential
[[Variables#HARRLvcs|HARRLvcs]]= 1000 mHa # [HA] Hartree RL components
[[Variables#EXXRLvcs|EXXRLvcs]]= 1000 mHa # [XX] Exchange RL components
[[Variables#CORRLvcs|CORRLvcs]]= 1000 mHa # [GW] Correlation RL components

With this input, we calculate the HARTREE plus SEX collisions integrals. Notice that the HARTREE term in principle can be calculated on the fly, but in this way it is more efficient especially for the non-linear response.
Here one has to converge the cutoff for the Hartree and the Screened Exchange. Around <code>5000 mHa</code> is a reasonable value for hBN. In this example we will use <code>1000 mHa</code> to speed up calculations. The collisions bands <code>COLLBands</code> have to be the same number of bands you want to use in the linear/nonlinear response.
Then you run
yambo_nl -F Inputs_nl/03_coll_hsex.in -J COLL_HSEX -C COLL_HSEX

The calculation may take sone time (about 1 minute in serial).
It produced many binary files <code>ndb.COLLISIONS_HXC_fragment_*</code> which will be needed to perform TD-HSEX simulations.

== Run the simulations ==
To generate the input for the real-time simulation you can run

yambo_nl -nl n -V qp -F Inputs_nl/04_td-sex.in

The input file is

nloptics # [R NL] Non-linear optics
% NLBands
 4 | 5 | # [NL] Bands
%
NLverbosity= "high" # [NL] Verbosity level (low | high)
NLtime= 20.00000 fs # [NL] Simulation Time
NLintegrator= "CRANKNIC" # [NL] Integrator ("EULEREXP/RK2/RK4/RK2EXP/HEUN/INVINT/CRANKNIC")
NLCorrelation= "SEX" # [NL] Correlation ("IPA/HARTREE/TDDFT/LRC/LRW/JGM/SEX")
NLLrcAlpha= 0.000000 # [NL] Long Range Correction
% NLEnRange
-1.000000 | -1.000000 | eV # [NL] Energy range
%
NLEnSteps= 1 # [NL] Energy steps
NLDamping= "0.10000 eV # [NL] Damping
#EvalCurrent # [NL] Evaluate the current
#FrPolPerdic # [DIP] Force periodicity of polarization respect to the external field
HARRLvcs= 1000 mHa # [HA] Hartree RL components
EXXRLvcs= 1000 mHa # [XX] Exchange RL components

The next section of the input is about external quasiparticle corrections (<code>EXTQP G</code>). There change only

% GfnQP_E
3.000000 | 1.000000 | 1.000000 | # [EXTQP G] E parameters (c/v) eV|adim|adim
%

The latter line introduces a scissor operator (a rigid shift of the conduction bands) of 3.0 eV. In principle, it is possible to perform a G0W0 calculation with Yambo and use the Quasi-particle band structure instead of the rigid shift.

The last section regards the applied field (<code>RT Field1</code>). There change these two lines:

% Field1_Dir
0.000000 | 1.000000 | 0.000000 | # [RT Field1] Versor
%
Field1_kind= "DELTA" # [RT Field1] Kind(SIN|SOFTSIN|RES|ANTIRES|GAUSS|DELTA|QSSIN)

Run the calculation:

yambo_nl -F Inputs_nl/04_td-sex.in -J "TD-SEX_nl,COLL_HSEX" -C TD-SEX_nl

== Output postprocessing: dielectric function ==

The post-processing of the output does not depend on the level of approximation. The same steps are taken as in the independent-particle approximation.

Use the command:
ypp_nl -F Inputs_nl/ypp_abs.in -J TD-SEX_nl -C TD-SEX_nl

Plot the results:
gnuplot> set xrange [3:10]
gnuplot> set yrange [0:12]
gnuplot> plot "TD-SEX_nl/o-TD-SEX_nl.YPP-eps_along_E" u 1:2 w l

[[File:HBN HSEX absorption.png|thumb|900px|center|In plane absorption of hBN at the TD-HSEX level compared with IP absorption]]

''For comparison, the linear response results at the BSE level can be obtained following the [[How to obtain an optical spectrum|BSE tutorial]]. ''

= References =
<references/>

=== Links for 2025 school ===
* Back to [[Modena 2025#Tutorials]]
* Next step [[Second-harmonic generation of 2D-hBN (5.3)|Second-harmonic generation of 2D-hBN]]

2025-05-19T11:04:24Z

Yambowikiadmin: Uploaded own work with UploadWizard