Lumen wiki - User contributions [en]

Linear response in velocity gauge

2026-07-22T21:46:29Z

Davide: References improved

[[File:Lresp.png|thumb|200px]]
In this tutorial we show how to get the linear response using the velocity gauge instead of the Berry phase approach in length gauge. 
We use the same DFT input files of the tutorial [[Linear response using Dynamical Berry Phase]], run the same setup and remove the symmetries.
 
== Real-time dynamics in velocity gauge==
Then in order to calculate the response in velocity gauge we use the command <code>yambo_nl -u p -V resp -F input_lr.in</code> to generate the input:

nlinear # [R NL] Non-linear optics
NL_Threads= 1 # [OPENMP/NL] Number of threads for nl-optics
% NLBands
3 | 6 | # [NL] Bands
%
NLstep= 0.0100 fs # [NL] Real Time step length
NLtime=55.000000 fs # [NL] Simulation Time
NLintegrator= "INVINT" # [NL] Integrator ("EULEREXP/RK4/RK2EXP/HEUN/INVINT/CRANKNIC")
NLCorrelation= "IPA" # [NL] Correlation ("IPA/HARTREE/TDDFT/LRC/JGM/SEX/HF")
NLLrcAlpha= 0.000000 # [NL] Long Range Correction
NLDamping= 0.000000 eV # [NL] Damping
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
Gauge= "velocity" # [BSE/X] Gauge (length|velocity)
Field1_Freq= 0.100000 eV # [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

Notice that in velocity gauge the current is correct while the polarization is calculated using only intra-band dipoles, that anyway is enough for the linear response.

==Analysis of the results using YamboPy using polarization ==
'''This part works only with Yambo 5.3 or the last version available on Github:''' [https://github.com/yambo-code/yambo https://github.com/yambo-code/yambo] 
You can get the dielectric constant along the field direction using the same procedure of the Berry Phase case, Fourier transforming the polarization:

import numpy as np
from yambopy import *
from yambopy.plot import *

NLDB=YamboNLDB()

pol =NLDB.Polarization[0]
time=NLDB.IO_TIME_points
t_initial=NLDB.Efield[0]["initial_time"]

pol_damped=np.empty_like(pol)
for i_d in range(3):
pol_damped[i_d,:]=damp_it(pol[i_d,:],time,t_initial,damp_type='LORENTZIAN',damp_factor=0.1/ha2ev)

Plot_Pol_or_Curr(time=time, pol=pol_damped, xlim=[0,55], save_file='polarization.pdf')

[[File:Polarization velocity.png|center|Polarization_velocity]]

and then calculate the dielectric constant as in the length gauge

Linear_Response(time=time,pol=-pol_damped,efield=NLDB.Efield[0],plot=False,plot_file='eps.pdf')

notice that in this case, only the linear response along the field is correct.

[[File:Eps velocity.png|center|Eps_velocity]]

==Analysis of the results using YamboPy using current ==

The same analysis can be performed using the current. In this case we remove the constant term then we reconstruct the polarization
and then the epsilon.

import numpy as np
from yambopy import *
from yambopy.plot import *
from scipy.interpolate import CubicSpline
import matplotlib.pyplot as plt

NLDB=YamboNLDB()

current =NLDB.Current[0]
time=NLDB.IO_TIME_points
t_initial=NLDB.Efield[0]["initial_time"]
#
# Remove the constant term
#
ave_current=np.zeros(3,dtype=np.double)
for i_d in range(3):
ave_current[i_d]=np.sum(current[i_d,:])/current.shape[1]
current[i_d,:]=current[i_d,:]-ave_current[i_d]

pol=np.empty_like(current)
for i_d in range(3):
spline_j=CubicSpline(time,current[i_d,:])
for idx, t in enumerate(time):
pol[i_d,idx]=spline_j.integrate(0.0,t)

#
# There is a wrong sign convention in the velocity gauge
# somewhere in the code
#
pol=-pol

Plot_Pol_or_Curr(time=time, pol=pol, xlim=[0,55], save_file='polarization.pdf')
#
pol_damped=np.empty_like(pol)
for i_d in range(3):
pol_damped[i_d,:]=damp_it(pol[i_d,:],time,t_initial,damp_type='LORENTZIAN',damp_factor=0.1/ha2ev)
#
#
# Calculate the linear response
#
Linear_Response(time=time,pol=pol_damped,efield=NLDB.Efield[0],plot=False,plot_file='eps.pdf' )

and you the get the spectra:

[[File:Eps vel.png|center|Yambo tutorial image]]

This procedure produce a small bump at <math>\omega=0</math> that is an artifact, different procedure has been discussed in the literature to remove it (see Refs. <ref name="Zhang2017"/>, <ref name="Sangalli2017"/>)

== References==
<references>
<ref name="Zhang2017">Manipulation of the dielectric properties of diamond by an ultrashort laser pulse, X. Zhang, F. Wang, L. Jiang, and Y. Yao Phys. Rev. '''B''' 95, 184301 (2017)</ref>
<ref name="Sangalli2017">Optical properties of periodic systems within the current-current response framework: Pitfalls and remedies,
D. Sangalli, J.A. Berger, C. Attaccalite, M. Grüning, P. Romaniello,
Phys. Rev. B '''95''', 155203 (2017)</ref>
</references>

Real time approach to non-linear response (SHG)

2026-07-22T21:12:04Z

Davide: removed useless footer

== Second Harmonic Generation in AlAs ==
In this tutorial, we will calculate the second harmonic generation of bulk AlAs, the Lumen databases can be downloaded here: [https://media.yambo-code.eu/educational/tutorials/files/AlAs_DBs.tar.gz AlAs_DBs.tar.gz] (10 MB). The first steps of this tutorials are the same as the one on [http://www.yambo-code.eu/wiki/index.php?title=Prerequisites_for_Real_Time_propagation_with_Yambo Prerequisites for Real-Time propagation with Yambo] with the only difference that we will work on AlAs and we will consider an external field in the direction <code> 1.000000 | 1.000000 | 0.000000 |</code>. The DFT input are available here: [https://media.yambo-code.eu/educational/tutorials/files/AlAs_abinit.tar.gz ABINIT] or [https://media.yambo-code.eu/educational/tutorials/files/AlAs_pwscf.tar.gz QuantumEspresso].
You can run the DFT calculation with ABINIT with the command:

abinit < AlAs.in > output_AlAs

or using QuantumEspresso:

pw.x -inp AlAs.scf.in > output_scf
pw.x -inp AlAs.nscf.in > output_scf

and then import the ABINIT wave-function with the command:

a2y -F AlAso_DS2_WFK.nc

and the QuantumEspresso one, in the folder AlAs.save with the command:

p2y
Now we consider an external field in the <code>[1,1,0]</code> direction and remove symmetries not compatible with this field, as explained in the tutorial [http://www.yambo-code.eu/wiki/index.php?title=Prerequisites_for_Real_Time_propagation_with_Yambo Prerequisites for Real-Time propagation with Yambo].

==Real-time simulation for the SHG==

You can generate the input file with the command <code>yambo_nl -u n</code>:

nloptics # [R NL] Non-linear optics
% NLBands
 3 | 6 | # [NL] Bands
%
NLstep= 0.0100 fs # [NL] Real Time step length
NLverbosity= "high" # [NL] Verbosity level (low | high)
NLtime=-1.000000 fs # [NL] Simulation Time
NLintegrator= "INVINT" # [NL] Integrator ("EULEREXP/RK4/RK2EXP/HEUN/INVINT/CRANKNIC")
NLCorrelation= "IPA" # [NL] Correlation ("IPA/HARTREE/TDDFT/LRC/JGM/SEX/HF")
NLLrcAlpha= 0.000000 # [NL] Long Range Correction
% NLEnRange
 1.000000 | 5.000000 | eV # [NL] Energy range
%
NLEnSteps= 10 # [NL] Energy steps
NLDamping= 0.150000 eV # [NL] Damping
% Field1_Dir
 1.000000 | 1.000000 | 0.000000 | # [RT Field1] Versor
%
Field1_kind= "SOFTSIN" # [NL ExtF] Kind(SIN|SOFTSIN|RES|ANTIRES|GAUSS|DELTA|QSSIN)

The line <code>|1.000000 | 1.000000 | 0.000000 | # [RT Field1] Versor</code> referees to the direction of the external field (x,y,0). 
The default parameters of Lumen are already tuned for second-harmonic generation, so the only thing you have to change is the band range, between 3 and 6 and the energy range between 1.0-5.0 eV and the number of energy steps in this interval that we set to 10. Notice that you cannot set to zero the lowest value of the energy range because this will requires a simulation that lasts infinite time, see below. Finally, consider that Yambo performs a separate calculation for each frequency, so if you set many energy steps the computational time grows linearly with this number. 

Notice that we set <code>NLverbosity= "high"</code> in this way the code will produce a file for each laser frequency containing the time dependent polarization. 

Run <code>yambo_nl</code>. Calculations will take about fifteen minutes on a single processor PC, you have time to study the next section that explains how non-linear response is extracted from the real-time simulations.
In order to speed up calculations, you can run them in parallel.

==Non-linear response with Yambo==

In order to calculate the non-linear response, the system is excited with different laser fields with a sinusoidal shape at frequencies <math>\omega_1, \omega_2, .... , \omega_n </math> determined by the parameters [[Variables#NLEnRangev|NLEnRange]] and [[Variables#NLEnSteps|NLEnSteps]]. A dephasing term is added to the Hamiltonian <math>\gamma = </math>[[Variables#NLDamping|NLDamping]] to simulate a finite broadening and to remove the eigenmodes that are excited by the sudden turning on of the external field. After the dephasing time the outgoing signal is analyzed to extract the non-linear coefficients as shown in the figure below:

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

In the report file <code>r_nlinear</code> you will find the length of the dephasing part and of the sampling one:

Dephasing Time [fs]: 52.65695
Sampling Time [fs]: 4.18566
Total simulation time [fs]: 56.84262
The length of the dephasing interval is inversely proportional to the damping term <math>T_{depth} \simeq 1/\gamma</math> while the length of the sampling is dictated by the smallest frequency we are interested in: <math>T_{samp} \simeq 1/\omega_1 </math>. For this reason if <math> \omega_1=0 </math> or <math> \gamma = 0 </math> simulation time goes to infinity. The response at zero frequency can be calculated as limit of small frequency perturbation.
You can have a look to the file <code>src/nloptics/NL_initialize.F</code> to see how simulation lengths are defined.
Calculations can take some time, run them in parallel, the best number of processors = number of frequencies or reduce the number of frequencies step in the SHG.

==Analysis of the results==

In the sampling region we suppose that the polarization can be written as: <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>. We sample the polarization signal at different times and invert the previous equation by truncating the sum at a finite order <ref>[https://arxiv.org/abs/1309.4012 Nonlinear optics from an ab initio approach by means of the dynamical Berry phase: Application to second- and third-harmonic generation in semiconductors], C. Attaccalite and M. Grüning, Phys. Rev. B 88, 235113(2013)</ref>.
This is done with the command <code>ypp_nl -u</code> that automatically produce an input with the correct values:

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

where <code>Xorder</code> is the order where the previous sum is truncated, and the TimeRange specifies the sampling region. Notice that differently from the first tutorial, in this case, we do not need Damping in <code>ypp_nl</code> because we already included it in the real-time dynamics. Run <code>ypp_nl</code> and it will produce a file called <code>o.YPP-X_probe_order_2</code>. This file contains the different components of <math>\chi^2_{xy*} </math> and in particular the columns 6 and 7 correspond to the imaginary and real part of <math>\chi^2_{xyz} </math>. You can plot the result with gnuplot and the command:

p 'o.YPP-X_probe_order_2' u 1:(sqrt($6**2+$7**2)) w lp ps 1.5 pt 7 lw 1.5

the result should look like

[[File:AlAs results.png|700px|center|SHG in AlAs]]

and compare it with the <math> \chi^2_{zxy} </math> calculated with more frequencies and with the converged result on a <math> 18x18x18</math> k-point grid and bands between 2 and 10. In the figure you can find also the comparison with the results of Ref. <ref>Luppi et al. [https://arxiv.org/abs/1006.2649 PRB B, '''82''', 235201(2010)]</ref>.
Notice that the QuantumEspresso results are slightly different from the Abinit ones. This is due to the different pseudo-potential employed in the calculations if pseudopotentials were the same calculation would have been identical. You can download the script to generate this plot and the converged results [http://www.attaccalite.com/lumen/tutorials/AlAs/AlAs_results.tgz here].

'''Notice that''' : the following parameters are not used in the non-linear response analysis: <code> EnRngeRt</code>, <code>ETStpsRt</code>, <code>DampMode</code>, <code>DampFactor</code>. The damping factor is set directly in the real-time simulation.

'''System of units in Non-linear Optics''': notice that in Yambo the non-linear response functions are in the gaussian system of units, for example the <math> \chi^2 (\omega) </math> is in cm/statvolt. In order to convert the non-linear coefficients between different system of units, you can have a look at the [https://www.attaccalite.com/tutorials_yambo/Appendix_C_Boyd.pdf Appendix C] of the book [https://www.elsevier.com/books/nonlinear-optics/boyd/978-0-12-369470-6 "Non-linear optics" by Robert W. Boyd].

More tutorials on non-linear response (Third harmonic generation, parallelization, spin-orbit etc..) with Yambo can be found here: [http://www.attaccalite.com/lumen/tutorials.html More Tutorials]. 

==Analysis of the results using YamboPy==
The analysis of the result performed in the previous section using <code>ypp_nl</code> can be performed using python script [https://github.com/yambo-code/yambopy YamboPy]. 
In YamboPy we implemented all necessary functions to reads the Yambo databases and post-process the results. Here we show how to get the SHG signal from the previous
simulations. We suppose you correctly installed <code>YamboPy</code> on your PC. Go in the folder where you ran non-linear calculation and type in the python:
 
from yambopy import *
NLDB=YamboNLDB()
pol =NLDB.Polarization[0]
time=NLDB.IO_TIME_points

in this way you read all the Non-linear databases. If your runs are in a different folder then the 'SAVE' one you can specify it using the command:
NLDB=YamboNLDB(calc='MYJOB')

now in the array <code>pol</code> you have all the polarization for all laser frequencies in the three Cartesian directions, while the variable <code>time</code> contains all the time series of your simulation.
Now you can get the non-linear response with the command:

SIN = Xn_from_sine(NLDB,X_order=5)
print(SIN)
OUT = SIN.perform_analysis()
SIN.output_analysis(OUT)

this command will perform a Fourier analysis of the results in the same way of <code>ypp_nl</code> and generate new files with all the requested harmonics ''o.YamboPy-X_probe_order_1'', ''o.YamboPy-X_probe_order_2'' etc...
equivalent to the one generated by <code>ypp_nl</code>. Results can be plot with gnuplot using the command:

p 'o.YamboPy-X_probe_order_2' u 1:(sqrt($6**2+$7**2)) w lp ps 1.5 pt 7 lw 1.5

The script that performs non-linear analysis can be found in <code>yambopy/nl/harmonic_analysis.py</code>. We strongly advice you to have a look to this script and modify it according
to your needs and in order to extract other non-linear response functions.

==How to choose the external field direction for the different X2==

In the X2abc response you have three field directions. The two components ''b'' and ''c'' are determined by the external field and the third ''a''
is the direction where you measure the response. 

The external field in Yambo is set using the variable
% Field1_Dir
1.000000 | 0.000000 | 0.000000 | # [NL ExtF] Versor .
%

For example if you set in Yambo field direction as:

1.0 | 0.0 | 0.0 | ---> you calculate all components X2*xx
0.0 | 1.0 | 0.0 | ---> you calculate all components X2*yy
0.0 | 0.0 | 1.0 | ---> you calculate all components X2*zz

the * correspond to the different direction in the response file <code>o.YPP-X_probe_order_2</code>.
In the <code>o.YPP-X_probe_order_2</code> you have seven columns:

The first column is the energy [[omega]],
the 2nd and 3rd columns are the imaginary and real part of X2x**
the 4th and 5th columns are the imaginary and real part of X2y**
the 6th and 7th columns are the imaginary and real part of X2z**

where the ** are the direction of the incoming field, see above.

For example if you want to calculate the <math>\chi^{(2)}_{zxx}</math> you put the external field in the direction <code> 1.0 | 0.0 | 0.0 |</code>
and then plot the column 6 and 7 of the <code>o.YPP-X_probe_order_2</code> file. 

In order to get the other components of the X2 the calculation are more involved. For example if you put field in the

1.0 | 1.0 | 0.0 |

this corresponds to <math>E(t)=\frac{1}{\sqrt{2}} [ E_x(t) + E_y(t) ]</math> and so at the second order you will have <math>E^2(t)=\frac{1}{2} \left [ E^2_x(t) + E^2_y(t) + E_x(t) E_y(t) + E_y(t) E_x(t) \right]</math> 
this means that the X2 in the <code>o.YPP-X_probe_order_2</code> will correspond to a linear combination of the following different response functions. 
For example in the the 2nd and 3rd columns are the imaginary and real part of: 

<math>\chi^{(2)}_{x**}=\frac{1}{2} \left [\chi^2_{xxx}+\chi^{(2)}_{xyy} + \chi^{(2)}_{xxy} + \chi^{(2)}_{xyx} \right ]</math>

the same for the other column replacing the first <code>x</code> for <code>y,z</code>.
If you know that some components are zero this allow you to extract the corresponding off-diagonal part of <math>\chi^{(2)}_{abc}</math> otherwise you have to perform different calculations to get the component you are interested in. 
Now we present '''a simple example of monolayer hBN (m-hBN)'''. In m-hBN there is only a non-zero component of the <math>\chi^{(2)}</math>and all the others are related by symmetry: 
<math>\chi^{(2)}_{yyy}=-\chi^{(2)}_{yxx}=-\chi^{(2)}_{xxy}=-\chi^{(2)}_{xyx}</math> 
so if you put the field in the <code>|1|1|0|</code> direction you will get: 
<math>\chi^{(2)}_{x**} = \frac{1}{2} \left[\chi^{(2)}_{xxy} + \chi^{(2)}_{xyx} \right] = - \chi^{(2)}_{yyy} </math> 
<math>\chi^{(2)}_{y**} =\frac{1}{2} \left[ \chi^{(2)}_{yyy} + \chi^{(2)}_{yxx} \right] = 0.0 </math>
 
In the figure below you can see the different <math>\chi^{(2)}</math> of monolayer hBN.

[[File:Xhi2 xy.png|center| 800px|xhi2_xy]]

as you can see in the second panel, when the field is in the <code>xy</code> you get a linear combination of <math>\chi^{(2)}</math> corresponding to the equation above, in this case the <code>y</code> is zero while the <code>x</code> component is equal to <math>-\chi^{(2)}_{yyy}</math>.

'''Nota bene:''' if you change the direction of the external field you must to remove the corresponding symmetries as explained in the tutorial [[Prerequisites for Real Time_propagation with Yambo#Reduce symmetries real-time linear response tutorial]]. If you have to calculate <math>\chi^{(2)}</math> along many directions it can be convenient to remove all symmetries at once.

==Non-linear response of low-dimensional structures: 2D or 1D ==

The Yambo code is a 3D periodic code, so when you want to study a low-dimensional system 2D or 1D, you have to use a supercell approach.
For example to simulate the non-linear response in a 2D-crystal one, can put the crystal in the xy plane and choose a larger distance in the z-direction
in order to reduce the interaction between the periodic replica. Clearly in this case k-point sampling will be only in the x,y directions.
For calculations on this type of systems, these two warnings must be taken into account:

* if you calculate GW corrections or include electron-hole interaction in the linear/non-linear response it is a good idea to use a '''Coulomb cutoff''', similar to the 2D BSE case [http://www.yambo-code.org/wiki/index.php?title=How_to_treat_low_dimensional_systems How to treat low dimensional systems]. You can add the cutoff just adding the "-r" in the input generation.

* the calculation of SHG and THG are always performed respect to the supercell, therefore you have to '''re-scale''' the result to the effective thickness of the layer, usually the inter-layer distance of the corresponding bulk material for 2D systems: <math>\chi_{rescaled}(\omega) = L_z/d_{eff} \cdot \chi(\omega) </math> where <math>L_z</math> is the z-dimension of the supercell, and <math>d_{eff}</math> is the effective thickness of the 2D system.

'''Nota bene''': Non-linear response along non-periodic directions is not implemented in Yambo, the code will print zero in these directions for all response functions. We plan to implement it in a future version of the code.

==References==
<references />

Real time approach to non-linear response (SHG)

2026-07-20T19:52:53Z

Davide: /* Real-time simulation for the SHG */

Real time approach to non-linear response (SHG)

2026-07-20T14:31:16Z

Davide: /* Real-time simulation for the SHG */

== Second Harmonic Generation in AlAs ==
In this tutorial, we will calculate the second harmonic generation of bulk AlAs, the Lumen databases can be downloaded here: [https://media.yambo-code.eu/educational/tutorials/files/AlAs_DBs.tar.gz AlAs_DBs.tar.gz] (10 MB). The first steps of this tutorials are the same as the one on [http://www.yambo-code.eu/wiki/index.php?title=Prerequisites_for_Real_Time_propagation_with_Yambo Prerequisites for Real-Time propagation with Yambo] with the only difference that we will work on AlAs and we will consider an external field in the direction <code> 1.000000 | 1.000000 | 0.000000 |</code>. The DFT input are available here: [https://media.yambo-code.eu/educational/tutorials/files/AlAs_abinit.tar.gz ABINIT] or [https://media.yambo-code.eu/educational/tutorials/files/AlAs_pwscf.tar.gz QuantumEspresso].
You can run the DFT calculation with ABINIT with the command:

abinit < AlAs.in > output_AlAs

or using QuantumEspresso:

pw.x -inp AlAs.scf.in > output_scf
pw.x -inp AlAs.nscf.in > output_scf

and then import the ABINIT wave-function with the command:

a2y -F AlAso_DS2_WFK.nc

and the QuantumEspresso one, in the folder AlAs.save with the command:

p2y
Now we consider an external field in the <code>[1,1,0]</code> direction and remove symmetries not compatible with this field, as explained in the tutorial [http://www.yambo-code.eu/wiki/index.php?title=Prerequisites_for_Real_Time_propagation_with_Yambo Prerequisites for Real-Time propagation with Yambo].

==Real-time simulation for the SHG==

You can generate the input file with the command <code>yambo_nl -u n</code>:

nloptics # [R NL] Non-linear optics
% NLBands
 3 | 6 | # [NL] Bands
%
NLstep= 0.0100 fs # [NL] Real Time step length
NLverbosity= "high" # [NL] Verbosity level (low | high)
NLtime=-1.000000 fs # [NL] Simulation Time
NLintegrator= "INVINT" # [NL] Integrator ("EULEREXP/RK4/RK2EXP/HEUN/INVINT/CRANKNIC")
NLCorrelation= "IPA" # [NL] Correlation ("IPA/HARTREE/TDDFT/LRC/JGM/SEX/HF")
NLLrcAlpha= 0.000000 # [NL] Long Range Correction
% NLEnRange
 1.000000 | 5.000000 | eV # [NL] Energy range
%
NLEnSteps= 5 # [NL] Energy steps
NLDamping= 0.150000 eV # [NL] Damping
% Field1_Dir
 1.000000 | 1.000000 | 0.000000 | # [RT Field1] Versor
%
Field1_kind= "SOFTSIN" # [NL ExtF] Kind(SIN|SOFTSIN|RES|ANTIRES|GAUSS|DELTA|QSSIN)

The line <code>|1.000000 | 1.000000 | 0.000000 | # [RT Field1] Versor</code> referees to the direction of the external field (x,y,0). 
The default parameters of Lumen are already tuned for second-harmonic generation, so the only thing you have to change is the band range, between 3 and 6 and the energy range between 1.0-5.0 eV and the number of energy steps in this interval that we set to 10. Notice that you cannot set to zero the lowest value of the energy range because this will requires a simulation that lasts infinite time, see below. Finally, consider that Yambo performs a separate calculation for each frequency, so if you set many energy steps the computational time grows linearly with this number. 

Notice that we set <code>NLverbosity= "high"</code> in this way the code will produce a file for each laser frequency containing the time dependent polarization. 

Run <code>yambo_nl</code>. Calculations will take about fifteen minutes on a single processor PC, you have time to study the next section that explains how non-linear response is extracted from the real-time simulations.
In order to speed up calculations, you can run them in parallel.

==Non-linear response with Yambo==

In order to calculate the non-linear response, the system is excited with different laser fields with a sinusoidal shape at frequencies <math>\omega_1, \omega_2, .... , \omega_n </math> determined by the parameters [[Variables#NLEnRangev|NLEnRange]] and [[Variables#NLEnSteps|NLEnSteps]]. A dephasing term is added to the Hamiltonian <math>\gamma = </math>[[Variables#NLDamping|NLDamping]] to simulate a finite broadening and to remove the eigenmodes that are excited by the sudden turning on of the external field. After the dephasing time the outgoing signal is analyzed to extract the non-linear coefficients as shown in the figure below:

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

In the report file <code>r_nlinear</code> you will find the length of the dephasing part and of the sampling one:

Dephasing Time [fs]: 52.65695
Sampling Time [fs]: 4.18566
Total simulation time [fs]: 56.84262
The length of the dephasing interval is inversely proportional to the damping term <math>T_{depth} \simeq 1/\gamma</math> while the length of the sampling is dictated by the smallest frequency we are interested in: <math>T_{samp} \simeq 1/\omega_1 </math>. For this reason if <math> \omega_1=0 </math> or <math> \gamma = 0 </math> simulation time goes to infinity. The response at zero frequency can be calculated as limit of small frequency perturbation.
You can have a look to the file <code>src/nloptics/NL_initialize.F</code> to see how simulation lengths are defined.
Calculations can take some time, run them in parallel, the best number of processors = number of frequencies or reduce the number of frequencies step in the SHG.

==Analysis of the results==

In the sampling region we suppose that the polarization can be written as: <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>. We sample the polarization signal at different times and invert the previous equation by truncating the sum at a finite order <ref>[https://arxiv.org/abs/1309.4012 Nonlinear optics from an ab initio approach by means of the dynamical Berry phase: Application to second- and third-harmonic generation in semiconductors], C. Attaccalite and M. Grüning, Phys. Rev. B 88, 235113(2013)</ref>.
This is done with the command <code>ypp_nl -u</code> that automatically produce an input with the correct values:

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

where <code>Xorder</code> is the order where the previous sum is truncated, and the TimeRange specifies the sampling region. Notice that differently from the first tutorial, in this case, we do not need Damping in <code>ypp_nl</code> because we already included it in the real-time dynamics. Run <code>ypp_nl</code> and it will produce a file called <code>o.YPP-X_probe_order_2</code>. This file contains the different components of <math>\chi^2_{xy*} </math> and in particular the columns 6 and 7 correspond to the imaginary and real part of <math>\chi^2_{xyz} </math>. You can plot the result with gnuplot and the command:

p 'o.YPP-X_probe_order_2' u 1:(sqrt($6**2+$7**2)) w lp ps 1.5 pt 7 lw 1.5

the result should look like

[[File:AlAs results.png|700px|center|SHG in AlAs]]

and compare it with the <math> \chi^2_{zxy} </math> calculated with more frequencies and with the converged result on a <math> 18x18x18</math> k-point grid and bands between 2 and 10. In the figure you can find also the comparison with the results of Ref. <ref>Luppi et al. [https://arxiv.org/abs/1006.2649 PRB B, '''82''', 235201(2010)]</ref>.
Notice that the QuantumEspresso results are slightly different from the Abinit ones. This is due to the different pseudo-potential employed in the calculations if pseudopotentials were the same calculation would have been identical. You can download the script to generate this plot and the converged results [http://www.attaccalite.com/lumen/tutorials/AlAs/AlAs_results.tgz here].

'''Notice that''' : the following parameters are not used in the non-linear response analysis: <code> EnRngeRt</code>, <code>ETStpsRt</code>, <code>DampMode</code>, <code>DampFactor</code>. The damping factor is set directly in the real-time simulation.

'''System of units in Non-linear Optics''': notice that in Yambo the non-linear response functions are in the gaussian system of units, for example the <math> \chi^2 (\omega) </math> is in cm/statvolt. In order to convert the non-linear coefficients between different system of units, you can have a look at the [https://www.attaccalite.com/tutorials_yambo/Appendix_C_Boyd.pdf Appendix C] of the book [https://www.elsevier.com/books/nonlinear-optics/boyd/978-0-12-369470-6 "Non-linear optics" by Robert W. Boyd].

More tutorials on non-linear response (Third harmonic generation, parallelization, spin-orbit etc..) with Yambo can be found here: [http://www.attaccalite.com/lumen/tutorials.html More Tutorials]. 

==Analysis of the results using YamboPy==
The analysis of the result performed in the previous section using <code>ypp_nl</code> can be performed using python script [https://github.com/yambo-code/yambopy YamboPy]. 
In YamboPy we implemented all necessary functions to reads the Yambo databases and post-process the results. Here we show how to get the SHG signal from the previous
simulations. We suppose you correctly installed <code>YamboPy</code> on your PC. Go in the folder where you ran non-linear calculation and type in the python:
 
from yambopy import *
NLDB=YamboNLDB()
pol =NLDB.Polarization[0]
time=NLDB.IO_TIME_points

in this way you read all the Non-linear databases. If your runs are in a different folder then the 'SAVE' one you can specify it using the command:
NLDB=YamboNLDB(calc='MYJOB')

now in the array <code>pol</code> you have all the polarization for all laser frequencies in the three Cartesian directions, while the variable <code>time</code> contains all the time series of your simulation.
Now you can get the non-linear response with the command:

SIN = Xn_from_sine(NLDB,X_order=5)
print(SIN)
OUT = SIN.perform_analysis()
SIN.output_analysis(OUT)

this command will perform a Fourier analysis of the results in the same way of <code>ypp_nl</code> and generate new files with all the requested harmonics ''o.YamboPy-X_probe_order_1'', ''o.YamboPy-X_probe_order_2'' etc...
equivalent to the one generated by <code>ypp_nl</code>. Results can be plot with gnuplot using the command:

p 'o.YamboPy-X_probe_order_2' u 1:(sqrt($6**2+$7**2)) w lp ps 1.5 pt 7 lw 1.5

The script that performs non-linear analysis can be found in <code>yambopy/nl/harmonic_analysis.py</code>. We strongly advice you to have a look to this script and modify it according
to your needs and in order to extract other non-linear response functions.

==How to choose the external field direction for the different X2==

In the X2abc response you have three field directions. The two components ''b'' and ''c'' are determined by the external field and the third ''a''
is the direction where you measure the response. 

The external field in Yambo is set using the variable
% Field1_Dir
1.000000 | 0.000000 | 0.000000 | # [NL ExtF] Versor .
%

For example if you set in Yambo field direction as:

1.0 | 0.0 | 0.0 | ---> you calculate all components X2*xx
0.0 | 1.0 | 0.0 | ---> you calculate all components X2*yy
0.0 | 0.0 | 1.0 | ---> you calculate all components X2*zz

the * correspond to the different direction in the response file <code>o.YPP-X_probe_order_2</code>.
In the <code>o.YPP-X_probe_order_2</code> you have seven columns:

The first column is the energy [[omega]],
the 2nd and 3rd columns are the imaginary and real part of X2x**
the 4th and 5th columns are the imaginary and real part of X2y**
the 6th and 7th columns are the imaginary and real part of X2z**

where the ** are the direction of the incoming field, see above.

For example if you want to calculate the <math>\chi^{(2)}_{zxx}</math> you put the external field in the direction <code> 1.0 | 0.0 | 0.0 |</code>
and then plot the column 6 and 7 of the <code>o.YPP-X_probe_order_2</code> file. 

In order to get the other components of the X2 the calculation are more involved. For example if you put field in the

1.0 | 1.0 | 0.0 |

this corresponds to <math>E(t)=\frac{1}{\sqrt{2}} [ E_x(t) + E_y(t) ]</math> and so at the second order you will have <math>E^2(t)=\frac{1}{2} \left [ E^2_x(t) + E^2_y(t) + E_x(t) E_y(t) + E_y(t) E_x(t) \right]</math> 
this means that the X2 in the <code>o.YPP-X_probe_order_2</code> will correspond to a linear combination of the following different response functions. 
For example in the the 2nd and 3rd columns are the imaginary and real part of: 

<math>\chi^{(2)}_{x**}=\frac{1}{2} \left [\chi^2_{xxx}+\chi^{(2)}_{xyy} + \chi^{(2)}_{xxy} + \chi^{(2)}_{xyx} \right ]</math>

the same for the other column replacing the first <code>x</code> for <code>y,z</code>.
If you know that some components are zero this allow you to extract the corresponding off-diagonal part of <math>\chi^{(2)}_{abc}</math> otherwise you have to perform different calculations to get the component you are interested in. 
Now we present '''a simple example of monolayer hBN (m-hBN)'''. In m-hBN there is only a non-zero component of the <math>\chi^{(2)}</math>and all the others are related by symmetry: 
<math>\chi^{(2)}_{yyy}=-\chi^{(2)}_{yxx}=-\chi^{(2)}_{xxy}=-\chi^{(2)}_{xyx}</math> 
so if you put the field in the <code>|1|1|0|</code> direction you will get: 
<math>\chi^{(2)}_{x**} = \frac{1}{2} \left[\chi^{(2)}_{xxy} + \chi^{(2)}_{xyx} \right] = - \chi^{(2)}_{yyy} </math> 
<math>\chi^{(2)}_{y**} =\frac{1}{2} \left[ \chi^{(2)}_{yyy} + \chi^{(2)}_{yxx} \right] = 0.0 </math>
 
In the figure below you can see the different <math>\chi^{(2)}</math> of monolayer hBN.

[[File:Xhi2 xy.png|center| 800px|xhi2_xy]]

as you can see in the second panel, when the field is in the <code>xy</code> you get a linear combination of <math>\chi^{(2)}</math> corresponding to the equation above, in this case the <code>y</code> is zero while the <code>x</code> component is equal to <math>-\chi^{(2)}_{yyy}</math>.

'''Nota bene:''' if you change the direction of the external field you must to remove the corresponding symmetries as explained in the tutorial [[Prerequisites for Real Time_propagation with Yambo#Reduce symmetries real-time linear response tutorial]]. If you have to calculate <math>\chi^{(2)}</math> along many directions it can be convenient to remove all symmetries at once.

==Non-linear response of low-dimensional structures: 2D or 1D ==

The Yambo code is a 3D periodic code, so when you want to study a low-dimensional system 2D or 1D, you have to use a supercell approach.
For example to simulate the non-linear response in a 2D-crystal one, can put the crystal in the xy plane and choose a larger distance in the z-direction
in order to reduce the interaction between the periodic replica. Clearly in this case k-point sampling will be only in the x,y directions.
For calculations on this type of systems, these two warnings must be taken into account:

* if you calculate GW corrections or include electron-hole interaction in the linear/non-linear response it is a good idea to use a '''Coulomb cutoff''', similar to the 2D BSE case [http://www.yambo-code.org/wiki/index.php?title=How_to_treat_low_dimensional_systems How to treat low dimensional systems]. You can add the cutoff just adding the "-r" in the input generation.

* the calculation of SHG and THG are always performed respect to the supercell, therefore you have to '''re-scale''' the result to the effective thickness of the layer, usually the inter-layer distance of the corresponding bulk material for 2D systems: <math>\chi_{rescaled}(\omega) = L_z/d_{eff} \cdot \chi(\omega) </math> where <math>L_z</math> is the z-dimension of the supercell, and <math>d_{eff}</math> is the effective thickness of the 2D system.

'''Nota bene''': Non-linear response along non-periodic directions is not implemented in Yambo, the code will print zero in these directions for all response functions. We plan to implement it in a future version of the code.

==References==
<references />

==Links for schools==
* Back to [[Rome 2023#Tutorials]]
* Back to [[ICTP2020]]

 
{| style="width:100%" border="1"
|style="width:25%; text-align:left"|Prev: [[Linear response from real time simulations]]
|style="width:45%; text-align:center"|Now: [[ICTP2020|ICTP Tutorials]] --> [[Real time approach to non-linear response|Non Linear Response]]
|style="width:50%; text-align:left"|Next: [[Correlation effects in the non-linear response]]
|-
|}

Compile the external libraries without access to internet connection

2026-07-07T14:13:24Z

Davide:

If some external libraries are not provided, the configure script will try to download these libraries and compile them automatically.
However, if you are on a machine without access to the internet, this procedure will fail. Below a solution to this problem, which is common on many HPC machiens.

First of all you have to download the source code on a local machine with access to internet.
You can take for example the bug-fixes from gitlab version
git clone https://gitlab.com/lumen-code/lumen.git lumen-src
cd lumen-src
git checkout bug-fixes
or you can download one of the lumen release via wget
wget https://gitlab.com/lumen-code/lumen/-/archive/2.1.0/lumen-2.1.0.tar.gz
tar -xzf lumen-2.1.0.tar.gz
mv lumen-2.1.0 lumen-src
cd lumen-src

Then you, to download all the libraries which lumen might need, simply use the command (no configure needed)
make download
(note: the download of some packages might fail, you can safely ignore these fails)

Finally you can compress the whole source and copy it on your hpc machine
cd ..
tar -zcvf lumen-src.tar.gz lumen-src
scp lumen-src.tar.gz USERNAME@YOUR.HPC.MACHINE:

After these steps, the automaitic compilation of external libraries will proceed also on the machine without internet connection.

Compile the external libraries without access to internet connection

2026-07-07T12:29:32Z

Davide: removed extra lines not needed

If some external libraries are not provided, the configure script will try to download these libraries and compile them automatically.
However, if you are on a machine without access to the internet, this procedure will fail. Below a solution to this problem, which is common on many HPC machiens.

First of all you have to download the source code on a local machine with access to internet.
You can take for example the bug-fixes from gitlab version
git clone https://gitlab.com/lumen-code/lumen.git lumen-src
cd lumen-src
git checkout bug-fixes
or you can download one of the lumen release via wget
wget https://gitlab.com/lumen-code/lumen/-/archive/2.1.0/lumen-2.1.0.tar.gz
tar -xzf lumen-2.1.0.tar.gz
mv lumen-2.1.0 lumen-src
cd lumen-src

Then you have to run configure and download all the necessary libraries, by doing
make download
(note: the download of some packages might fail, you can safely ignore these fails)

Finally you can compress the whole source and copy it on your hpc machine
cd ..
tar -zcvf lumen-src.tar.gz lumen-src
scp lumen-src.tar.gz USERNAME@YOUR.HPC.MACHINE:

After these steps, the automaitic compilation of external libraries will proceed also on the machine without internet connection.

Install Lumen on IRENE machine (TGCC Joliot Curie)

2026-07-07T12:28:59Z

Davide: library download info moved in a dedicated page

In this page we explain how to configure compile and run Lumen on the IRENE machine on [https://www-hpc.cea.fr/en/complexe/tgcc-JoliotCurie.htm TGCC Joliot Curie].

Since Lumen needs some libraries, and these libraries cannot be downloaded directly from IRENE, first of all you to download them on a local machine.
Please follow the instructions here: [[Compile the external libraries without access to internet connection]]

Now you can copy the file lumen-src.tar.gz on the Irene cluster and use it for the compilation.
Then you load the following modules:

ml feature/intelmpi/mpi_compiler/inteloneapi
module load mpi/intelmpi/20.0.4
module load inteloneapi/24.0.0

and then configure Lambo with the command:

./configure FC=ifx F77=ifx CC=icx MPIF77=mpif90 MPICC=mpicc MPIFC=mpif90 --enable-open-mp --enable-par-linalg --with-scalapack-libs=mkl --with-blacs-libs=mkl

and now you can compile with

make core

You can also compile specific projects for example

make rt-project
make nl-project

then you can submite the job:

ccc_msub job.sh

a typical job file is:

#!/bin/bash
#MSUB -r TEST_lumen
#MSUB -q rome
#MSUB -n 32
#MSUB -A gen13493
#MSUB -m scratch
#MSUB -T 86400
#MSUB -x
#MSUB -@ your_email@your_istitution.com

OMP_NUM_THREADS=1

ml feature/intelmpi/mpi_compiler/inteloneapi
module load mpi/intelmpi/20.0.4
module load inteloneapi/24.0.0

ccc_mprun ./yambo -F my_yambo_input.in

'''N.B.''' We are working on SLEPC and PETSC configuration with Lumen on IRENE, we will update this page soon.

Compile the external libraries without access to internet connection

2026-07-07T12:26:35Z

Davide: Created page with instructions for libraries download

If some external libraries are not provided, the configure script will try to download these libraries and compile them automatically.
However, if you are on a machine without access to the internet, this procedure will fail. Below a solution to this problem, which is common on many HPC machiens.

First of all you have to download the source code on a local machine with access to internet.
You can take for example the bug-fixes from gitlab version
git clone https://gitlab.com/lumen-code/lumen.git lumen-src
cd lumen-src
git checkout bug-fixes
or you can download one of the lumen release via wget
wget https://gitlab.com/lumen-code/lumen/-/archive/2.1.0/lumen-2.1.0.tar.gz
tar -xzf lumen-2.1.0.tar.gz
mv lumen-2.1.0 lumen-src
cd lumen-src

Then you have to run configure and download all the necessary libraries, by doing
make download
(note: the download of some packages might fail, you can safely ignore these fails)

Finally you can compress the whole source and copy it on your hpc machine
cd ..
tar -zcvf lumen-src.tar.gz lumen-src
scp lumen-src.tar.gz USERNAME@YOUR.HPC.MACHINE:

After these steps, the automaitic compilation of external libraries will proceed also on the machine without internet connection

finally you can copy the

Installation

2026-07-07T12:14:30Z

Davide: /* Libraries */

In this page we explain how to install Lumen on different machine with different compiler and configurations.

===General installations===
* [[Install Lumen on Ubuntu/LinuxMint]]
* [[Install Lumen on Ubuntu/LinuxMint with Intel compiler]]
* [[Install Lumen on Ubuntu/LinuxMint with NVfortran compiler]]
* [[Install Lumenon MacOS (ARM)]]
* [[Install Lumen on a Debian machine with some libraries provided via modules]]
* [[Install Lumen on MacOsX (from Yuncheng Mao)]]

===Libraries===
* [[Compile the external libraries without access to internet connection]]
* [[Compile the libraries in Lumen independently]]
* [[Configure library externally]] (this seems a duplicate of the previous one, to check)
* [[Install Lumen with external HDF5 and NetCDF libraries]]
* [[Install Lumen with external chase libraries|Install Lumen with external ChASE libraries]]
* [[Install Lumen with external slepc and petsc libraries|Install Lumen with external Slepc and Petsc libraries]]

===Installations on specific machines===
* [[Install Lumen Claudia cluster in Marseille|Install Lumen on Claudia cluster (Marseille]])
* [[Install Lumen on Mesocentre (Marseille)]]
* [[Install Lumen on IRENE machine (TGCC Joliot Curie)]]
* [[Install Lumen on ISM cluster (Rome)]]
* [[Install Lumen on Leonardo|Install Lumen on Leonardo (Cineca, Bologna)]]
* [[Install Lumen on Tirant|Install Lumen/LetzElPhC/Yambopy on Tirant (Univ. of Valencia Cluster)]]
* [[Install Lumen/LetzElPhC/Yambopy on Lynx(Univ. of Valencia)| Install Lumen on Lynx (Univ. of Valencia Cluster)]]
* [[Install Lumen on Viper (MPCDF, Munich)]]

Installation

2026-07-07T12:14:01Z

Davide: /* Libraries */

In this page we explain how to install Lumen on different machine with different compiler and configurations.

===General installations===
* [[Install Lumen on Ubuntu/LinuxMint]]
* [[Install Lumen on Ubuntu/LinuxMint with Intel compiler]]
* [[Install Lumen on Ubuntu/LinuxMint with NVfortran compiler]]
* [[Install Lumenon MacOS (ARM)]]
* [[Install Lumen on a Debian machine with some libraries provided via modules]]
* [[Install Lumen on MacOsX (from Yuncheng Mao)]]

===Libraries===
* [[Compile the external libraries on machines without access to external links]]
* [[Compile the libraries in Lumen independently]]
* [[Configure library externally]] (this seems a duplicate of the previous one, to check)
* [[Install Lumen with external HDF5 and NetCDF libraries]]
* [[Install Lumen with external chase libraries|Install Lumen with external ChASE libraries]]
* [[Install Lumen with external slepc and petsc libraries|Install Lumen with external Slepc and Petsc libraries]]

===Installations on specific machines===
* [[Install Lumen Claudia cluster in Marseille|Install Lumen on Claudia cluster (Marseille]])
* [[Install Lumen on Mesocentre (Marseille)]]
* [[Install Lumen on IRENE machine (TGCC Joliot Curie)]]
* [[Install Lumen on ISM cluster (Rome)]]
* [[Install Lumen on Leonardo|Install Lumen on Leonardo (Cineca, Bologna)]]
* [[Install Lumen on Tirant|Install Lumen/LetzElPhC/Yambopy on Tirant (Univ. of Valencia Cluster)]]
* [[Install Lumen/LetzElPhC/Yambopy on Lynx(Univ. of Valencia)| Install Lumen on Lynx (Univ. of Valencia Cluster)]]
* [[Install Lumen on Viper (MPCDF, Munich)]]

Linear response from real time simulations (density matrix only)

2026-07-07T10:27:13Z

Davide: /* Circular fields */

== Introduction ==

In this example, we will consider a single layer of hexagonal boron nitride (hBN).
If you didn't before you can 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].
and/or follow the instructions to generate the databases here: [[Prerequisites for Real Time propagation with Yambo]]

You should now be inside the folder <code>hBN-2D-RT/YAMBO/</code>
Before proceeding with the real-time simulations it is useful to compute the Independent Particles (IP) absorption spectrum of hBN along the y direction.
If you didn't before create an Inputs folder and then the input file with the commands
mkdir Inputs
yambo_rt -o c -F Inputs/01_ip.in
and set the proper input parameters
optics # [R OPT] Optics
chi # [R CHI] Dyson equation for Chi.
dipoles # [R ] Compute the dipoles
Chimod= "IP" # [X] IP/Hartree/ALDA/LRC/PF/BSfxc
% QpntsRXd
1 | 1 | # [Xd] Transferred momenta
%
% BndsRnXd
3 | 6 | # [Xd] Polarization function bands
%
% EnRngeXd
0.00000 | 20.00000 | eV # [Xd] Energy range
%
% DmRngeXd
0.10000 | 0.10000 | eV # [Xd] Damping range
%
ETStpsXd= 2001 # [Xd] Total Energy steps
% LongDrXd
0.000000 | 1.000000 | 0.000000 | # [Xd] [cc] Electric Field
%
and then run the code
yambo_rt -F Inputs/01_ip.in -J CHI_IP -C CHI_IP

The resulting spectrum will give an idea of the energy involved in the real-time simulations
gnuplot
gnuplot> plot "CHI_IP/o-CHI_IP.eps_q1_ip" u 1:2 w l
[[File:Independent Particles absorption for hBN.png|thumb|center|900px|Independent Particles absorptions for hBN comuting using 4 bands and with LDA eigenvalues]]

==Real-time dynamics at the independent particles level ==

Now you should enter the folder <code>hBN-2D-RT/YAMBO/FixSymm</code>
Let us first create a folder for the input files to be run with <code>yambo_rt</code>.

cd FixSymm
mkdir Inputs_rt

In order to calculate linear-response in real-time, we will perturb the system with a delta function in time external field.

Use the command
yambo_rt -n p -potential ip -F Inputs_rt/01_td_ip.in

to generate the input:

negf # [R] Real-Time dynamics
RT_Threads=0 # [OPENMP/RT] Number of threads for real-time
HXC_Potential= "IP" # [SC] SC HXC Potential
% RTBands
3 | 6 | # [RT] Bands
%
Integrator= "EULER RK2" # [RT] Integrator. Use keywords space separated ( "EULER/EXPn/INV" "SIMPLE/RK2/RK4/HEUN" "RWA")
PhLifeTime= 0.000000 fs # [RT] Dephasing Time
RTstep=10.000000 as # [RT] Real Time step length
NETime= 55.00000 fs # [RT] Simulation Time
% IOtime
1.00 | 5.00 | 0.10 | fs # [RT] Time between to consecutive I/O (OBSERVABLEs,CARRIERs - GF - OUTPUT)
%
Field1_Freq= 0.00 eV # [RT Field1] Frequency
Field1_Int=1.E3 kWLm2 # [RT Field1] Intensity
Field1_Width= 0.000000 fs # [RT Field1] Width
Field1_kind= "DELTA" # [RT Field1] Kind(SIN|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.000000fs # [RT Field1] Initial Time

Set the field direction along y, the field type to <code>DELTA</code>, the length of the simulation to 55 fs, number of bands from 3 to 6, dephasing to zero and the field intensity to 1.E3 [kW/cm2].
A small intensity is needed to ensure that we remain in the perturbative regime and that the response is dominated by linear term.
The <code>yambo_rt</code> is optimized for TD-SEX (or even more sophisticated calculations).

Now run the simulation:
yambo_rt -F Inputs_rt/01_td_ip.in -J TD-IP_rt -C TD-IP_rt
The run produces, besides the standard report and (eventually) log files of yambo, 3 output files:
TD-IP_rt/01_td_ip.in_TD-IP_rt
TD-IP_rt/r-TD-IP_rt_negf
TD-IP_rt/o-TD-IP_rt.polarization
TD-IP_rt/o-TD-IP_rt.external_field
TD-IP_rt/o-TD-IP_rt.current
Moreover it also makes a copy of the input file. If the run is taking to long you can open the copy of the input file <code>TD-IP_rt/01_td_ip.in_TD-IP_rt</code> and add there the string
STOP_NOW
It will finish the simulation in a proper way.

You can now plot the resulting time-dependent polarization and obtain something like this
gnuplot
plot "TD-IP_rt/o-TD-IP_rt.polarization" u 1:3 w l

[[File:IP Polarization along y direction.png|thumb|center|900px|Time dependent polarization generated with a TD-IP run using yambo_rt]]

The polarization is oscillating very quickly and with the time step for the output file we have chosen (100 as) we can hardly resolve the oscillations. The slowest oscillation is expected at the first absorption peak, which is located slightly above 4 eV. This corresponds to a period of about 1 fs. Higher energy transitions are however also activated by the delta like pulse, which spans all frequencies (remember that the Fourier transform of a delta function is a constant).

The I/O time, which is negligible in such calculations, is less optimized and becomes the most time demanding step here.
We can overcome this setting the IOtime for the polarization to 50 as (i.e. each 5 time steps). After the run you can have a look at the timing report and you will see that indeed most of the timing was spent in the I/O:
RT databases IO : 6.8745 s CPU ( 5501 calls, 0.0012 s avg)



One can test the performance of different integrators. For example, by setting in input <code>Integrator= "INV RK2"</code>, <code>yambo_rt</code> uses a similar integration scheme as <code>yambo_nl</code> (the latter is based on a Berry-phase approach).

=== Circularly polarized fields ===

You can also perform a real time simulation under the action of a circularly polarized field. You just need to modify the input as follows

Field1_pol= "circular" # [RT Field1] Pol(linear|circular)
% Field1_Dir
0.000000 | 1.000000 | 0.000000 | # [RT Field1] Versor
%
% Field1_Dir_Cir
0.000000 | 0.000000 | 1.000000 | # [RT Field1] Versor
%

Note. This will translate in the superposition of two fields shifted in time by 90 degree. It will result into a circularly polarized field if the two are orthogonal, as in the above input. You can change the sign of the Field1_Dir_Cir variable to switch between plus and minus polarization. Since the simulation is performed in presence of two external fields, the SAVE folder must be compatible with both fields, i.e. you must put both field directions in input when removing symmetries.

==Results Analysis==

We can now proceed to the Fourier transform of the polarization to obtain the dielectric function

We can use
ypp_rt -n X -F Inputs_rt/ypp_abs.in

to generate the input file

TDplots # [R] TD observables plot
RT_X # [R] Response functions Post-Processing
RTtime # [R] Post-Processing kind: function of time
Xorder= 1 # Max order of the response/exc functions
% EnRngeRt
0.00000 | 10.00000 | eV # Energy range
%
ETStpsRt= 1001 # Total Energy steps
% TimeRange
-1.000000 |-1.000000 | fs # Time-window where processing is done
%
DampMode= "LORENTZIAN" # Damping type ( NONE | LORENTZIAN | GAUSSIAN )
DampFactor= 0.10000 eV # Damping parameter
#SkipOBS_IO # Do not dump on file the RT observables (P(t),J(t),D(t)...)

where we set a Lorentzian smearing corresponding to 0.1 eV. Notice that due to the finite time of our simulation smearing is always necessary to Fourier transform the result.
Then we run
ypp_rt -F Inputs_rt/ypp_abs.in -J TD-IP_rt -C TD-IP_rt
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_rt.YPP-eps_along_E
o-TD-IP_rt.YPP-eels_along_E
o-TD-IP_rt.YPP-polarization
o-TD-IP_rt.YPP-current
o-TD-IP_rt.YPP-E_frequency
inside the folder <code>TD-IP_rt</code>

First let's have a look to the polarization.
gnuplot
plot "TD-IP_rt/o-TD-IP_rt.polarization" u 1:3 w l, "TD-IP_rt/o-TD-IP_rt.YPP-polarization" u 1:3 w l
We clearly see the effect of the Lorentzian damping added to the time dependent polarization
[[File:IP Polarization along y direction damped.png|thumb|center|900px|Time dependent polarization generated with the ypp_rt code. Comparing the damped polarization with the one produced by yambo_rt]]

Finally, we can plot the imaginary part of the dielectric function, which is obtained from the Fourier transform of the polarization.

Now we can plot the dielectric constant and compare it with the linear response:
gnuplot
plot "../CHI_IP/o-CHI_IP.eps_q1_ip" u 1:2 w l, "TD-IP_rt/o-TD-IP_rt.YPP-eps_along_E" u 1:2 w l

[[File:Bn optics.png|thumb|900px|center|Imaginary part of the dielectric constant. Comparison between standard IP calculation (yambo), linear response based on the density matrix (yambo_rt,this tutorial) and Berry-phase approach (yambo_nl).]]

==Links==
* Back to [[Rome 2023#Tutorials]]
* On to [[Real time Bethe-Salpeter Equation (density matrix only)|TD-BSE tutorial]]

Linear response from real time simulations (density matrix only)

2026-07-07T10:26:46Z

Davide: Added subsection on circluarly polarized fields

== Introduction ==

In this example, we will consider a single layer of hexagonal boron nitride (hBN).
If you didn't before you can 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].
and/or follow the instructions to generate the databases here: [[Prerequisites for Real Time propagation with Yambo]]

You should now be inside the folder <code>hBN-2D-RT/YAMBO/</code>
Before proceeding with the real-time simulations it is useful to compute the Independent Particles (IP) absorption spectrum of hBN along the y direction.
If you didn't before create an Inputs folder and then the input file with the commands
mkdir Inputs
yambo_rt -o c -F Inputs/01_ip.in
and set the proper input parameters
optics # [R OPT] Optics
chi # [R CHI] Dyson equation for Chi.
dipoles # [R ] Compute the dipoles
Chimod= "IP" # [X] IP/Hartree/ALDA/LRC/PF/BSfxc
% QpntsRXd
1 | 1 | # [Xd] Transferred momenta
%
% BndsRnXd
3 | 6 | # [Xd] Polarization function bands
%
% EnRngeXd
0.00000 | 20.00000 | eV # [Xd] Energy range
%
% DmRngeXd
0.10000 | 0.10000 | eV # [Xd] Damping range
%
ETStpsXd= 2001 # [Xd] Total Energy steps
% LongDrXd
0.000000 | 1.000000 | 0.000000 | # [Xd] [cc] Electric Field
%
and then run the code
yambo_rt -F Inputs/01_ip.in -J CHI_IP -C CHI_IP

The resulting spectrum will give an idea of the energy involved in the real-time simulations
gnuplot
gnuplot> plot "CHI_IP/o-CHI_IP.eps_q1_ip" u 1:2 w l
[[File:Independent Particles absorption for hBN.png|thumb|center|900px|Independent Particles absorptions for hBN comuting using 4 bands and with LDA eigenvalues]]

==Real-time dynamics at the independent particles level ==

Now you should enter the folder <code>hBN-2D-RT/YAMBO/FixSymm</code>
Let us first create a folder for the input files to be run with <code>yambo_rt</code>.

cd FixSymm
mkdir Inputs_rt

In order to calculate linear-response in real-time, we will perturb the system with a delta function in time external field.

Use the command
yambo_rt -n p -potential ip -F Inputs_rt/01_td_ip.in

to generate the input:

negf # [R] Real-Time dynamics
RT_Threads=0 # [OPENMP/RT] Number of threads for real-time
HXC_Potential= "IP" # [SC] SC HXC Potential
% RTBands
3 | 6 | # [RT] Bands
%
Integrator= "EULER RK2" # [RT] Integrator. Use keywords space separated ( "EULER/EXPn/INV" "SIMPLE/RK2/RK4/HEUN" "RWA")
PhLifeTime= 0.000000 fs # [RT] Dephasing Time
RTstep=10.000000 as # [RT] Real Time step length
NETime= 55.00000 fs # [RT] Simulation Time
% IOtime
1.00 | 5.00 | 0.10 | fs # [RT] Time between to consecutive I/O (OBSERVABLEs,CARRIERs - GF - OUTPUT)
%
Field1_Freq= 0.00 eV # [RT Field1] Frequency
Field1_Int=1.E3 kWLm2 # [RT Field1] Intensity
Field1_Width= 0.000000 fs # [RT Field1] Width
Field1_kind= "DELTA" # [RT Field1] Kind(SIN|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.000000fs # [RT Field1] Initial Time

Set the field direction along y, the field type to <code>DELTA</code>, the length of the simulation to 55 fs, number of bands from 3 to 6, dephasing to zero and the field intensity to 1.E3 [kW/cm2].
A small intensity is needed to ensure that we remain in the perturbative regime and that the response is dominated by linear term.
The <code>yambo_rt</code> is optimized for TD-SEX (or even more sophisticated calculations).

Now run the simulation:
yambo_rt -F Inputs_rt/01_td_ip.in -J TD-IP_rt -C TD-IP_rt
The run produces, besides the standard report and (eventually) log files of yambo, 3 output files:
TD-IP_rt/01_td_ip.in_TD-IP_rt
TD-IP_rt/r-TD-IP_rt_negf
TD-IP_rt/o-TD-IP_rt.polarization
TD-IP_rt/o-TD-IP_rt.external_field
TD-IP_rt/o-TD-IP_rt.current
Moreover it also makes a copy of the input file. If the run is taking to long you can open the copy of the input file <code>TD-IP_rt/01_td_ip.in_TD-IP_rt</code> and add there the string
STOP_NOW
It will finish the simulation in a proper way.

You can now plot the resulting time-dependent polarization and obtain something like this
gnuplot
plot "TD-IP_rt/o-TD-IP_rt.polarization" u 1:3 w l

[[File:IP Polarization along y direction.png|thumb|center|900px|Time dependent polarization generated with a TD-IP run using yambo_rt]]

The polarization is oscillating very quickly and with the time step for the output file we have chosen (100 as) we can hardly resolve the oscillations. The slowest oscillation is expected at the first absorption peak, which is located slightly above 4 eV. This corresponds to a period of about 1 fs. Higher energy transitions are however also activated by the delta like pulse, which spans all frequencies (remember that the Fourier transform of a delta function is a constant).

The I/O time, which is negligible in such calculations, is less optimized and becomes the most time demanding step here.
We can overcome this setting the IOtime for the polarization to 50 as (i.e. each 5 time steps). After the run you can have a look at the timing report and you will see that indeed most of the timing was spent in the I/O:
RT databases IO : 6.8745 s CPU ( 5501 calls, 0.0012 s avg)



One can test the performance of different integrators. For example, by setting in input <code>Integrator= "INV RK2"</code>, <code>yambo_rt</code> uses a similar integration scheme as <code>yambo_nl</code> (the latter is based on a Berry-phase approach).

=== Circular fields ===

You can also perform a real time simulation under the action of a circular field. You just need to modify the input as follows

Field1_pol= "circular" # [RT Field1] Pol(linear|circular)
% Field1_Dir
0.000000 | 1.000000 | 0.000000 | # [RT Field1] Versor
%
% Field1_Dir_Cir
0.000000 | 0.000000 | 1.000000 | # [RT Field1] Versor
%

Note. This will translate in the superposition of two fields shifted in time by 90 degree. It will result into a circularly polarized field if the two are orthogonal, as in the above input. You can change the sign of the Field1_Dir_Cir variable to switch between plus and minus polarization. Since the simulation is performed in presence of two external fields, the SAVE folder must be compatible with both fields, i.e. you must put both field directions in input when removing symmetries.

==Results Analysis==

We can now proceed to the Fourier transform of the polarization to obtain the dielectric function

We can use
ypp_rt -n X -F Inputs_rt/ypp_abs.in

to generate the input file

TDplots # [R] TD observables plot
RT_X # [R] Response functions Post-Processing
RTtime # [R] Post-Processing kind: function of time
Xorder= 1 # Max order of the response/exc functions
% EnRngeRt
0.00000 | 10.00000 | eV # Energy range
%
ETStpsRt= 1001 # Total Energy steps
% TimeRange
-1.000000 |-1.000000 | fs # Time-window where processing is done
%
DampMode= "LORENTZIAN" # Damping type ( NONE | LORENTZIAN | GAUSSIAN )
DampFactor= 0.10000 eV # Damping parameter
#SkipOBS_IO # Do not dump on file the RT observables (P(t),J(t),D(t)...)

where we set a Lorentzian smearing corresponding to 0.1 eV. Notice that due to the finite time of our simulation smearing is always necessary to Fourier transform the result.
Then we run
ypp_rt -F Inputs_rt/ypp_abs.in -J TD-IP_rt -C TD-IP_rt
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_rt.YPP-eps_along_E
o-TD-IP_rt.YPP-eels_along_E
o-TD-IP_rt.YPP-polarization
o-TD-IP_rt.YPP-current
o-TD-IP_rt.YPP-E_frequency
inside the folder <code>TD-IP_rt</code>

First let's have a look to the polarization.
gnuplot
plot "TD-IP_rt/o-TD-IP_rt.polarization" u 1:3 w l, "TD-IP_rt/o-TD-IP_rt.YPP-polarization" u 1:3 w l
We clearly see the effect of the Lorentzian damping added to the time dependent polarization
[[File:IP Polarization along y direction damped.png|thumb|center|900px|Time dependent polarization generated with the ypp_rt code. Comparing the damped polarization with the one produced by yambo_rt]]

Finally, we can plot the imaginary part of the dielectric function, which is obtained from the Fourier transform of the polarization.

Now we can plot the dielectric constant and compare it with the linear response:
gnuplot
plot "../CHI_IP/o-CHI_IP.eps_q1_ip" u 1:2 w l, "TD-IP_rt/o-TD-IP_rt.YPP-eps_along_E" u 1:2 w l

[[File:Bn optics.png|thumb|900px|center|Imaginary part of the dielectric constant. Comparison between standard IP calculation (yambo), linear response based on the density matrix (yambo_rt,this tutorial) and Berry-phase approach (yambo_nl).]]

==Links==
* Back to [[Rome 2023#Tutorials]]
* On to [[Real time Bethe-Salpeter Equation (density matrix only)|TD-BSE tutorial]]

Real-time dynamics with external field read from file

2026-07-07T10:18:25Z

Davide:

[[File:Modulated-Gaussian-pulse.png|right|500px |external field from file]]

'''In this tutorial we will see how to generate an external field in a format compatible with Lumen and then use it to run real-time simulation (see for instance Ref. <ref name="sangalli"/>).'''

We will also prove a simple python script to do it, that you can easily modify to create your own external field.

The tutorial is divided in two sections: 1) generate an external field, either with ''YamboPy'' (1.a) or with ''ypp_rt'' (1.b); 2) use the external field in the real-time calculations.

== Generate external field ==

The idea is to generate a file, which we call here <code>EXTFIELD.time</code>, containing the time dependent vector . <code>EXTFIELD.time</code> will be later used in a real-time simulation.
The structure of the file will be the following

8001 5.00000035E-03
0.00000000 0.00000000 8.00802198E-08 5.33216726E-03
5.00000035E-03 0.00000000 1.10453682E-03 5.34489006E-03
1.00000007E-02 0.00000000 2.20612506E-03 5.30774519E-03
1.50000006E-02 0.00000000 3.29522160E-03 5.22070238E-03
2.00000014E-02 0.00000000 4.36120294E-03 5.08422730E-03
2.50000022E-02 0.00000000 5.39377611E-03 4.89927875E-03
3.00000012E-02 0.00000000 6.38338085E-03 4.66719037E-03
3.50000039E-02 0.00000000 7.32029136E-03 4.38982388E-03
4.00000028E-02 0.00000000 8.19511060E-03 4.06956300E-03
4.50000018E-02 0.00000000 8.99977796E-03 3.70889902E-03
5.00000045E-02 0.00000000 9.72598325E-03 3.31101404E-03
5.49999997E-02 0.00000000 1.03663849E-02 2.87936884E-03
6.00000024E-02 0.00000000 1.09142270E-02 2.41791131E-03
6.50000051E-02 0.00000000 1.13641499E-02 1.93037838E-03
7.00000077E-02 0.00000000 1.17109315E-02 1.42144389E-03
7.50000030E-02 0.00000000 1.19506177E-02 8.95809906E-04
8.00000057E-02 0.00000000 1.20803779E-02 3.57942627E-04
8.50000083E-02 0.00000000 1.20981205E-02 -1.87048965E-04
9.00000036E-02 0.00000000 1.20029533E-02 -7.33944471E-04
9.50000063E-02 0.00000000 1.17948940E-02 -1.27798819E-03
0.100000009 0.00000000 1.14751300E-02 -1.81389856E-03
...

The first column contains the time, then in order <math>a(t)</math>, <math>a'(t)</math>, and <math>a''(t)</math>. Befor this info there is a row indicating the length of the file (number of time steps), and the time step duration (<math>dt_f</math>).
Since in the simulation a discrete time step (<math>dt</math>) is also used, <code>EXTFIELD.time</code> must contain all the time steps needed. When generating an external field, care must be taken in choosing the time step to ensure compatibility with the integrator.

* EULER and INVINT are single-time integrators; therefore, the field time step must be equal to the one used in the simulation or an integer fraction of it: <math>dt_f=dt/n</math>.
* CRANKNIC, RK2EXP and RK2 are two-time integrators; the field time step should be half the one used in the simulation or an integer fraction of it: <math>dt_f=dt/(2n)</math>.
* RK4 and RK4EXP are four-time integrators that requires 1/2 and 1/3 of the time-step, therefore the field time-step should be one sixth of the one used in the simulation or and integer fraction of it: <math>dt_f=dt/(6n)</math>.

The code will check that these prescriptions are respected before start simulations. In practice, if you generate a field with <math>dt_f=dt/(6n)</math> this will be ok for all available integrators.

=== Generate external field with YamboPy ===

In YamboPy you can find a tutorial to generate external field in the folder <code>yambopy/tutorial/external_field/generate_field.py</code>You can define your external field in this way:

def gauss_field(t_steps):
# In this example we generate a guassian external field
sigma=5.0*fs2aut # Field width
T_0 =3.0*sigma
Expf =np.exp(-(t_steps-T_0)**2/(2.0*sigma**2) )
a_t =sigma*np.sqrt(math.pi/2.0)*(special.erf((t_steps-T_0)/(sigma*np.sqrt(2.0)))+1.0)
ap_t =Expf
app_t=-Expf*(t_steps-T_0)/sigma**2
return a_t,ap_t,app_t

Just a function that generate <math>a(t)</math>, <math>a'(t)</math>, and <math>a''(t)</math> that will be then read from Yambo_rt/Yambo_nl. You can add you function to <code>generate_field.py</code>scritp or use the two that are already implemented

SIN and GAUSS fields. Notice that the electric field correspond to <math>E(t) = -E_{amp} a'(t)</math>.

The field amplitude (intensity) and versor are not required because it will be read from the Yambo input

'''Warning:''' only linear polarized external field can generated with YamboPY (to be generlized)

=== Generate external field with ypp_rt ===
With <code>ypp_rt</code> one can generate all external fields that implemented in the Yambo code in the file <code>src/modules/mod_fields.F</code>. This can be usefull to visualize external fields used in Yambo without the need of running a real simulation. However if you want to add a new external field you have to code it in the file <code>src/modules/mod_fields.F</code>. The command to generate external field on file is: <code>ypp_rt -rtplot f</code>. This will generate the following input file:

TDplots # [R] TD observables plot
RTtime # [R] Post-Processing kind: function of time
RTfields # [R] Analize time-dependent fields
TimeStep= 1.6666667 as # Time step
% TimeRange
0.000000 |100.000000 | fs # Time-window where processing is done
%
% EnRngeRt
0.00000 | 20.00000 | eV # Energy range
%
ETStpsRt= 200 # Total Energy steps
ChirpFac= 0.000000 fs # Prefactor for linear chirping in frequency space
Field1_Freq= 2.000000 eV # [RT Field1] Frequency
Field1_Int= 1000.000000 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 | 0.000000 | 0.000000 | # [RT Field1] Versor
%
Field1_Tstart= 0.010000 fs # [RT Field1] Initial Time
In the input file you have to specify a series of parameters, in red: the time step (<code>TimeSpep</code>), the time range (<code>TimeRange</code>), the field kind (<code>Field1_kind</code>), its parameters (<code>Field1_Width</code>, <code>Field1_Freq</code>, etc...) and finally the starting time (<code>Field1_Tstart</code>), usually the first step of the dynamics. Notice that Field direction (<code>Field1_Dir</code>) field intensity (<code>Field1_Int</code>) are not used in the external field generation but read from Yambo input. 
Then you run <code>ypp_rt</code> with this input and it will produce a series of files: <code>YPP_EXTFIELD1_P1.time</code>, <code>YPP_EXTFIELD1_P2.time</code>, <code>YPP_EXTFIELD1_P1.freq</code>, <code>YPP_EXTFIELD1_P2.freq</code><code>.</code>

These files contain the external field and its Fourier transform. In the next section we will see how to use them in real-time dynamics.

== Use generated external field with yambo_rt/yambo_nl ==
In yambo_rt/yambo_nl you can tell to the code to read an external field from file. For example in yambo_nl you generate your input wil the command <code>yambo_nl -u p</code>

nloptics # [R] Non-linear spectroscopy
% NLBands
7 | 10 | # [NL] Bands range
%
NLverbosity= "high" # [NL] Verbosity level (low | high)
NLstep= 10.0 as # [NL] Real Time step length
NLtime=120.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/LSEX/LHF")
NLLrcAlpha= 0.000000 # [NL] Long Range Correction
NLDamping= 0.200000 eV # [NL] Damping (or dephasing)
RADLifeTime=-1.000000 fs # [RT] Radiative life-time (if negative RADLifeTime=Phase_LifeTime)
#EvalCurrent # [NL] Evaluate the current
HARRLvcs= 11529 RL # [HA] Hartree RL components
EXXRLvcs= 11529 RL # [XX] Exchange RL components
CORRLvcs= 11529 RL # [GW] Correlation RL components
Field1_Freq= 0.100000 eV # [RT Field1] Frequency
Field1_Int= 1000.00 kWLm2 # [RT Field1] Intensity
Field1_Width= 0.000000 fs # [RT Field1] Width
Field1_kind= "FROM_FILE ./YPP_EXTFIELD1_P1.time" # [RT Field1] Kind(SIN|COS|RES|ANTIRES|GAUSS|DELTA|QSSIN)
Field1_pol= "linear" # [RT Field1] Pol(linear|circular)
% Field1_Dir
 1.000000 | 0.000000 | 0.000000 | # [RT Field1] Versor
%
You have to specify the file kind <code>FROM_FILE</code> and give the path of the file containing the external field. Do not forget also to use the same <code>NLtime</code> and in the <code>ypp</code> and <code>yambo</code> inputs. Notice that field intensity and field versor is specified in the <code>Yambo</code> input only.

Notice that <code>NLtime</code> in the simulation can be larger then the one used in the external field generation, the code will assume that field that after a given time the field is zero.

== References ==
<references>
<ref name="sangalli" >D. Sangalli, ''Excitons and carriers in transient absorption and time-resolved ARPES spectroscopy: An ab initio approach'', [https://arxiv.org/pdf/2106.03592 Phys. Rev. Materials 5, 083803, (2021)]. </ref>
</references>

Real-time dynamics with external field read from file

2026-07-07T10:17:18Z

Davide: Removed "with" from section title

[[File:Modulated-Gaussian-pulse.png|right|500px |external field from file]]

'''In this tutorial we will see how to generate an external field in a format compatible with Lumen and then use it to run real-time simulation (see for instance Ref. <ref name="sangalli"/>).'''

We will also prove a simple python script to do it, that you can easily modify to create your own external field.

The tutorial is divided in three sections: 1) generate external field with ''YamboPy''; 2) generate external field with ''ypp_rt''; 3) use the external field in the real-time calculations.

== Generate external field ==

The idea is to generate a file, which we call here <code>EXTFIELD.time</code>, containing the time dependent vector . <code>EXTFIELD.time</code> will be later used in a real-time simulation.
The structure of the file will be the following

8001 5.00000035E-03
0.00000000 0.00000000 8.00802198E-08 5.33216726E-03
5.00000035E-03 0.00000000 1.10453682E-03 5.34489006E-03
1.00000007E-02 0.00000000 2.20612506E-03 5.30774519E-03
1.50000006E-02 0.00000000 3.29522160E-03 5.22070238E-03
2.00000014E-02 0.00000000 4.36120294E-03 5.08422730E-03
2.50000022E-02 0.00000000 5.39377611E-03 4.89927875E-03
3.00000012E-02 0.00000000 6.38338085E-03 4.66719037E-03
3.50000039E-02 0.00000000 7.32029136E-03 4.38982388E-03
4.00000028E-02 0.00000000 8.19511060E-03 4.06956300E-03
4.50000018E-02 0.00000000 8.99977796E-03 3.70889902E-03
5.00000045E-02 0.00000000 9.72598325E-03 3.31101404E-03
5.49999997E-02 0.00000000 1.03663849E-02 2.87936884E-03
6.00000024E-02 0.00000000 1.09142270E-02 2.41791131E-03
6.50000051E-02 0.00000000 1.13641499E-02 1.93037838E-03
7.00000077E-02 0.00000000 1.17109315E-02 1.42144389E-03
7.50000030E-02 0.00000000 1.19506177E-02 8.95809906E-04
8.00000057E-02 0.00000000 1.20803779E-02 3.57942627E-04
8.50000083E-02 0.00000000 1.20981205E-02 -1.87048965E-04
9.00000036E-02 0.00000000 1.20029533E-02 -7.33944471E-04
9.50000063E-02 0.00000000 1.17948940E-02 -1.27798819E-03
0.100000009 0.00000000 1.14751300E-02 -1.81389856E-03
...

The first column contains the time, then in order <math>a(t)</math>, <math>a'(t)</math>, and <math>a''(t)</math>. Befor this info there is a row indicating the length of the file (number of time steps), and the time step duration (<math>dt_f</math>).
Since in the simulation a discrete time step (<math>dt</math>) is also used, <code>EXTFIELD.time</code> must contain all the time steps needed. When generating an external field, care must be taken in choosing the time step to ensure compatibility with the integrator.

* EULER and INVINT are single-time integrators; therefore, the field time step must be equal to the one used in the simulation or an integer fraction of it: <math>dt_f=dt/n</math>.
* CRANKNIC, RK2EXP and RK2 are two-time integrators; the field time step should be half the one used in the simulation or an integer fraction of it: <math>dt_f=dt/(2n)</math>.
* RK4 and RK4EXP are four-time integrators that requires 1/2 and 1/3 of the time-step, therefore the field time-step should be one sixth of the one used in the simulation or and integer fraction of it: <math>dt_f=dt/(6n)</math>.

The code will check that these prescriptions are respected before start simulations. In practice, if you generate a field with <math>dt_f=dt/(6n)</math> this will be ok for all available integrators.

=== Generate external field with YamboPy ===

In YamboPy you can find a tutorial to generate external field in the folder <code>yambopy/tutorial/external_field/generate_field.py</code>You can define your external field in this way:

def gauss_field(t_steps):
# In this example we generate a guassian external field
sigma=5.0*fs2aut # Field width
T_0 =3.0*sigma
Expf =np.exp(-(t_steps-T_0)**2/(2.0*sigma**2) )
a_t =sigma*np.sqrt(math.pi/2.0)*(special.erf((t_steps-T_0)/(sigma*np.sqrt(2.0)))+1.0)
ap_t =Expf
app_t=-Expf*(t_steps-T_0)/sigma**2
return a_t,ap_t,app_t

Just a function that generate <math>a(t)</math>, <math>a'(t)</math>, and <math>a''(t)</math> that will be then read from Yambo_rt/Yambo_nl. You can add you function to <code>generate_field.py</code>scritp or use the two that are already implemented

SIN and GAUSS fields. Notice that the electric field correspond to <math>E(t) = -E_{amp} a'(t)</math>.

The field amplitude (intensity) and versor are not required because it will be read from the Yambo input

'''Warning:''' only linear polarized external field can generated with YamboPY (to be generlized)

=== Generate external field with ypp_rt ===
With <code>ypp_rt</code> one can generate all external fields that implemented in the Yambo code in the file <code>src/modules/mod_fields.F</code>. This can be usefull to visualize external fields used in Yambo without the need of running a real simulation. However if you want to add a new external field you have to code it in the file <code>src/modules/mod_fields.F</code>. The command to generate external field on file is: <code>ypp_rt -rtplot f</code>. This will generate the following input file:

TDplots # [R] TD observables plot
RTtime # [R] Post-Processing kind: function of time
RTfields # [R] Analize time-dependent fields
TimeStep= 1.6666667 as # Time step
% TimeRange
0.000000 |100.000000 | fs # Time-window where processing is done
%
% EnRngeRt
0.00000 | 20.00000 | eV # Energy range
%
ETStpsRt= 200 # Total Energy steps
ChirpFac= 0.000000 fs # Prefactor for linear chirping in frequency space
Field1_Freq= 2.000000 eV # [RT Field1] Frequency
Field1_Int= 1000.000000 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 | 0.000000 | 0.000000 | # [RT Field1] Versor
%
Field1_Tstart= 0.010000 fs # [RT Field1] Initial Time
In the input file you have to specify a series of parameters, in red: the time step (<code>TimeSpep</code>), the time range (<code>TimeRange</code>), the field kind (<code>Field1_kind</code>), its parameters (<code>Field1_Width</code>, <code>Field1_Freq</code>, etc...) and finally the starting time (<code>Field1_Tstart</code>), usually the first step of the dynamics. Notice that Field direction (<code>Field1_Dir</code>) field intensity (<code>Field1_Int</code>) are not used in the external field generation but read from Yambo input. 
Then you run <code>ypp_rt</code> with this input and it will produce a series of files: <code>YPP_EXTFIELD1_P1.time</code>, <code>YPP_EXTFIELD1_P2.time</code>, <code>YPP_EXTFIELD1_P1.freq</code>, <code>YPP_EXTFIELD1_P2.freq</code><code>.</code>

These files contain the external field and its Fourier transform. In the next section we will see how to use them in real-time dynamics.

== Use generated external field with yambo_rt/yambo_nl ==
In yambo_rt/yambo_nl you can tell to the code to read an external field from file. For example in yambo_nl you generate your input wil the command <code>yambo_nl -u p</code>

nloptics # [R] Non-linear spectroscopy
% NLBands
7 | 10 | # [NL] Bands range
%
NLverbosity= "high" # [NL] Verbosity level (low | high)
NLstep= 10.0 as # [NL] Real Time step length
NLtime=120.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/LSEX/LHF")
NLLrcAlpha= 0.000000 # [NL] Long Range Correction
NLDamping= 0.200000 eV # [NL] Damping (or dephasing)
RADLifeTime=-1.000000 fs # [RT] Radiative life-time (if negative RADLifeTime=Phase_LifeTime)
#EvalCurrent # [NL] Evaluate the current
HARRLvcs= 11529 RL # [HA] Hartree RL components
EXXRLvcs= 11529 RL # [XX] Exchange RL components
CORRLvcs= 11529 RL # [GW] Correlation RL components
Field1_Freq= 0.100000 eV # [RT Field1] Frequency
Field1_Int= 1000.00 kWLm2 # [RT Field1] Intensity
Field1_Width= 0.000000 fs # [RT Field1] Width
Field1_kind= "FROM_FILE ./YPP_EXTFIELD1_P1.time" # [RT Field1] Kind(SIN|COS|RES|ANTIRES|GAUSS|DELTA|QSSIN)
Field1_pol= "linear" # [RT Field1] Pol(linear|circular)
% Field1_Dir
 1.000000 | 0.000000 | 0.000000 | # [RT Field1] Versor
%
You have to specify the file kind <code>FROM_FILE</code> and give the path of the file containing the external field. Do not forget also to use the same <code>NLtime</code> and in the <code>ypp</code> and <code>yambo</code> inputs. Notice that field intensity and field versor is specified in the <code>Yambo</code> input only.

Notice that <code>NLtime</code> in the simulation can be larger then the one used in the external field generation, the code will assume that field that after a given time the field is zero.

== References ==
<references>
<ref name="sangalli" >D. Sangalli, ''Excitons and carriers in transient absorption and time-resolved ARPES spectroscopy: An ab initio approach'', [https://arxiv.org/pdf/2106.03592 Phys. Rev. Materials 5, 083803, (2021)]. </ref>
</references>

Lumen release notes

2026-07-06T13:29:03Z

Davide: /* Lumen Releases */

This page contains the list of the releases of the lumen code.

Starting from version 2.0, the release notes are available in a dedicated page on the wiki, while the source code is available on gitlab.

== Lumen Releases ==

* Lumen 2.1: [[Lumen 2.1|Release Notes]]
** [https://gitlab.com/lumen-code/lumen/-/releases/2.1.0 Version 2.1.0] (06/07/2026)

* Lumen 2.0: [[Lumen 2.0|Release Notes]]
** [https://gitlab.com/lumen-code/lumen/-/releases/2.0.1 Version 2.0.1] (29/01/2026); [https://zenodo.org/records/20121198 Zenodo 20121198]
** [https://gitlab.com/lumen-code/lumen/-/releases/2.0.0 Version 2.0.0] (04/12/2025);

== Previous versions ==

* [https://www.attaccalite.com/lumen/ Lumen 1.X] (2018)

Lumen release notes

2026-07-06T13:26:00Z

Davide: /* Lumen Releases */ tag 2.1.0 created

This page contains the list of the releases of the lumen code.

Starting from version 2.0, the release notes are available in a dedicated page on the wiki, while the source code is available on gitlab.

== Lumen Releases ==

* Lumen 2.1: [[Lumen 2.1|Release Notes]]
** [https://gitlab.com/lumen-code/lumen/-/releases/2.1.0 Tag 2.1.0] (06/07/2026)

* Lumen 2.0: [[Lumen 2.0|Release Notes]]
** [https://gitlab.com/lumen-code/lumen/-/releases/2.0.1 Tag 2.0.1] (29/01/2026); [https://zenodo.org/records/20121198 Zenodo 20121198]
** [https://gitlab.com/lumen-code/lumen/-/releases/2.0.0 Tag 2.0.0] (04/12/2025);

== Previous versions ==

* [https://www.attaccalite.com/lumen/ Lumen 1.X] (2018)

Tutorials

2026-06-24T06:18:33Z

Davide: /* Electron-phonon and exciton-phonon coupling */

== 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 standalone|BSE on h-BN]] (standalone version); [[BSE tutorial on hBN|BSE on h-BN]] (modular version);
* [[How to analyse excitons]]

== Advanced tutorials on GW and Quasi-particles ==
* [[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)

== Advanced tutorials Linear response and BSE ==
=== Optical properties and more in extended systems ===
* [[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)
* [[Magneto-electric reponse]]
* [[Magnons with BSE formalism]] (todo)
* [[Bands selection in BSE and application to semi-core states]] (todo)
* [[Surface spectroscopy]] (to restore from older version)
* [[Fast convergence of dielectric constant with double grid]]
=== Optical properties and more in lower dimensional systems (1D/2D) ===
* [[Hydrogen chain|TDDFT Failure and long range correlations]] (1D)
* [[Si wire|Linear Response from a silicon wire (1D)]] (to restore from old version)
* [[Si Surface|Linear Response of a silicon surface (2D)]]

=== Optical properties and more in molecules ===
* [[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)
* [[Dichroism in molecules]]

== Post Processing with ypp and YamboPy ==
* [[Yambo Post Processing (ypp)]]
* [[First steps in YamboPy]]

== Electron-phonon and exciton-phonon coupling ==
* [[Electron Phonon Coupling (old tutorial)]]
* [[Electron Phonon Coupling (new)]]
* [[Electron and exciton life-time]]
* [[Optical properties at finite temperature]]
* [[Phonon-assisted luminescence by finite atomic displacements]]
* [[Exciton-phonon coupling and luminescence hBN]]
* [[Exciton-phonon coupling and luminescence]]
* [[Thermal lines and special displacements with YamboPy]]

== Real time propagation ==

=== Linear response from real time propagation ===
* [[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 from real time propagation ===

* [[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)]]
* [[Angular dependence of non-linear response]]
* [[Third Harmonic Generation (THG)]]
* [[Spin-orbit coupling and non-linear response]]
* [[Two-photon absorption]]
* [[Sum frequency generation]]
* [[Shift current]]

=== Non-linear response in Floquet formulation ===
* [[SHG from Floquet formulation]]
=== 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]]

=== More on real time propagation ===

* [[Real-time dynamics with external field read from file]]
* [[Parallelization for non-linear response calculations]]
* [[A fast approach to excitonic effects in linear/non-linear response]]
*

== Developing Lumen ==

Main Page

2026-06-24T06:17:52Z

Davide: Removed link to unesed page

Lumen wiki

* [[Installation]]
* [[Tutorials]]
* [https://enccs.github.io/max-coe-workshop/day5-tutorials/ Yambo tutorials on Markdown]
* [[Lumen release notes|Lumen Releases]]
* [https://wiki.yambo-code.eu/wiki/index.php?title=First_steps_in_Yambopy YamboPy]
* [https://www.lumen-code.org/ford/index.html Ford documentation]
* [[Test-suite]]

Internal wiki

2026-06-24T06:15:31Z

Davide: Page made empty, links moved somewhere else

Tutorials

2026-06-24T06:14:37Z

Davide: /* Electron-phonon and exciton-phonon coupling */

== 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 standalone|BSE on h-BN]] (standalone version); [[BSE tutorial on hBN|BSE on h-BN]] (modular version);
* [[How to analyse excitons]]

== Advanced tutorials on GW and Quasi-particles ==
* [[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)

== Advanced tutorials Linear response and BSE ==
=== Optical properties and more in extended systems ===
* [[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)
* [[Magneto-electric reponse]]
* [[Magnons with BSE formalism]] (todo)
* [[Bands selection in BSE and application to semi-core states]] (todo)
* [[Surface spectroscopy]] (to restore from older version)
* [[Fast convergence of dielectric constant with double grid]]
=== Optical properties and more in lower dimensional systems (1D/2D) ===
* [[Hydrogen chain|TDDFT Failure and long range correlations]] (1D)
* [[Si wire|Linear Response from a silicon wire (1D)]] (to restore from old version)
* [[Si Surface|Linear Response of a silicon surface (2D)]]

=== Optical properties and more in molecules ===
* [[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)
* [[Dichroism in molecules]]

== Post Processing with ypp and YamboPy ==
* [[Yambo Post Processing (ypp)]]
* [[First steps in YamboPy]]

== Electron-phonon and exciton-phonon coupling ==
* [[Electron Phonon Coupling (old tutorial)]]
* [[Electron Phonon Coupling (new)]]
* [[Electron_and_exciton_life-time]]
* [[Optical properties at finite temperature]]
* [[Phonon-assisted luminescence by finite atomic displacements]]
* [[Exciton-phonon_coupling_and_luminescence_hBN]]
* [[Exciton-phonon coupling and luminescence]]
* [[Thermal lines and special displacements with YamboPy]]

== Real time propagation ==

=== Linear response from real time propagation ===
* [[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 from real time propagation ===

* [[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)]]
* [[Angular dependence of non-linear response]]
* [[Third Harmonic Generation (THG)]]
* [[Spin-orbit coupling and non-linear response]]
* [[Two-photon absorption]]
* [[Sum frequency generation]]
* [[Shift current]]

=== Non-linear response in Floquet formulation ===
* [[SHG from Floquet formulation]]
=== 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]]

=== More on real time propagation ===

* [[Real-time dynamics with external field read from file]]
* [[Parallelization for non-linear response calculations]]
* [[A fast approach to excitonic effects in linear/non-linear response]]
*

== Developing Lumen ==

Configure library externally

2026-06-21T21:45:08Z

Davide: removed extra newline

== In this page we explain how to configure HDF5 and NetCDFc and NetCDF-Fortran ==

'''Gfortran Installation'''

export INSTALL_DIR={your installation directory}

Configuration of HDF5 1.14.6 with gfortran:

./configure --prefix=${INSTALL_DIR} --enable-fortran --disable-shared --enable-static --without-pic --disable-cxx --enable-parallel CC=mpicc CPP="mpicc -E" CFLAGS="-O2 -std=gnu99 -D_C_US -D_FORTRAN_US" FC=mpif90 --enable-build-mode=production

Configure NetCDFC-4.9.0 with gfortran

./configure --prefix=${INSTALL_DIR} --without-pic --enable-static --disable-shared --disable-dap --enable-netcdf-4 CC=gcc CPP="gcc -E" CFLAGS="-O2 -std=gnu99 -D_C_US -D_FORTRAN_US" FC=mpif90 CPPFLAGS="-DpgiFortran -I${INSTALLDIR}" LIBS="-L${INSTALL_DIR} -lhdf5hl_fortran -lhdf5_fortran -lhdf5_hl -lhdf5" F77=mpif90 LDFLAGS="-Wl,--copy-dt-needed-entries"

Configure NetCDF-fortran 4.6.2 with gfortran

./configure --prefix=${INSTALL_DIR} --libdir="${INSTALL_DIR}/lib" --without-pic --enable-static --disable-shared CC=gcc CFLAGS="-O2 -std=gnu99 -D_C_US -D_FORTRAN_US" CPP="mpif90 -E" FC=mpif90 CPPFLAGS="-DpgiFortran -I${INSTALLDIR}/include" LIBS="-L${INSTALL_DIR}/lib/ -lnetcdf -lhdf5hl_fortran -lhdf5_fortran -lhdf5_hl -hdf5" F77=mpif90

'''Intel ONEAPI Installation'''

export INSTALLDIR={your installation directory}

Configuration of HDF5 1.14.6 with intel:

/configure --prefix=${INSTALL_DIR} --enable-fortran --enable-parallel CC=mpiicx CPP="mpiicx -E" FC=mpiifx --enable-build-mode=production --enable-fortran --disable-shared --enable-static --without-pic --disable-cxx

Configure NetCDFC-4.9.0 with intel

./configure --prefix=${INSTALL_DIR} --without-pic --enable-static --disable-shared --disable-dap --disable-libxml2 --enable-netcdf-4 CC=mpiicx CPP="mpiicx -E" CFLAGS="-O2 -std=gnu99 -D_C_US -D_FORTRAN_US" FC=mpiifort CPPFLAGS="-DpgiFortran -I${INSTALL_DIR}/include" LIBS="-L${INSTALL_DIR}/lib -lhdf5hl_fortran -lhdf5_fortran -lhdf5_hl -lhdf5" F77=mpiifx

Configure NetCDF-fortran 4.6.2 with intel

./configure --prefix=${INSTALL_DIR} --libdir=/home/attacc/local_intel/lib --without-pic --enable-static --disable-shared CC=mpiicx CFLAGS="-O2 -std=gnu99 -D_C_US -D_FORTRAN_US" CPP="mpiicx -E" FC=mpiifx CPPFLAGS="-DpgiFortran -I ${INSTALL_DIR}/include" LIBS="-L${INSTALL_DIR}/lib/ -lnetcdf -lhdf5hl_fortran -lhdf5_fortran -lhdf5_hl -lhdf5" F77=mpiifx

'''Linking libraries to Lumen'''

./configure --with-netcdf-path="${INSTALL_DIR}" --with-netcdff-path="${INSTALL_DIR}" --with-hdf5-path="${INSTALL_DIR}"

Configure library externally

2026-06-21T21:44:35Z

Davide: Style improved

== In this page we explain how to configure HDF5 and NetCDFc and NetCDF-Fortran ==

'''Gfortran Installation'''

export INSTALL_DIR={your installation directory}

Configuration of HDF5 1.14.6 with gfortran:

./configure --prefix=${INSTALL_DIR} --enable-fortran --disable-shared --enable-static --without-pic --disable-cxx --enable-parallel CC=mpicc CPP="mpicc -E" CFLAGS="-O2 -std=gnu99 -D_C_US -D_FORTRAN_US" FC=mpif90 --enable-build-mode=production

Configure NetCDFC-4.9.0 with gfortran

./configure --prefix=${INSTALL_DIR} --without-pic --enable-static --disable-shared --disable-dap --enable-netcdf-4 CC=gcc CPP="gcc -E" CFLAGS="-O2 -std=gnu99 -D_C_US -D_FORTRAN_US" FC=mpif90 CPPFLAGS="-DpgiFortran -I${INSTALLDIR}" LIBS="-L${INSTALL_DIR} -lhdf5hl_fortran -lhdf5_fortran -lhdf5_hl -lhdf5" F77=mpif90 LDFLAGS="-Wl,--copy-dt-needed-entries"

Configure NetCDF-fortran 4.6.2 with gfortran

./configure --prefix=${INSTALL_DIR} --libdir="${INSTALL_DIR}/lib" --without-pic --enable-static --disable-shared CC=gcc CFLAGS="-O2 -std=gnu99 -D_C_US -D_FORTRAN_US" CPP="mpif90 -E" FC=mpif90 CPPFLAGS="-DpgiFortran
-I${INSTALLDIR}/include" LIBS="-L${INSTALL_DIR}/lib/ -lnetcdf -lhdf5hl_fortran -lhdf5_fortran -lhdf5_hl -hdf5" F77=mpif90

'''Intel ONEAPI Installation'''

export INSTALLDIR={your installation directory}

Configuration of HDF5 1.14.6 with intel:

/configure --prefix=${INSTALL_DIR} --enable-fortran --enable-parallel CC=mpiicx CPP="mpiicx -E" FC=mpiifx --enable-build-mode=production --enable-fortran --disable-shared --enable-static --without-pic --disable-cxx

Configure NetCDFC-4.9.0 with intel

./configure --prefix=${INSTALL_DIR} --without-pic --enable-static --disable-shared --disable-dap --disable-libxml2 --enable-netcdf-4 CC=mpiicx CPP="mpiicx -E" CFLAGS="-O2 -std=gnu99 -D_C_US -D_FORTRAN_US" FC=mpiifort CPPFLAGS="-DpgiFortran -I${INSTALL_DIR}/include" LIBS="-L${INSTALL_DIR}/lib -lhdf5hl_fortran -lhdf5_fortran -lhdf5_hl -lhdf5" F77=mpiifx

Configure NetCDF-fortran 4.6.2 with intel

./configure --prefix=${INSTALL_DIR} --libdir=/home/attacc/local_intel/lib --without-pic --enable-static --disable-shared CC=mpiicx CFLAGS="-O2 -std=gnu99 -D_C_US -D_FORTRAN_US" CPP="mpiicx -E" FC=mpiifx CPPFLAGS="-DpgiFortran -I ${INSTALL_DIR}/include" LIBS="-L${INSTALL_DIR}/lib/ -lnetcdf -lhdf5hl_fortran -lhdf5_fortran -lhdf5_hl -lhdf5" F77=mpiifx

'''Linking libraries to Lumen'''

./configure --with-netcdf-path="${INSTALL_DIR}" --with-netcdff-path="${INSTALL_DIR}" --with-hdf5-path="${INSTALL_DIR}"

Installation

2026-06-21T21:42:37Z

Davide: moved here link from main page

In this page we explain how to install Lumen on different machine with different compiler and configurations.

===General installations===
* [[Install Lumen on Ubuntu/LinuxMint]]
* [[Install Lumen on Ubuntu/LinuxMint with Intel compiler]]
* [[Install Lumen on Ubuntu/LinuxMint with NVfortran compiler]]
* [[Install Lumenon MacOS (ARM)]]
* [[Install Lumen on a Debian machine with some libraries provided via modules]]
* [[Install Lumen on MacOsX (from Yuncheng Mao)]]

===Libraries===
* [[Compile the libraries in Lumen independently]]
* [[Configure library externally]] (this seems a duplicate of the previous one, to check)
* [[Install Lumen with external HDF5 and NetCDF libraries]]
* [[Install Lumen with external chase libraries|Install Lumen with external ChASE libraries]]
* [[Install Lumen with external slepc and petsc libraries|Install Lumen with external Slepc and Petsc libraries]]

===Installations on specific machines===
* [[Install Lumen Claudia cluster in Marseille|Install Lumen on Claudia cluster (Marseille]])
* [[Install Lumen on Mesocentre (Marseille)]]
* [[Install Lumen on IRENE machine (TGCC Joliot Curie)]]
* [[Install Lumen on ISM cluster (Rome)]]
* [[Install Lumen on Leonardo|Install Lumen on Leonardo (Cineca, Bologna)]]
* [[Install Lumen on Tirant|Install Lumen/LetzElPhC/Yambopy on Tirant (Univ. of Valencia Cluster)]]
* [[Install Lumen/LetzElPhC/Yambopy on Lynx(Univ. of Valencia)| Install Lumen on Lynx (Univ. of Valencia Cluster)]]
* [[Install Lumen on Viper (MPCDF, Munich)]]

Main Page

2026-06-21T21:41:33Z

Davide: Link moved under installation

Lumen wiki

* [[Installation]]
* [[Tutorials]]<nowiki/>
* [https://enccs.github.io/max-coe-workshop/day5-tutorials/ Yambo tutorials on Markdown]
* [[Lumen release notes|Lumen Releases]]
* [https://wiki.yambo-code.eu/wiki/index.php?title=First_steps_in_Yambopy YamboPy]
* [https://www.lumen-code.org/ford/index.html Ford documentation]
* [[Internal wiki]]
* [[Test-suite]]

Lumen 2.1

2026-06-16T11:26:30Z

Davide:

[[File:Times_logo_red.png|right|thumb|250px|TIMES logo red|Lumen 2.1 release was supported by the EU through the [https://times.uv.es/ TIMES MSCA-DN]]]

==Release Notes==

This is the second release of the lumen code, after the reboot based on yambo 5.3. A Partial sync with branch 5.4 of the yambo code was done.

===Main new features===

* Non-linear response with the Floquet formalism: in Lumen 2.1 it is possible to calculate non-linear response functions using Floquet formalism, for the moment this implementation is limited to the independent particle approximation with or without quasi-particle correction. Tutorial: [[SHG from Floquet formulation]]

* Exciton symmetries: Lumen 2.1 is now able to detect the point group of the system based on the list of symmetry operations. It is then possible to compute exciton symmetries via ypp after solving BSE with diagonalization solver. Tutorials: [[LiF#Detection of excitons symmetries (requires lumen 2.1 or better)| Excitons symmetries in LiF]]; [[How to analyse excitons#Detection of excitons symmetries (requires lumen 2.1 or better)| Excitons symmetries in 2D-hBN]]

* Exciton radiative lifetimes: old implementation of exciton radiative lifetimes in ypp has been recovered.
* Excitonic TR-ARPES: with ypp it is possible to compute TR-ARPES from BSE post-processing
* Exciton-exciton transitions: with ypp it is possible to compute exciton-exciton transitions from BSE post-processing
* Coded symmetry average of the dielectric constant for q=0 momentum
* Added interface to the [https://www.lumen-code.org/wiki/index.php?title=Install_Lumen_with_external_chase_libraries CHASE library] as [[BSE solvers overview|BSE solver]]
* In spin polarized systems, it is now possible to get the contribution of a given spin-channel to absorption
* Added Gaussian smearing for the phonon-assisted response functions
* Added preliminary support for system with non-symmorphic symmetry operations. At present the database conversion can be done, and later the non-symmorphic symmetry operation can be removed via ypp
* Added support to Covariant dipoles in shifted grids with multiple shifts
* Plotting of exciton wave-function: it is now possible to plot both fixed electron and fixed hole case. Moreover it is possible to plot the coefficients of the excitonic WF in k-space. File bxsf is produced.

===Code optimization===
* Memory is now distributed in the old diagonalization solver (without using the ldiago library) interfaced with scalapack
* Fixed memory duplication in Slepc solver
* I/O in BSE can now be controlled by input as none/low/medium/high

===Fixes===

* Fixed compilation on MacOsX
* Fixed reading imaginary part from ndb.QP in yambo_nl
* Different compilation issues fixed: with intel compilers, hdf5, openACC, etc....
* Collisions table: added missing n_sp_pol variable

===More===
* added support to llvm compiler
* upgraded internal version of the following libraries: libxc, slepc, petsc, devxlib
* upgraded interface with libxc to work with version 7.0, and to handle different versions (with different API)
* upgraded interface with slepc/petsc to handle different versions (with different API)
* ExcPh databases now can be loaded also if requested band ranges is smaller than available one
* support for BSE with frozen states improved

BSE solvers overview

2026-06-16T11:23:39Z

Davide: /* Summary */

This tutorial explores the options available in Yambo for the numerical solution of the Bethe-Salpeter equation. Before running this tutorial, you should have carried out the tutorial '''[[Calculating optical spectra including excitonic effects: a step-by-step guide]]'''. We will use the same example system (bulk hBN). The various options correspond to different numerical methods for solving (large) eigenvalue problems. The numerical methods are not covered in this tutorial and the interested user is directed to the references.

After completion of this tutorial, you are expected to recognize which solver to use depending on the system, computational load and desired output. Also, for each of the solver, you should be able to choose the values for the key input parameters.

Three main solvers are available in Yambo:
__TOC__

==Full diagonalization solver (-y d)==
Diagonalization of the full Hamiltonian can be done in serial, using the standard LAPACK library. Also, it is possible to do full diagonalization asking only for a subset of eigenvalues of eigenvectors.

This solver outputs the macroscopic dielectric function over a range of energies chosen by the user. It also calculates all exciton energies and composition in terms of electron-hole pairs. Upon user's choice, the latter is stored on disk for use in postprocessing for exciton analysis and plotting as seen in the tutorial '''[[How to analyse excitons - CECAM 2021 school|How to analyse excitons]]'''.

The computational cost of full diagonalization is cubic in the number of electron-hole pairs so this solver can be used for small system only (few light atoms per unit cell). However, thanks to modern libraries like scalapack and ELPA, it is nowadays possible to fully diagonalize very large BSE matrices exploiting MPI and GPUs<ref name="Milev2025" />.

In Lumen we have two working interfaces with Scalapack and one with ELPA. The modern one, which requires the interface with the [https://gitlab.com/lumen-code/Ldiago ldiago library] developed by M. Nalabothula, is controlled by
BSSmod="d"
The old one, which was recently revised to include memory distribution, but it is still significantly less efficient, is controlled by
BSSmod="o"

An example of a tutorial doing full diagonalization can be found here: '''[[Bethe-Salpeter solver: diagonalization]]'''

==Lanczos-Haydock solver (-y h)==
Subspace iterative [https://en.wikipedia.org/wiki/Lanczos_algorithm | Lanczos algorithm] which by-pass diagonalization with the Haydock approach<ref name="Haydock1980"/>.

This solver outputs the macroscopic dielectric function over a range of energies chosen by the user. It can be applied to medium-large systems as it is substantially cheaper than the full diagonalization (quadratic with the number of electron-hole pairs). However, this approach ''does not provide'' the individual exciton energies and composition in terms of electron-hole pairs.

This solver is controlled by
BSSmod="h"

An example of a tutorial using Haydock can be found here: '''[[Bethe-Salpeter solver: Lanczos-Haydock]]''' 
This solver is also compatible with the double-grid<ref name="Alliati2022" />, see tutorial: [[Fast convergence of dielectric constant with double grid]].

==SLEPc solver (-y s) ==
Subspace iterative algorithm using the [https://slepc.upv.es SLEPc library]<ref name="Hernandez2005" /> (SLEPC solver).

This solver outputs the individual exciton energies and composition about an energy value chosen by the user. It can be applied to medium-large systems but does not provide the overall spectrum. This solver is typically used after obtaining the spectra from the Lanczos-Haydock solver to either analyze a specific peak or getting the energy and composition of a few low-lying excitons. It can be used for large-medium systems. The computational cost grows with the number of excitons one asks to compute.

This solver is controlled by
BSSmod="s"

'''Important:'''
To exploit the slepc solver, you need to link the code with Slepc and Petsc library.
See here for more info: [Install Lumen with external slepc and petsc libraries].

An example of tutorial using Slepc can be found here: '''[[Bethe-Salpeter solver: SLEPC]]'''

==ChASE solver (-y c) ==
This solver is controlled by
BSSmod="c"

'''Important:'''
To exploit the chase solver, you need to link the code with the ChASE library.
See here for more info: [Install Lumen with external chase libraries].

TODO: add a tutorial and few more info

==INVERSION solver (-y i) ==
The inversion solver was designed to invert the Dyson equation frequency by frequency. For each frequency the inversion (which can be done with Lapack) has roughly the same cost of a diagonalization. Needing to do so for every frequency, this solver easily becomes very expensive. It was designed to deal with frequency dependent kernel, and also extended with an old implementation of the double-grid see Ref.<ref name="Kammerland2012" />. However, it is not used much anymore.

This solver is controlled by
BSSmod="i"

==Summary==
A first summary about, the output (dielectric function, exciton energies and composition) and cost of the three solvers are summarized in this table.

[[File:Solver.png|Yambo tutorial image]]

Below some performances on the performance of the solvers able to extract the eigenvectors on two BSE matrices:
* Matrix of size 25920 x 25920 (5,006 GB), Requested eigenvalues: 500, runs on 72 CPU (3 nodesx24 cores), BLACS 8x9
* Matrix of size 61440 x 61440 (28,125 GB). Requested eigenvalues: 500. Runs on 100 CPU (4 nodesx25 cores, 500 GB), BLACS 10x10

Solver : Slepc - Chase - Ldiago(ELPA2) - Ldiago(ELPA1) - Ldiago(heevr) - Odiago(heevx)
25920 : 104.7396s - 127.4608s - 53.8467s - 212.3914s - 240.1606s - 1045.s
61440 : 499.0757s - 497.0590s - 361.0130s - 1269.s - 2134.s - 12660.s

As you can see, the new interface with diagonalization library (Ldiago) using ELPA2 solver on CPU (unfortunately ELPA2 is not GPU ported) is even more efficient that iterative solvers like Slepc and Chase. However, it has a larger memory footprint (peak memory is about 6 times the BSE matrix).
The timing of the Haydock solver is not reported, but it remains the fastest method.

More info about the performances of different solvers can be found in this Merge Request:
[https://gitlab.com/lumen-code/lumen/-/merge_requests/304#note_3161882534 Performances of BSE solvers] 

==Navigate==
* Back to [[BSE tutorial on hBN]] module
* Back to [[Tutorials]] main page

=References =
<references>
<ref name="Haydock1980">R. Haydock, ''Solid State Phys.'', '''35''' 215 (1980) edited by H. Ehrenfest, F. Seitz, and D. Turnbull, Academic Press</ref>
<ref name="Hernandez2005">V. Hernandez, J.E. Roman and V. Vidal in '' ACM Transactions on Mathematical Software'', '''31''' 315 (2005)</ref>
<ref name="Kammerland2012" >David Kammerlander, Silvana Botti, Miguel A. L Marques, Andrea Marini, and Claudio Attaccalite, Phys. Rev. B 86, 125203 (2012) </ref>
<ref name="Alliati2022">I.M. Alliati, D. Sangalli, M. Gruning, ''Double k-Grid Method for Solving the Bethe-Salpeter Equation via Lanczos Approaches'', [https://www.frontiersin.org/journals/chemistry/articles/10.3389/fchem.2021.763946/full Frontiers in Chemistry, '''9''', p.763946 (2022)]</ref>
<ref name="Milev2025">P. Milev, et al. TO FIX</ref>
</references>

BSE solvers overview

2026-06-16T11:21:43Z

Davide: /* Full diagonalization solver (-y d) */

This tutorial explores the options available in Yambo for the numerical solution of the Bethe-Salpeter equation. Before running this tutorial, you should have carried out the tutorial '''[[Calculating optical spectra including excitonic effects: a step-by-step guide]]'''. We will use the same example system (bulk hBN). The various options correspond to different numerical methods for solving (large) eigenvalue problems. The numerical methods are not covered in this tutorial and the interested user is directed to the references.

After completion of this tutorial, you are expected to recognize which solver to use depending on the system, computational load and desired output. Also, for each of the solver, you should be able to choose the values for the key input parameters.

Three main solvers are available in Yambo:
__TOC__

==Full diagonalization solver (-y d)==
Diagonalization of the full Hamiltonian can be done in serial, using the standard LAPACK library. Also, it is possible to do full diagonalization asking only for a subset of eigenvalues of eigenvectors.

This solver outputs the macroscopic dielectric function over a range of energies chosen by the user. It also calculates all exciton energies and composition in terms of electron-hole pairs. Upon user's choice, the latter is stored on disk for use in postprocessing for exciton analysis and plotting as seen in the tutorial '''[[How to analyse excitons - CECAM 2021 school|How to analyse excitons]]'''.

The computational cost of full diagonalization is cubic in the number of electron-hole pairs so this solver can be used for small system only (few light atoms per unit cell). However, thanks to modern libraries like scalapack and ELPA, it is nowadays possible to fully diagonalize very large BSE matrices exploiting MPI and GPUs<ref name="Milev2025" />.

In Lumen we have two working interfaces with Scalapack and one with ELPA. The modern one, which requires the interface with the [https://gitlab.com/lumen-code/Ldiago ldiago library] developed by M. Nalabothula, is controlled by
BSSmod="d"
The old one, which was recently revised to include memory distribution, but it is still significantly less efficient, is controlled by
BSSmod="o"

An example of a tutorial doing full diagonalization can be found here: '''[[Bethe-Salpeter solver: diagonalization]]'''

==Lanczos-Haydock solver (-y h)==
Subspace iterative [https://en.wikipedia.org/wiki/Lanczos_algorithm | Lanczos algorithm] which by-pass diagonalization with the Haydock approach<ref name="Haydock1980"/>.

This solver outputs the macroscopic dielectric function over a range of energies chosen by the user. It can be applied to medium-large systems as it is substantially cheaper than the full diagonalization (quadratic with the number of electron-hole pairs). However, this approach ''does not provide'' the individual exciton energies and composition in terms of electron-hole pairs.

This solver is controlled by
BSSmod="h"

An example of a tutorial using Haydock can be found here: '''[[Bethe-Salpeter solver: Lanczos-Haydock]]''' 
This solver is also compatible with the double-grid<ref name="Alliati2022" />, see tutorial: [[Fast convergence of dielectric constant with double grid]].

==SLEPc solver (-y s) ==
Subspace iterative algorithm using the [https://slepc.upv.es SLEPc library]<ref name="Hernandez2005" /> (SLEPC solver).

This solver outputs the individual exciton energies and composition about an energy value chosen by the user. It can be applied to medium-large systems but does not provide the overall spectrum. This solver is typically used after obtaining the spectra from the Lanczos-Haydock solver to either analyze a specific peak or getting the energy and composition of a few low-lying excitons. It can be used for large-medium systems. The computational cost grows with the number of excitons one asks to compute.

This solver is controlled by
BSSmod="s"

'''Important:'''
To exploit the slepc solver, you need to link the code with Slepc and Petsc library.
See here for more info: [Install Lumen with external slepc and petsc libraries].

An example of tutorial using Slepc can be found here: '''[[Bethe-Salpeter solver: SLEPC]]'''

==ChASE solver (-y c) ==
This solver is controlled by
BSSmod="c"

'''Important:'''
To exploit the chase solver, you need to link the code with the ChASE library.
See here for more info: [Install Lumen with external chase libraries].

TODO: add a tutorial and few more info

==INVERSION solver (-y i) ==
The inversion solver was designed to invert the Dyson equation frequency by frequency. For each frequency the inversion (which can be done with Lapack) has roughly the same cost of a diagonalization. Needing to do so for every frequency, this solver easily becomes very expensive. It was designed to deal with frequency dependent kernel, and also extended with an old implementation of the double-grid see Ref.<ref name="Kammerland2012" />. However, it is not used much anymore.

This solver is controlled by
BSSmod="i"

==Summary==
A first summary about, the output (dielectric function, exciton energies and composition) and cost of the three solvers are summarized in this table.

[[File:Solver.png|Yambo tutorial image]]

Also some performances on the performance of solvers able to extract the eigenvectors on a matrix of size 61440 x 61440 (28,125 GB). Requested eigenvalues: 500. Runs on 100 CPU (4 nodesx25 cores, 500 GB), BLACS 10x10

Solver : Slepc - Chase - Ldiago(ELPA2) - Ldiago(ELPA1) - Ldiago(heevr) - Odiago(heevx)
25920 : 104.7396s - 127.4608s - 53.8467s - 212.3914s - 240.1606s - 1045.s
61440 : 499.0757s - 497.0590s - 361.0130s - 1269.s - 2134.s - 12660.s

As you can see, the new interface with diagonalization library (Ldiago) using ELPA2 solver on CPU (unfortunately ELPA2 is not GPU ported) is even more efficient that iterative solvers like Slepc and Chase. However, it has a larger memory footprint (peak memory is about 6 times the BSE matrix).
The timing of the Haydock solver is not reported, but it remains the fastest method.

More info about the performances of different solvers can be found in this Merge Request:
[https://gitlab.com/lumen-code/lumen/-/merge_requests/304#note_3161882534 Performances of BSE solvers] 

==Navigate==
* Back to [[BSE tutorial on hBN]] module
* Back to [[Tutorials]] main page

=References =
<references>
<ref name="Haydock1980">R. Haydock, ''Solid State Phys.'', '''35''' 215 (1980) edited by H. Ehrenfest, F. Seitz, and D. Turnbull, Academic Press</ref>
<ref name="Hernandez2005">V. Hernandez, J.E. Roman and V. Vidal in '' ACM Transactions on Mathematical Software'', '''31''' 315 (2005)</ref>
<ref name="Kammerland2012" >David Kammerlander, Silvana Botti, Miguel A. L Marques, Andrea Marini, and Claudio Attaccalite, Phys. Rev. B 86, 125203 (2012) </ref>
<ref name="Alliati2022">I.M. Alliati, D. Sangalli, M. Gruning, ''Double k-Grid Method for Solving the Bethe-Salpeter Equation via Lanczos Approaches'', [https://www.frontiersin.org/journals/chemistry/articles/10.3389/fchem.2021.763946/full Frontiers in Chemistry, '''9''', p.763946 (2022)]</ref>
<ref name="Milev2025">P. Milev, et al. TO FIX</ref>
</references>

BSE solvers overview

2026-06-16T11:20:05Z

Davide: /* Summary */

This tutorial explores the options available in Yambo for the numerical solution of the Bethe-Salpeter equation. Before running this tutorial, you should have carried out the tutorial '''[[Calculating optical spectra including excitonic effects: a step-by-step guide]]'''. We will use the same example system (bulk hBN). The various options correspond to different numerical methods for solving (large) eigenvalue problems. The numerical methods are not covered in this tutorial and the interested user is directed to the references.

After completion of this tutorial, you are expected to recognize which solver to use depending on the system, computational load and desired output. Also, for each of the solver, you should be able to choose the values for the key input parameters.

Three main solvers are available in Yambo:
__TOC__

==Full diagonalization solver (-y d)==
Diagonalization of the full Hamiltonian can be done in serial, using the standard LAPACK library. Also, it is possible to do full diagonalization asking only for a subset of eigenvalues of eigenvectors.

This solver outputs the macroscopic dielectric function over a range of energies chosen by the user. It also calculates all exciton energies and composition in terms of electron-hole pairs. Upon user's choice, the latter is stored on disk for use in postprocessing for exciton analysis and plotting as seen in the tutorial '''[[How to analyse excitons - CECAM 2021 school|How to analyse excitons]]'''.

The computational cost of full diagonalization is cubic in the number of electron-hole pairs so this solver can be used for small system only (few light atoms per unit cell). However, thanks to modern libraries like scalapack and ELPA, it is nowadays possible to fully diagonalize very large BSE matrices exploiting MPI and GPUs<ref name="Milev2025" />. In Lumen we have two working interfaces with Scalapack and one with ELPA.

The modern one, which requires the interface with the ldiago library, is controlled by
BSSmod="d"
The old one, which was recently revised to include memory distribution, but it is still significantly less efficient, is controlled by
BSSmod="o"

An example of a tutorial doing full diagonalization can be found here: '''[[Bethe-Salpeter solver: diagonalization]]'''

==Lanczos-Haydock solver (-y h)==
Subspace iterative [https://en.wikipedia.org/wiki/Lanczos_algorithm | Lanczos algorithm] which by-pass diagonalization with the Haydock approach<ref name="Haydock1980"/>.

This solver outputs the macroscopic dielectric function over a range of energies chosen by the user. It can be applied to medium-large systems as it is substantially cheaper than the full diagonalization (quadratic with the number of electron-hole pairs). However, this approach ''does not provide'' the individual exciton energies and composition in terms of electron-hole pairs.

This solver is controlled by
BSSmod="h"

An example of a tutorial using Haydock can be found here: '''[[Bethe-Salpeter solver: Lanczos-Haydock]]''' 
This solver is also compatible with the double-grid<ref name="Alliati2022" />, see tutorial: [[Fast convergence of dielectric constant with double grid]].

==SLEPc solver (-y s) ==
Subspace iterative algorithm using the [https://slepc.upv.es SLEPc library]<ref name="Hernandez2005" /> (SLEPC solver).

This solver outputs the individual exciton energies and composition about an energy value chosen by the user. It can be applied to medium-large systems but does not provide the overall spectrum. This solver is typically used after obtaining the spectra from the Lanczos-Haydock solver to either analyze a specific peak or getting the energy and composition of a few low-lying excitons. It can be used for large-medium systems. The computational cost grows with the number of excitons one asks to compute.

This solver is controlled by
BSSmod="s"

'''Important:'''
To exploit the slepc solver, you need to link the code with Slepc and Petsc library.
See here for more info: [Install Lumen with external slepc and petsc libraries].

An example of tutorial using Slepc can be found here: '''[[Bethe-Salpeter solver: SLEPC]]'''

==ChASE solver (-y c) ==
This solver is controlled by
BSSmod="c"

'''Important:'''
To exploit the chase solver, you need to link the code with the ChASE library.
See here for more info: [Install Lumen with external chase libraries].

TODO: add a tutorial and few more info

==INVERSION solver (-y i) ==
The inversion solver was designed to invert the Dyson equation frequency by frequency. For each frequency the inversion (which can be done with Lapack) has roughly the same cost of a diagonalization. Needing to do so for every frequency, this solver easily becomes very expensive. It was designed to deal with frequency dependent kernel, and also extended with an old implementation of the double-grid see Ref.<ref name="Kammerland2012" />. However, it is not used much anymore.

This solver is controlled by
BSSmod="i"

==Summary==
A first summary about, the output (dielectric function, exciton energies and composition) and cost of the three solvers are summarized in this table.

[[File:Solver.png|Yambo tutorial image]]

Also some performances on the performance of solvers able to extract the eigenvectors on a matrix of size 61440 x 61440 (28,125 GB). Requested eigenvalues: 500. Runs on 100 CPU (4 nodesx25 cores, 500 GB), BLACS 10x10

Solver : Slepc - Chase - Ldiago(ELPA2) - Ldiago(ELPA1) - Ldiago(heevr) - Odiago(heevx)
25920 : 104.7396s - 127.4608s - 53.8467s - 212.3914s - 240.1606s - 1045.s
61440 : 499.0757s - 497.0590s - 361.0130s - 1269.s - 2134.s - 12660.s

As you can see, the new interface with diagonalization library (Ldiago) using ELPA2 solver on CPU (unfortunately ELPA2 is not GPU ported) is even more efficient that iterative solvers like Slepc and Chase. However, it has a larger memory footprint (peak memory is about 6 times the BSE matrix).
The timing of the Haydock solver is not reported, but it remains the fastest method.

More info about the performances of different solvers can be found in this Merge Request:
[https://gitlab.com/lumen-code/lumen/-/merge_requests/304#note_3161882534 Performances of BSE solvers] 

==Navigate==
* Back to [[BSE tutorial on hBN]] module
* Back to [[Tutorials]] main page

=References =
<references>
<ref name="Haydock1980">R. Haydock, ''Solid State Phys.'', '''35''' 215 (1980) edited by H. Ehrenfest, F. Seitz, and D. Turnbull, Academic Press</ref>
<ref name="Hernandez2005">V. Hernandez, J.E. Roman and V. Vidal in '' ACM Transactions on Mathematical Software'', '''31''' 315 (2005)</ref>
<ref name="Kammerland2012" >David Kammerlander, Silvana Botti, Miguel A. L Marques, Andrea Marini, and Claudio Attaccalite, Phys. Rev. B 86, 125203 (2012) </ref>
<ref name="Alliati2022">I.M. Alliati, D. Sangalli, M. Gruning, ''Double k-Grid Method for Solving the Bethe-Salpeter Equation via Lanczos Approaches'', [https://www.frontiersin.org/journals/chemistry/articles/10.3389/fchem.2021.763946/full Frontiers in Chemistry, '''9''', p.763946 (2022)]</ref>
<ref name="Milev2025">P. Milev, et al. TO FIX</ref>
</references>

BSE solvers overview

2026-06-16T10:44:35Z

Davide: /* Summary */

This tutorial explores the options available in Yambo for the numerical solution of the Bethe-Salpeter equation. Before running this tutorial, you should have carried out the tutorial '''[[Calculating optical spectra including excitonic effects: a step-by-step guide]]'''. We will use the same example system (bulk hBN). The various options correspond to different numerical methods for solving (large) eigenvalue problems. The numerical methods are not covered in this tutorial and the interested user is directed to the references.

After completion of this tutorial, you are expected to recognize which solver to use depending on the system, computational load and desired output. Also, for each of the solver, you should be able to choose the values for the key input parameters.

Three main solvers are available in Yambo:
__TOC__

==Full diagonalization solver (-y d)==
Diagonalization of the full Hamiltonian can be done in serial, using the standard LAPACK library. Also, it is possible to do full diagonalization asking only for a subset of eigenvalues of eigenvectors.

This solver outputs the macroscopic dielectric function over a range of energies chosen by the user. It also calculates all exciton energies and composition in terms of electron-hole pairs. Upon user's choice, the latter is stored on disk for use in postprocessing for exciton analysis and plotting as seen in the tutorial '''[[How to analyse excitons - CECAM 2021 school|How to analyse excitons]]'''.

The computational cost of full diagonalization is cubic in the number of electron-hole pairs so this solver can be used for small system only (few light atoms per unit cell). However, thanks to modern libraries like scalapack and ELPA, it is nowadays possible to fully diagonalize very large BSE matrices exploiting MPI and GPUs<ref name="Milev2025" />. In Lumen we have two working interfaces with Scalapack and one with ELPA.

The modern one, which requires the interface with the ldiago library, is controlled by
BSSmod="d"
The old one, which was recently revised to include memory distribution, but it is still significantly less efficient, is controlled by
BSSmod="o"

An example of a tutorial doing full diagonalization can be found here: '''[[Bethe-Salpeter solver: diagonalization]]'''

==Lanczos-Haydock solver (-y h)==
Subspace iterative [https://en.wikipedia.org/wiki/Lanczos_algorithm | Lanczos algorithm] which by-pass diagonalization with the Haydock approach<ref name="Haydock1980"/>.

This solver outputs the macroscopic dielectric function over a range of energies chosen by the user. It can be applied to medium-large systems as it is substantially cheaper than the full diagonalization (quadratic with the number of electron-hole pairs). However, this approach ''does not provide'' the individual exciton energies and composition in terms of electron-hole pairs.

This solver is controlled by
BSSmod="h"

An example of a tutorial using Haydock can be found here: '''[[Bethe-Salpeter solver: Lanczos-Haydock]]''' 
This solver is also compatible with the double-grid<ref name="Alliati2022" />, see tutorial: [[Fast convergence of dielectric constant with double grid]].

==SLEPc solver (-y s) ==
Subspace iterative algorithm using the [https://slepc.upv.es SLEPc library]<ref name="Hernandez2005" /> (SLEPC solver).

This solver outputs the individual exciton energies and composition about an energy value chosen by the user. It can be applied to medium-large systems but does not provide the overall spectrum. This solver is typically used after obtaining the spectra from the Lanczos-Haydock solver to either analyze a specific peak or getting the energy and composition of a few low-lying excitons. It can be used for large-medium systems. The computational cost grows with the number of excitons one asks to compute.

This solver is controlled by
BSSmod="s"

'''Important:'''
To exploit the slepc solver, you need to link the code with Slepc and Petsc library.
See here for more info: [Install Lumen with external slepc and petsc libraries].

An example of tutorial using Slepc can be found here: '''[[Bethe-Salpeter solver: SLEPC]]'''

==ChASE solver (-y c) ==
This solver is controlled by
BSSmod="c"

'''Important:'''
To exploit the chase solver, you need to link the code with the ChASE library.
See here for more info: [Install Lumen with external chase libraries].

TODO: add a tutorial and few more info

==INVERSION solver (-y i) ==
The inversion solver was designed to invert the Dyson equation frequency by frequency. For each frequency the inversion (which can be done with Lapack) has roughly the same cost of a diagonalization. Needing to do so for every frequency, this solver easily becomes very expensive. It was designed to deal with frequency dependent kernel, and also extended with an old implementation of the double-grid see Ref.<ref name="Kammerland2012" />. However, it is not used much anymore.

This solver is controlled by
BSSmod="i"

==Summary==
A first summary about, the output (dielectric function, exciton energies and composition) and cost of the three solvers are summarized in this table.

[[File:Solver.png|Yambo tutorial image]]

Here more info about the performances of different solvers:
[https://gitlab.com/lumen-code/lumen/-/merge_requests/304#note_3161882534 Performances of BSE solvers] 
As you can see, the ELPA2 solver on CPU is even more efficient that iterative solvers like Slepc and Chase

==Navigate==
* Back to [[BSE tutorial on hBN]] module
* Back to [[Tutorials]] main page

=References =
<references>
<ref name="Haydock1980">R. Haydock, ''Solid State Phys.'', '''35''' 215 (1980) edited by H. Ehrenfest, F. Seitz, and D. Turnbull, Academic Press</ref>
<ref name="Hernandez2005">V. Hernandez, J.E. Roman and V. Vidal in '' ACM Transactions on Mathematical Software'', '''31''' 315 (2005)</ref>
<ref name="Kammerland2012" >David Kammerlander, Silvana Botti, Miguel A. L Marques, Andrea Marini, and Claudio Attaccalite, Phys. Rev. B 86, 125203 (2012) </ref>
<ref name="Alliati2022">I.M. Alliati, D. Sangalli, M. Gruning, ''Double k-Grid Method for Solving the Bethe-Salpeter Equation via Lanczos Approaches'', [https://www.frontiersin.org/journals/chemistry/articles/10.3389/fchem.2021.763946/full Frontiers in Chemistry, '''9''', p.763946 (2022)]</ref>
<ref name="Milev2025">P. Milev, et al. TO FIX</ref>
</references>

BSE solvers overview

2026-06-16T10:43:36Z

Davide:

This tutorial explores the options available in Yambo for the numerical solution of the Bethe-Salpeter equation. Before running this tutorial, you should have carried out the tutorial '''[[Calculating optical spectra including excitonic effects: a step-by-step guide]]'''. We will use the same example system (bulk hBN). The various options correspond to different numerical methods for solving (large) eigenvalue problems. The numerical methods are not covered in this tutorial and the interested user is directed to the references.

After completion of this tutorial, you are expected to recognize which solver to use depending on the system, computational load and desired output. Also, for each of the solver, you should be able to choose the values for the key input parameters.

Three main solvers are available in Yambo:
__TOC__

==Full diagonalization solver (-y d)==
Diagonalization of the full Hamiltonian can be done in serial, using the standard LAPACK library. Also, it is possible to do full diagonalization asking only for a subset of eigenvalues of eigenvectors.

This solver outputs the macroscopic dielectric function over a range of energies chosen by the user. It also calculates all exciton energies and composition in terms of electron-hole pairs. Upon user's choice, the latter is stored on disk for use in postprocessing for exciton analysis and plotting as seen in the tutorial '''[[How to analyse excitons - CECAM 2021 school|How to analyse excitons]]'''.

The computational cost of full diagonalization is cubic in the number of electron-hole pairs so this solver can be used for small system only (few light atoms per unit cell). However, thanks to modern libraries like scalapack and ELPA, it is nowadays possible to fully diagonalize very large BSE matrices exploiting MPI and GPUs<ref name="Milev2025" />. In Lumen we have two working interfaces with Scalapack and one with ELPA.

The modern one, which requires the interface with the ldiago library, is controlled by
BSSmod="d"
The old one, which was recently revised to include memory distribution, but it is still significantly less efficient, is controlled by
BSSmod="o"

An example of a tutorial doing full diagonalization can be found here: '''[[Bethe-Salpeter solver: diagonalization]]'''

==Lanczos-Haydock solver (-y h)==
Subspace iterative [https://en.wikipedia.org/wiki/Lanczos_algorithm | Lanczos algorithm] which by-pass diagonalization with the Haydock approach<ref name="Haydock1980"/>.

This solver outputs the macroscopic dielectric function over a range of energies chosen by the user. It can be applied to medium-large systems as it is substantially cheaper than the full diagonalization (quadratic with the number of electron-hole pairs). However, this approach ''does not provide'' the individual exciton energies and composition in terms of electron-hole pairs.

This solver is controlled by
BSSmod="h"

An example of a tutorial using Haydock can be found here: '''[[Bethe-Salpeter solver: Lanczos-Haydock]]''' 
This solver is also compatible with the double-grid<ref name="Alliati2022" />, see tutorial: [[Fast convergence of dielectric constant with double grid]].

==SLEPc solver (-y s) ==
Subspace iterative algorithm using the [https://slepc.upv.es SLEPc library]<ref name="Hernandez2005" /> (SLEPC solver).

This solver outputs the individual exciton energies and composition about an energy value chosen by the user. It can be applied to medium-large systems but does not provide the overall spectrum. This solver is typically used after obtaining the spectra from the Lanczos-Haydock solver to either analyze a specific peak or getting the energy and composition of a few low-lying excitons. It can be used for large-medium systems. The computational cost grows with the number of excitons one asks to compute.

This solver is controlled by
BSSmod="s"

'''Important:'''
To exploit the slepc solver, you need to link the code with Slepc and Petsc library.
See here for more info: [Install Lumen with external slepc and petsc libraries].

An example of tutorial using Slepc can be found here: '''[[Bethe-Salpeter solver: SLEPC]]'''

==ChASE solver (-y c) ==
This solver is controlled by
BSSmod="c"

'''Important:'''
To exploit the chase solver, you need to link the code with the ChASE library.
See here for more info: [Install Lumen with external chase libraries].

TODO: add a tutorial and few more info

==INVERSION solver (-y i) ==
The inversion solver was designed to invert the Dyson equation frequency by frequency. For each frequency the inversion (which can be done with Lapack) has roughly the same cost of a diagonalization. Needing to do so for every frequency, this solver easily becomes very expensive. It was designed to deal with frequency dependent kernel, and also extended with an old implementation of the double-grid see Ref.<ref name="Kammerland2012" />. However, it is not used much anymore.

This solver is controlled by
BSSmod="i"

==Summary==
A first summary about, the output (dielectric function, exciton energies and composition) and cost of the three solvers are summarized in this table.

[[File:Solver.png|Yambo tutorial image]]

Here more info about the performances of different solvers:
[https://gitlab.com/lumen-code/lumen/-/merge_requests/304#note_3161882534 Performances of BSE solvers]

==Navigate==
* Back to [[BSE tutorial on hBN]] module
* Back to [[Tutorials]] main page

=References =
<references>
<ref name="Haydock1980">R. Haydock, ''Solid State Phys.'', '''35''' 215 (1980) edited by H. Ehrenfest, F. Seitz, and D. Turnbull, Academic Press</ref>
<ref name="Hernandez2005">V. Hernandez, J.E. Roman and V. Vidal in '' ACM Transactions on Mathematical Software'', '''31''' 315 (2005)</ref>
<ref name="Kammerland2012" >David Kammerlander, Silvana Botti, Miguel A. L Marques, Andrea Marini, and Claudio Attaccalite, Phys. Rev. B 86, 125203 (2012) </ref>
<ref name="Alliati2022">I.M. Alliati, D. Sangalli, M. Gruning, ''Double k-Grid Method for Solving the Bethe-Salpeter Equation via Lanczos Approaches'', [https://www.frontiersin.org/journals/chemistry/articles/10.3389/fchem.2021.763946/full Frontiers in Chemistry, '''9''', p.763946 (2022)]</ref>
<ref name="Milev2025">P. Milev, et al. TO FIX</ref>
</references>

BSE solvers overview

2026-06-16T10:35:30Z

Davide: Improving info on solvers

This tutorial explores the options available in Yambo for the numerical solution of the Bethe-Salpeter equation. Before running this tutorial, you should have carried out the tutorial '''[[Calculating optical spectra including excitonic effects: a step-by-step guide]]'''. We will use the same example system (bulk hBN). The various options correspond to different numerical methods for solving (large) eigenvalue problems. The numerical methods are not covered in this tutorial and the interested user is directed to the references.

After completion of this tutorial, you are expected to recognize which solver to use depending on the system, computational load and desired output. Also, for each of the solver, you should be able to choose the values for the key input parameters.

Three main solvers are available in Yambo:
__TOC__

==Full diagonalization solver (-y d)==
Diagonalization of the full Hamiltonian can be done in serial, using the standard LAPACK library. Also, it is possible to do full diagonalization asking only for a subset of eigenvalues of eigenvectors.

This solver outputs the macroscopic dielectric function over a range of energies chosen by the user. It also calculates all exciton energies and composition in terms of electron-hole pairs. Upon user's choice, the latter is stored on disk for use in postprocessing for exciton analysis and plotting as seen in the tutorial '''[[How to analyse excitons - CECAM 2021 school|How to analyse excitons]]'''.

The computational cost of full diagonalization is cubic in the number of electron-hole pairs so this solver can be used for small system only (few light atoms per unit cell). However, thanks to modern libraries like scalapack and ELPA, it is nowadays possible to fully diagonalize very large BSE matrices exploiting MPI and GPUs<ref name="Milev2025" />. In Lumen we have two working interfaces with Scalapack and one with ELPA.

An example of a tutorial doing full diagonalization can be found here: '''[[Bethe-Salpeter solver: diagonalization]]'''

==Lanczos-Haydock solver (-y h)==
Subspace iterative [https://en.wikipedia.org/wiki/Lanczos_algorithm | Lanczos algorithm] which by-pass diagonalization with the Haydock approach<ref name="Haydock1980"/>.

This solver outputs the macroscopic dielectric function over a range of energies chosen by the user. It can be applied to medium-large systems as it is substantially cheaper than the full diagonalization (quadratic with the number of electron-hole pairs). However, this approach ''does not provide'' the individual exciton energies and composition in terms of electron-hole pairs.

An example of a tutorial using Haydock can be found here: '''[[Bethe-Salpeter solver: Lanczos-Haydock]]''' 
This solver is also compatible with the double-grid<ref name="Alliati2022" />, see tutorial: [[Fast convergence of dielectric constant with double grid]].

==SLEPc solver (-y s) ==
Subspace iterative algorithm using the [https://slepc.upv.es SLEPc library]<ref name="Hernandez2005" /> (SLEPC solver).

This solver outputs the individual exciton energies and composition about an energy value chosen by the user. It can be applied to medium-large systems but does not provide the overall spectrum. This solver is typically used after obtaining the spectra from the Lanczos-Haydock solver to either analyze a specific peak or getting the energy and composition of a few low-lying excitons. It can be used for large-medium systems. The computational cost grows with the number of excitons one asks to compute.

'''Important:'''
To exploit the slepc solver, you need to link the code with Slepc and Petsc library.
See here for more infoL [Install Lumen with external slepc and petsc libraries].

An example of tutorial using Slepc can be found here: '''[[Bethe-Salpeter solver: SLEPC]]'''

==ChASE solver (-y c) ==
TO ADD INFO

==INVERSION solver (-y i) ==
The inversion solver was designed to invert the Dyson equation frequency by frequency. For each frequency the inversion (which can be done with Lapack) has roughly the same cost of a diagonalization. Needing to do so for every frequency, this solver easily becomes very expensive. It was designed to deal with frequency dependent kernel, and also extended with an old implementation of the double-grid see Ref.<ref name="Kammerland2012" />. However, it is not used much anymore.

==Summary==
To summarise, the output (dielectric function, exciton energies and composition) and cost of the three solvers are summarized in this table.

[[File:Solver.png|Yambo tutorial image]]

==Navigate==
* Back to [[BSE tutorial on hBN]] module
* Back to [[Tutorials]] main page

=References =
<references>
<ref name="Haydock1980">R. Haydock, ''Solid State Phys.'', '''35''' 215 (1980) edited by H. Ehrenfest, F. Seitz, and D. Turnbull, Academic Press</ref>
<ref name="Hernandez2005">V. Hernandez, J.E. Roman and V. Vidal in '' ACM Transactions on Mathematical Software'', '''31''' 315 (2005)</ref>
<ref name="Kammerland2012" >David Kammerlander, Silvana Botti, Miguel A. L Marques, Andrea Marini, and Claudio Attaccalite, Phys. Rev. B 86, 125203 (2012) </ref>
<ref name="Alliati2022">I.M. Alliati, D. Sangalli, M. Gruning, ''Double k-Grid Method for Solving the Bethe-Salpeter Equation via Lanczos Approaches'', [https://www.frontiersin.org/journals/chemistry/articles/10.3389/fchem.2021.763946/full Frontiers in Chemistry, '''9''', p.763946 (2022)]</ref>
<ref name="Milev2025">P. Milev, et al. TO FIX</ref>
</references>

Install Lumen with external slepc and petsc libraries

2026-06-16T10:32:36Z

Davide:

You need available Slepc and Petsc libraries compiled (eventually with gpu support).

To get the Petsc and Slepc libraries see here: TODO

Alternatively, you can try the automatic compilation of the internal Petsc and Slepc library by adding to the configure script:
--enable-slepc-linalg

...

Install Lumen with external slepc and petsc libraries

2026-06-16T10:30:26Z

Davide: First version of the page created

Install Lumen with external chase libraries

2026-06-16T10:28:38Z

Davide:

Available since Lumen 2.1.
You need an available Chase library compiled (eventually with gpu support).
To get the chase library see here: https://github.com/ChASE-library/ChASE

For tests and discussion on the chase interface with lumen see here: https://gitlab.com/lumen-code/lumen/-/merge_requests/304

== Compile chase with nvfortran ==
module load nvidia_profile

cmake \
-DCMAKE_CXX_COMPILER=nvc++ \
-DCMAKE_C_COMPILER=nvc \
-DCMAKE_Fortran_COMPILER=nvfortran \
-DCMAKE_BUILD_TYPE=RelWithDebInfo \
-DCMAKE_CXX_STANDARD_REQUIRED=ON \
-DCMAKE_CXX_STANDARD=17 \
-DCMAKE_CUDA_ARCHITECTURES=120 \
-DMPI_C_COMPILER=mpicc \
-DMPI_CXX_COMPILER=mpicxx \
-DMPI_Fortran_COMPILER=mpifort \
-DSCALAPACK_DIR=/opt/nvidia/hpc_sdk/Linux_x86_64/2025/comm_libs/13.0/hpcx/hpcx-2.25.1/ompi/lib -DSCALAPACK_VERBOSE=ON \
-DCMAKE_INSTALL_PREFIX="SET_PATH_HERE" \
../../

== Compile chase with gfortran ==
module load gnu_profile

cmake \
-DCMAKE_CXX_COMPILER=mpicxx \
-DCMAKE_C_COMPILER=mpicc \
-DCMAKE_Fortran_COMPILER=mpif90 \
-DCMAKE_BUILD_TYPE=RelWithDebInfo \
-DMPI_C_COMPILER=mpicc \
-DMPI_CXX_COMPILER=mpicxx \
-DMPI_Fortran_COMPILER=mpif90 \
-DSCALAPACK_DIR=/opt/gcc/15.2.0/scalapack-2.2.2/ -DSCALAPACK_VERBOSE=ON \
-DCMAKE_INSTALL_PREFIX=/home/sangalli/data/local-libs/ChASE-github/gfortran15.2_openmpi5.0/install \
../../

== Link chase ==
--enable-chase \
--with-chase-libdir="${CHASE_LIBS}/lib" \
--with-chase-includedir="${CHASE_LIBS}/include" \

Installation

2026-06-16T10:28:06Z

Davide:

In this page we explain how to install Lumen on different machine with different compiler and configurations.

===General installations===
* [[Install Lumen on Ubuntu/LinuxMint]]
* [[Install Lumen on Ubuntu/LinuxMint with Intel compiler]]
* [[Install Lumen on Ubuntu/LinuxMint with NVfortran compiler]]
* [[Install Lumenon MacOS (ARM)]]
* [[Install Lumen on a Debian machine with some libraries provided via modules]]
* [[Install Lumen on MacOsX (from Yuncheng Mao)]]

===Libraries===
* [[Compile the libraries in Lumen independently]]
* [[Install Lumen with external HDF5 and NetCDF libraries]]
* [[Install Lumen with external chase libraries|Install Lumen with external ChASE libraries]]
* [[Install Lumen with external slepc and petsc libraries|Install Lumen with external Slepc and Petsc libraries]]

===Installations on specific machines===
* [[Install Lumen Claudia cluster in Marseille|Install Lumen on Claudia cluster (Marseille]])
* [[Install Lumen on Mesocentre (Marseille)]]
* [[Install Lumen on IRENE machine (TGCC Joliot Curie)]]
* [[Install Lumen on ISM cluster (Rome)]]
* [[Install Lumen on Leonardo|Install Lumen on Leonardo (Cineca, Bologna)]]
* [[Install Lumen on Tirant|Install Lumen/LetzElPhC/Yambopy on Tirant (Univ. of Valencia Cluster)]]
* [[Install Lumen/LetzElPhC/Yambopy on Lynx(Univ. of Valencia)| Install Lumen on Lynx (Univ. of Valencia Cluster)]]
* [[Install Lumen on Viper (MPCDF, Munich)]]

BSE solvers overview

2026-06-16T08:54:02Z

Davide: updating solvers info

Lumen 2.1

2026-06-16T08:51:00Z

Davide: Added link to Exciton symmetries tutorial

[[File:Times_logo_red.png|right|thumb|250px|TIMES logo red|Lumen 2.1 release was supported by the EU through the [https://times.uv.es/ TIMES MSCA-DN]]]

==Release Notes==

This is the second release of the lumen code, after the reboot based on yambo 5.3. A Partial sync with branch 5.4 of the yambo code was done.

===Main new features===

* Non-linear response with the Floquet formalism: in Lumen 2.1 it is possible to calculate non-linear response functions using Floquet formalism, for the moment this implementation is limited to the independent particle approximation with or without quasi-particle correction. Tutorial: [[SHG from Floquet formulation]]

* Exciton symmetries: Lumen 2.1 is now able to detect the point group of the system based on the list of symmetry operations. It is then possible to compute exciton symmetries via ypp after solving BSE with diagonalization solver. Tutorials: [[LiF#Detection_of_excitons_symmetries_(requires_lumen_2.1_or_better)| Excitons symmetries in LiF]]; [[How_to_analyse_excitons#Detection_of_excitons_symmetries_(requires_lumen_2.1_or_better)| Excitons symmetries in 2D-hBN]]

* Exciton radiative lifetimes: old implementation of exciton radiative lifetimes in ypp has been recovered.
* Excitonic TR-ARPES: with ypp it is possible to compute TR-ARPES from BSE post-processing
* Exciton-exciton transitions: with ypp it is possible to compute exciton-exciton transitions from BSE post-processing
* Coded symmetry average of the dielectric constant for q=0 momentum
* Added interface to the CHASE library as BSE solver
* In spin polarized systems, it is now possible to get the contribution of a given spin-channel to absorption
* Added Gaussian smearing for the phonon-assisted response functions
* Added preliminary support for system with non-symmorphic symmetry operations. At present the database conversion can be done, and later the non-symmorphic symmetry operation can be removed via ypp
* Added support to Covariant dipoles in shifted grids with multiple shifts
* Plotting of exciton wave-function: it is now possible to plot both fixed electron and fixed hole case. Moreover it is possible to plot the coefficients of the excitonic WF in k-space. File bxsf is produced.

===Code optimization===
* Memory is now distributed in the old diagonalization solver (without using the ldiago library) interfaced with scalapack
* Fixed memory duplication in Slepc solver
* I/O in BSE can now be controlled by input as none/low/medium/high

===Fixes===

* Fixed compilation on MacOsX
* Fixed reading imaginary part from ndb.QP in yambo_nl
* Different compilation issues fixed: with intel compilers, hdf5, openACC, etc....
* Collisions table: added missing n_sp_pol variable

===More===
* added support to llvm compiler
* upgraded internal version of the following libraries: libxc, slepc, petsc, devxlib
* upgraded interface with libxc to work with version 7.0, and to handle different versions (with different API)
* upgraded interface with slepc/petsc to handle different versions (with different API)
* ExcPh databases now can be loaded also if requested band ranges is smaller than available one
* support for BSE with frozen states improved

How to analyse excitons

2026-06-16T08:49:32Z

Davide:

In this tutorial you will learn (for a 2D-hBN) how to:
* analyze a BSE optical spectrum in terms of excitonic eigenvectors and eigenvalues
* look at the spatial distribution of the exciton
* For a similar tutorial focusing on 3D-hBN, go to [[How to analyse excitons - CECAM 2021 school|this page]].

== Prerequisites ==
'''Previous modules'''
* You must have completed the [[How to treat low dimensional systems]] tutorial
'''You will need''':
* <code>ypp </code> executable
* <code>xcrysden</code> executable
* <code>gnuplot or xmgrace</code> executable
'''

== YAMBO calculations ==
If you have completed the tutorials of 2D hBN you should have all the databases required to do this tutorial in your <code>YAMBO_TUTORIALS/hBN-2D/SAVE</code> and <code>2D_WR_WC</code> (databases generated with RIM and cutoff) directories

$ ls ./SAVE
ndb.gops ndb.kindx ns.db1 ns.kb_pp_pwscf_fragment_1 ....
$ ls ./2D_WR_WC
ndb.BS_Q1_CPU_0 ndb.cutoff ndb.dip_iR_and_P_fragment_1 ndb.pp_fragment_1 ...

==Sort the excitonic eigenvalues==

$ ypp -J 2D_WR_WC -e s -b 1

We are sorting the excitons for the q-index = 1 (optical limit q=0).
The new generated file ''o-2D_WR_WC.exc_qpt1_E_sorted '' (''o-2D_WR_WC.exc_qpt1_I_sorted '') reports
the energies of the excitons and their Dipole Oscillator Strengths sorted by energy (Index).
[[File:strengh.png|none|x120px|]]
Open the first file and look inside. The first exciton is at 4.83 eV and the second one has the highest strength (normalized to 1)

Or you can make a plot
$ gnuplot
gnuplot> set style line 2 lc rgb 'black' pt 7 # circle
gnuplot> plot 'o-2D_WR_WC.exc_qpt1_E_sorted' with points ls 2 title 'Strenghts'

[[File:strenght.png|none|600px|]]

Attention the convergence of these results with different k-points grids is mandatory!

== Calculate the exciton oscillator strenght and amplitude ==

We can now analyze the excitons in terms of single-particle states, to do that create the appropriate input
$ ypp -F ypp_AMPL.in -J 2D_WR_WC -e a 1

Suppose you wish to analyze the first 5 excitons then change this line as:
States= "1 - 5" # Index of the BS state(s)

Close the input and run ypp

$ ypp -F ypp_AMPL.in -J 2D_WR_WC

$ls o*exc*at*
o-2D_WR_WC.exc_qpt1_amplitude_at_1 o-2D_WR_WC.exc_qpt1_weights_at_1 ...

For an exciton <math>|\lambda></math> , ''o-2D_WR_WC.exc_qpt1_weights_at_*'' report the Weights
[[File:Weights.png|none|x60px|]]
and ''o-2D_WR_WC.exc_qpt1_amplitude_**'' report the amplitudes
[[File:Ampl.png|none|x70px|]]

Open the file ''o-2D_WR_WC.exc_weights_at_1''
# Band_V Band_C K ibz Symm. Weight Energy
#
4.000000 5.000000 7.000000 2.000000 0.922095 4.401093
4.000000 5.000000 7.000000 1.000000 0.922086 4.401093

The first exciton is essentially done of only single particle transitions from VBM to CBM at K (last k-point of the grid).

'''Note: the transistion Energy reported in the weight file are Kohn-Sham energy difference and not Quasiparticle energies.'''

[[File:Amplitude_plot.png|none|600px|]]

== Plot the exciton spatial distribution ==

To see the spatial character of the exciton YPP writes the exciton spatial distribution, in other words the probability to find the electron somewhere in the space when the hole is fixed in a give position.
Different output formats can be selected and 1D,2D,3D plots done.
Create the input and change the size of the cell where to see the exciton.
Note that If the k-grid of the BSE simulation is a NxNx1 the exciton has an induced fictitious periodicity every Nx Nx1 Cell of the simulation.
For hBN-2D this is not a problem because the exciton is strongly localized but in other systems with more delocalized excitons to look at the real exciton size it is necessary to use
very large k-grids in the BSE
$ ypp -F ypp_WF.in -J 2D_WR_WC -e w 1

excitons # [R] Excitons
wavefunction # [R] Wavefunction
Format= "x" # Output format [(c)ube/(g)nuplot/(x)crysden]
Direction= "12" # [rlu] [1/2/3] for 1d or [12/13/23] for 2d [123] for 3D
FFTGvecs= 3951 RL # [FFT] Plane-waves
States= "1 - 1" # Index of the BS state(s)
Degen_Step= 0.0100 eV # Maximum energy separation of two degenerate states
% Cells
5 | 5 | 1 | # Number of cell repetitions in each direction (odd or 1)
%
% Hole
2.4 | 1.400 | 0.00 | # [cc] Hole position in unit cell

Close the input and run ypp
$ ypp -F ypp_WF.in -J 2D_WR_WC

$ xcrysden --xsf o-2D_WR_WC.exc_2d_1.xsf

[[File:exc_BN2D.png|none|400px|]]



== Plot electron/hole average density==

Another way to analyze excitons, it is the possibility to plot the average electron/hole densities defined as:

[[File:Electron hole density.png|400px]]

to generate the corresponding input just type

ypp -F ypp_WF.in -e w -avehole

and choose the exciton you want to plot. The electron/hole average densities correspond to generalized valence/conduction orbitals for a given exciton. They are interesting in particular for molecular crystals because they allow to distinguishing charge-transfer exctions from Frenkel or Wannier ones, by looking the relative position of the electron/hole densities.
For example see the electron/hole average density calculated in a Metal Organic Framework (MOF) that contains azobenzene:

[[File:Average electron hole.png|center |600px|Yambo tutorial image]]

this figure is taken from ref. <ref>Strongly Bound Excitons in Metal-Organic Framework MOF-5: A Many-Body Perturbation Theory Study, A. R. Kshirsagar et al., [https://doi.org/10.26434/chemrxiv.14034917.v1 preprint ChemRxiv]</ref>.

== Interpolate exciton dispersion (only in Yambo 5.x)==

Starting from Yambo 5.x it is possible to calculate excitons at finite momentum '''q'''. The momentum can be specified with the variable
%BSEQptR
iq_start | iq_end | # [BSK] Transferred momenta range
%
where iq_start and iq_end is the first and last index of momentum in the irreducible Brillouin zone (IBZ). You can find the list of q-vectors in the r_setup report.
If you calculate the BSE for all q-points of the IBZ, then you can interpolate exciton dispersion along any direction of the full Brillouin zone(BZ) using the
interpolation scheme of ref. <ref>Warren E. Pickett, Henry Krakauer, and Philip B. Allen [https://journals.aps.org/prb/abstract/10.1103/PhysRevB.38.2721 PRB 38 p2721 (1988)]</ref>, by doing:

ypp -e i

and then setting:

excitons # [R] Excitonic properties
interpolate # [R] Interpolate
States= "1 - 4" # Index of the BS state(s)
INTERP_mode= "BOLTZ" # Interpolation mode (NN=nearest point, BOLTZ=boltztrap aproach)
% INTERP_Grid
-1 |-1 |-1 | # Interpolation BZ Grid
%
#PrtDOS # Print Exciton Density of States
BANDS_steps= 100 # Number of divisions
cooIn= "rlu" # Points coordinates (in) cc/rlu/iku/alat
cooOut= "rlu" # Points coordinates (out) cc/rlu/iku/alat
%BANDS_kpts # K points of the bands circuit
0 | 0 | 0 |
0.333333333333333 | 0.33333333333333 | 0 |
%

Running ypp, it will interpolate the first 4 excitons energies on 100 points along with the line Gamma->K. If you plot the output file "o.excitons_interpolated",
you can compare the interpolated results with the original points calculated in the BSE, that are available in the report "r_excitons_interpolate", and you get something like:

[[File:Excitons dispersion.png|center | Yambo exciton dispersion]]

It is also possible to interpolate excitons on a given regular grid by setting <code>INTERP_Grid</code>, for example in the figure below we interpolate the dispersion of the lowest exciton in MoS2 calculated on a 39x39x1 and interpolated in a 117x117x1 grid in all the Brillouin zone:

[[File:Screenshot from 2021-03-24 11-38-46.png|center| 700px|Yambo tutorial image]]

Finally if you turn on the flag <code>PrtDOS</code> the code calculate the excitonic density of states on the excitons interpolated on the regular grid given by <code>INTERP_Grid</code>.

= Detection of excitons symmetries (requires lumen 2.1 or better) =

Here we show how to detect excitons symmetries in hBN<ref name="Nalabothula2026"/> from the BSE eigenvectors.

TO DO ...

== References ==

<references>
<ref name="Nalabothula2026"> M. Nalabothula, D. Sangalli, F. Paleari, S. Reichardt, L. Wirtz, [https://doi.org/10.1103/7h2w-ys9h Phys. Rev. B. '''113''', 205130 (2026)] </ref>
</references>

How to analyse excitons

2026-06-16T08:49:05Z

Davide:

In this tutorial you will learn (for a 2D-hBN) how to:
* analyze a BSE optical spectrum in terms of excitonic eigenvectors and eigenvalues
* look at the spatial distribution of the exciton
* For a similar tutorial focusing on 3D-hBN, go to [[How to analyse excitons - CECAM 2021 school|this page]].

== Prerequisites ==
'''Previous modules'''
* You must have completed the [[How to treat low dimensional systems]] tutorial
'''You will need''':
* <code>ypp </code> executable
* <code>xcrysden</code> executable
* <code>gnuplot or xmgrace</code> executable
'''

== YAMBO calculations ==
If you have completed the tutorials of 2D hBN you should have all the databases required to do this tutorial in your <code>YAMBO_TUTORIALS/hBN-2D/SAVE</code> and <code>2D_WR_WC</code> (databases generated with RIM and cutoff) directories

$ ls ./SAVE
ndb.gops ndb.kindx ns.db1 ns.kb_pp_pwscf_fragment_1 ....
$ ls ./2D_WR_WC
ndb.BS_Q1_CPU_0 ndb.cutoff ndb.dip_iR_and_P_fragment_1 ndb.pp_fragment_1 ...

==Sort the excitonic eigenvalues==

$ ypp -J 2D_WR_WC -e s -b 1

We are sorting the excitons for the q-index = 1 (optical limit q=0).
The new generated file ''o-2D_WR_WC.exc_qpt1_E_sorted '' (''o-2D_WR_WC.exc_qpt1_I_sorted '') reports
the energies of the excitons and their Dipole Oscillator Strengths sorted by energy (Index).
[[File:strengh.png|none|x120px|]]
Open the first file and look inside. The first exciton is at 4.83 eV and the second one has the highest strength (normalized to 1)

Or you can make a plot
$ gnuplot
gnuplot> set style line 2 lc rgb 'black' pt 7 # circle
gnuplot> plot 'o-2D_WR_WC.exc_qpt1_E_sorted' with points ls 2 title 'Strenghts'

[[File:strenght.png|none|600px|]]

Attention the convergence of these results with different k-points grids is mandatory!

== Calculate the exciton oscillator strenght and amplitude ==

We can now analyze the excitons in terms of single-particle states, to do that create the appropriate input
$ ypp -F ypp_AMPL.in -J 2D_WR_WC -e a 1

Suppose you wish to analyze the first 5 excitons then change this line as:
States= "1 - 5" # Index of the BS state(s)

Close the input and run ypp

$ ypp -F ypp_AMPL.in -J 2D_WR_WC

$ls o*exc*at*
o-2D_WR_WC.exc_qpt1_amplitude_at_1 o-2D_WR_WC.exc_qpt1_weights_at_1 ...

For an exciton <math>|\lambda></math> , ''o-2D_WR_WC.exc_qpt1_weights_at_*'' report the Weights
[[File:Weights.png|none|x60px|]]
and ''o-2D_WR_WC.exc_qpt1_amplitude_**'' report the amplitudes
[[File:Ampl.png|none|x70px|]]

Open the file ''o-2D_WR_WC.exc_weights_at_1''
# Band_V Band_C K ibz Symm. Weight Energy
#
4.000000 5.000000 7.000000 2.000000 0.922095 4.401093
4.000000 5.000000 7.000000 1.000000 0.922086 4.401093

The first exciton is essentially done of only single particle transitions from VBM to CBM at K (last k-point of the grid).

'''Note: the transistion Energy reported in the weight file are Kohn-Sham energy difference and not Quasiparticle energies.'''

[[File:Amplitude_plot.png|none|600px|]]

== Plot the exciton spatial distribution ==

To see the spatial character of the exciton YPP writes the exciton spatial distribution, in other words the probability to find the electron somewhere in the space when the hole is fixed in a give position.
Different output formats can be selected and 1D,2D,3D plots done.
Create the input and change the size of the cell where to see the exciton.
Note that If the k-grid of the BSE simulation is a NxNx1 the exciton has an induced fictitious periodicity every Nx Nx1 Cell of the simulation.
For hBN-2D this is not a problem because the exciton is strongly localized but in other systems with more delocalized excitons to look at the real exciton size it is necessary to use
very large k-grids in the BSE
$ ypp -F ypp_WF.in -J 2D_WR_WC -e w 1

excitons # [R] Excitons
wavefunction # [R] Wavefunction
Format= "x" # Output format [(c)ube/(g)nuplot/(x)crysden]
Direction= "12" # [rlu] [1/2/3] for 1d or [12/13/23] for 2d [123] for 3D
FFTGvecs= 3951 RL # [FFT] Plane-waves
States= "1 - 1" # Index of the BS state(s)
Degen_Step= 0.0100 eV # Maximum energy separation of two degenerate states
% Cells
5 | 5 | 1 | # Number of cell repetitions in each direction (odd or 1)
%
% Hole
2.4 | 1.400 | 0.00 | # [cc] Hole position in unit cell

Close the input and run ypp
$ ypp -F ypp_WF.in -J 2D_WR_WC

$ xcrysden --xsf o-2D_WR_WC.exc_2d_1.xsf

[[File:exc_BN2D.png|none|400px|]]



== Plot electron/hole average density (only in Yambo 5.x)==

Another way to analyze excitons, it is the possibility to plot the average electron/hole densities defined as:

[[File:Electron hole density.png|400px]]

to generate the corresponding input just type

ypp -F ypp_WF.in -e w -avehole

and choose the exciton you want to plot. The electron/hole average densities correspond to generalized valence/conduction orbitals for a given exciton. They are interesting in particular for molecular crystals because they allow to distinguishing charge-transfer exctions from Frenkel or Wannier ones, by looking the relative position of the electron/hole densities.
For example see the electron/hole average density calculated in a Metal Organic Framework (MOF) that contains azobenzene:

[[File:Average electron hole.png|center |600px|Yambo tutorial image]]

this figure is taken from ref. <ref>Strongly Bound Excitons in Metal-Organic Framework MOF-5: A Many-Body Perturbation Theory Study, A. R. Kshirsagar et al., [https://doi.org/10.26434/chemrxiv.14034917.v1 preprint ChemRxiv]</ref>.

== Interpolate exciton dispersion (only in Yambo 5.x)==

Starting from Yambo 5.x it is possible to calculate excitons at finite momentum '''q'''. The momentum can be specified with the variable
%BSEQptR
iq_start | iq_end | # [BSK] Transferred momenta range
%
where iq_start and iq_end is the first and last index of momentum in the irreducible Brillouin zone (IBZ). You can find the list of q-vectors in the r_setup report.
If you calculate the BSE for all q-points of the IBZ, then you can interpolate exciton dispersion along any direction of the full Brillouin zone(BZ) using the
interpolation scheme of ref. <ref>Warren E. Pickett, Henry Krakauer, and Philip B. Allen [https://journals.aps.org/prb/abstract/10.1103/PhysRevB.38.2721 PRB 38 p2721 (1988)]</ref>, by doing:

ypp -e i

and then setting:

excitons # [R] Excitonic properties
interpolate # [R] Interpolate
States= "1 - 4" # Index of the BS state(s)
INTERP_mode= "BOLTZ" # Interpolation mode (NN=nearest point, BOLTZ=boltztrap aproach)
% INTERP_Grid
-1 |-1 |-1 | # Interpolation BZ Grid
%
#PrtDOS # Print Exciton Density of States
BANDS_steps= 100 # Number of divisions
cooIn= "rlu" # Points coordinates (in) cc/rlu/iku/alat
cooOut= "rlu" # Points coordinates (out) cc/rlu/iku/alat
%BANDS_kpts # K points of the bands circuit
0 | 0 | 0 |
0.333333333333333 | 0.33333333333333 | 0 |
%

Running ypp, it will interpolate the first 4 excitons energies on 100 points along with the line Gamma->K. If you plot the output file "o.excitons_interpolated",
you can compare the interpolated results with the original points calculated in the BSE, that are available in the report "r_excitons_interpolate", and you get something like:

[[File:Excitons dispersion.png|center | Yambo exciton dispersion]]

It is also possible to interpolate excitons on a given regular grid by setting <code>INTERP_Grid</code>, for example in the figure below we interpolate the dispersion of the lowest exciton in MoS2 calculated on a 39x39x1 and interpolated in a 117x117x1 grid in all the Brillouin zone:

[[File:Screenshot from 2021-03-24 11-38-46.png|center| 700px|Yambo tutorial image]]

Finally if you turn on the flag <code>PrtDOS</code> the code calculate the excitonic density of states on the excitons interpolated on the regular grid given by <code>INTERP_Grid</code>.

= Detection of excitons symmetries (requires lumen 2.1 or better) =

Here we show how to detect excitons symmetries in hBN<ref name="Nalabothula2026"/> from the BSE eigenvectors.

TO DO ...

== References ==

<references>
<ref name="Nalabothula2026"> M. Nalabothula, D. Sangalli, F. Paleari, S. Reichardt, L. Wirtz, [https://doi.org/10.1103/7h2w-ys9h Phys. Rev. B. '''113''', 205130 (2026)] </ref>
</references>

LiF

2026-06-16T08:42:52Z

Davide:

= Prerequisites =
== The material ==

[[File:LiF-crystal.jpg|x300px|frameless|LiF crystal]]
[[File:LiF_lif_structure.jpg|x300px|frameless|LiF structure]]
[[File:lif_rholf.png|x350px|frameless|Optical absorption of solid LiF. Dots: experiment. Continuous line: BSE. Dashed line: RPA-QP]]

LiF is an isulator with a wide band-ga and the following characteristics:
* FCC lattice
* Two atoms per cell, Li and F (8 electrons)
* Lattice constant 7.61 [a.u.]
* Plane waves cutoff 40 Hartree (1800 RL vectors)

In this tutorial we will use LiF to show the importance of excitonic effects in the optical properites of materials with weak screening.
We will start the TDDFT tutorial exactly where standard local approximations for fxc (like ALDA) do not work at all. The optical spectrum of solid
[[LiF]] <ref name="wikipedia/> is dominated by a strongly bound exciton (3 eV binding energy) <ref name="Rohlfing"/><ref name="Marini"/><ref name="Adragna"/><ref name="Sottile"/><ref name="Reining"/><ref name="Botti"/><ref name="Sangalli"/> and will immediately see that any independent particle approximation drastically fails to describe the complex structures that are observed experimentally.

Follow the instructions in [[Tutorials_download]] and download/unpack the [https://media.yambo-code.eu/educational/tutorials/files/LiF.tar.gz LiF.tar.gz]. 
'''Important''': Increase the number of bands in the NSCF DFT up to 50

== Computing the ground state properties with PWscf ==

$ cd YAMBO_TUTORIALS/LiF/PWSCF
$ ls
inputs outputs psps

First run the SCF calculation to generate the ground-state charge density, occupations, Fermi level, and so on:
$ pw.x < inputs/scf.in > scf.out

Next run a non-SCF calculation to generate a set of Kohn-Sham eigenvalues and eigenvectors for both occupied and unoccupied states (100 bands):
$ pw.x < inputs/nscf.in > nscf.out

Note the presence of the following flags in the input file:
wf_collect=.true.
force_symmorphic=.true.
diago_full_acc=.true.
which are needed for generating the Yambo databases accurately. Full explanations of these variables are given on the [http://www.quantum-espresso.org/wp-content/uploads/Doc/INPUT_PW.html quantum-ESPRESSO input variables page].

After these two runs, you should have a ''LiF.save'' directory:
$ ls LiF.save
charge-density.dat data-file-schema.xml F.LDA.cpi.UPF Li.LDA.cpi.UPF wfc10.dat wfc1.dat wfc2.dat wfc3.dat wfc4.dat wfc5.dat wfc6.dat wfc7.dat wfc8.dat wfc9.dat

== Conversion to Yambo format ==
The PWscf ''LiF.save'' output is converted to the Yambo format using the <code>p2y</code> executable (pwscf to yambo), found in the yambo ''bin'' directory.
Enter ''LiF.save'' and launch <code>p2y</code>:

$ cd LiF.save
$ p2y
...
<---> DBs path set to .
<---> Index file set to data-file.xml
<---> Header/K-points/Energies... done
...
<---> == DB1 (Gvecs and more) ...
<---> ... Database done
<---> == DB2 (wavefunctions) ... done ==
<---> == DB3 (PseudoPotential) ... done ==
<---> == P2Y completed ==

This output repeats some information about the system and generates a ''SAVE'' directory:
$ ls SAVE
ns.db1 ns.wf ns.kb_pp_pwscf
ns.wf_fragments_1_1 ...
ns.kb_pp_pwscf_fragment_1 ...
These files, with an ''n'' prefix, indicate that they are in netCDF format, and thus not human readable. However, they are perfectly transferable across different architectures.

In practice we suggest to move the ''SAVE'' folder into a new clean folder.

In this tutorial however, we ask instead that you continue using a ''SAVE'' folder that we prepared previously:
$ cd ../../Optics/YAMBO
$ ls
SAVE

In the following we will describe in detail the different TDDFT related calculations. For each example you can find a link to the Input file with a brief explanation of the variables meaning.

We will also use this first system to explain in detail some of the most relevant technical aspects of yambo.

= Initialization (yambo -i -V RL) =

This section describes the standard initialization procedure. As this initialization is common to all yambo runs it will be described extensively only in this case.

To run this example enter the ''YAMBO'' directory (where you have previously created the SAVE folder) and type:

<pre> localhost:>yambo -F Inputs/01_init -J 01_init
</pre>
and you will see
<pre> <---> [01] Game Summary
<---> [02] Input variables setup
<---> RL-shells | | [000%] --(E) --(X)
<---> RL-shells |####################| [100%] --(E) --(X)
<---> [02.01] K-grid lattice
<---> [02.02] Input (E)nergies[ev] & Occupations
<---> [03] Transferred momenta grid: Indexes
<---> [04] Game Over & Game summary
</pre>
If now you edit the ''r-01_setup'' you will find a lot of useless ([[File:LiF_smile.gif|60x23px]]) informations about your system like, the Fermi level

<pre> [02.02] Input (E)nergies[ev] & Occupations
========

Fermi Energy[ev] - T[ev/K] :-0.241627 0.000100 1.160400
Bands summary : Full Empty
0001-0004 0005-0050
Indirect Gap [ev][min-max]: 9.11330 14.79631</pre>
confirming that, yes, our bulk LiF is a '''wide''' gap insulator. Then you can see the effect of using
<pre> MaxGvecs= 800 RL # [INI] Max number of G-vectors planned to use
</pre>
in the ''Inputs/01_init'' file. In the report file we notice that
<pre> G-vectors [RL space]: 1807
[wavefunctions]: 1807</pre>
whereas the Reciprocal Space shells found by yambo are
<pre> G-vector (S)hells. Format: [Snn] G`s
[S0030]: 773 [S0029]: 749 [S0028]: 725 [S0027]: 701 [S0026]: 645
[S0025]: 609 [S0024]: 561 [S0023]: 537 [S0022]: 531 [S0021]: 459
[S0020]: 411 [S0019]: 387 [S0018]: 339 [S0017]: 331 [S0016]: 307
[S0015]: 283 [S0014]: 259 [S0013]: 229 [S0012]: 181 [S0011]: 169
[S0010]: 137 [S0009]: 113 [S0008]: 89 [S0007]: 65 [S0006]: 59
[S0005]: 51 [S0004]: 27 [S0003]: 15 [S0002]: 9 [S0001]: 1</pre>
so that the highest number of RL vectors we can use is 773, lower than the charge number given by Abinit and different from the value for ''MaxGvecs'' we provided. This happens because yambo has calculated the Reciprocal Space shells up to 800 to find that the highest closed shell corresponds to 773 RL vectors. Note that, independently of the value you set for any variable that defines a dimension in the RL space, yambo redefines it, in order to match it to the nearest closed shell.



=Random-Phase approximation (yambo -o c) (yambo -o c -k hartree) (yambo -o c -k hartree -V qp) =

The simplest approximation that can be used to calculate the absorption spectrum of LiF is the independent particle approximation. This is done in ''Inputs/02_RPA_no_LF'', that is generated directly by "yambo -o c" command line without specifying an approximation for the kernel (in fact there is no kernel!), To run this example type

<pre> localhost:>yambo -F Inputs/02_RPA_no_LF -J 02_RPA_no_LF
</pre>
The optional ''-J'' flag is used to label the output/report/log files.
In ''Inputs/03_RPA_LF'' we include a kernel in the Hartree approximation. This correspond to the Random-Phase-Approximation (RPA). The response function size is changed to 51 RL. This means including Local Field Effects corresponding to charge oscillations expanded up to the 51st RL component. The converged value for ''NGsBlkXd'' shall be determined by doing several calculations with different values and checking the effect on the final physical result (the absorption spectrum, in this case).

After having run

<pre> localhost:>yambo -F Inputs/03_RPA_LF -J 03_RPA_LF
</pre>
we can compare the result of our calculations with the experimental curve (that you find in the ''Experiment'' folder).

[[File:LiF_02_03.png|512x384px
]]

The comparison with the experiment clearly show that the RPA is not at all adequate in the case of LiF. We notice two main discrepancies between theory and experiment: the RPA absorption onset is too low and the shape lacks the sharp peak observed experimentally at around 12 eV. The onset position can be artificially corrected using a scissor operator (representing the QP gap correction) of 5.19 eV that simply opens rigidly the LDA gap. This is shown in ''Inputs/03_RPA_LF_QP'' with the line

<pre>% XfnQP_E
5.190000 | 1.000000 | 1.000000 | # [EXTQP Xd] E parameters (c/v)
%</pre>
The resulting spectrum fits better with the experiment but now there seems to be an even more serious problem:

[[File:LiF_03.png|512x384px]]

there is an experimental peak well below the QP absorption onset:

an [http://en.wikipedia.org/wiki/Exciton EXCITON!]

[[File:LiF_exciton.jpg|160x202px]]


= ALDA (yambo -o c -k alda) =

Our first attempt to go beyond RPA is using TDDFT in the Adiabatic LDA approximation. Depending on the exchange-correlation(xc) functional used in the ground state calculation yambo will produce a corresponding xc-kernel. In this case, as reported by the ''a2y'' command we are using a Perdew Wang parametrization.
Generate input file with the command: <code> yambo -o c -k alda -F tdlda.in</code>

optics # [R] Linear Response optical properties
chi # [R][CHI] Dyson equation for Chi.
tddft # [R][K] Use TDDFT kernel
#DegFix # Force the code to impose the energy levels to respect their degeneracy
Chimod= "ALDA" # [X] IP/Hartree/ALDA/LRC/PF/BSfxc
FxcGRLc= 51 RL # [TDDFT] XC-kernel RL size
NGsBlkXd= 51 RL # [Xd] Response block size
% QpntsRXd
 1 | 1 | # [Xd] Transferred momenta
%
% BndsRnXd
 1 | 10 | # [Xd] Polarization function bands
%
% EnRngeXd
 7.50000 | 25.00000 | eV # [Xd] Energy range
%
% DmRngeXd
0.100000 | 0.100000 | eV # [Xd] Damping range
%
ETStpsXd= 300 # [Xd] Total Energy steps
% LongDrXd
1.000000 | 0.000000 | 0.000000 | # [Xd] [cc] Electric Field
%

In this input we set to 51 RL (reciprocal lattice vectors) the size of the dielectric constant. This is also the number of RL vectors
used for the exchange correlation functional (FxcGRLc).
You can run the input with the command: <code>yambo -F tdlda.in</code>

Output will be in the file <code>o.eps_q1_inv_alda_dyson</code>, where the first column is the energy while the second and third are the imaginary
and real parts of the dielectric constant.

[[File:LiF_04.png|512x384px]]

Unfortunately ALDA only slightly changes the RPA result. The situation will change in more isolated system (reduced dimensionality).



=The Statically screened Electron-electron interaction (yambo -X s) =

As simple RPA and ALDA approximations have not described correctly the experimental spectrum we have to move towards more elaborate techniques like the Bethe-Salpeter (BS) equation.

A key ingredient in the BS kernel is the electron-electron interaction commonly evaluated in the static approximation. The input file ''Inputs/05_W'' describes how to calculate it. The variables used in this input file have the same physical meaning of those used in the optical absorption calculation. The only difference is that, in general, the response function dimension obtained in the examples (03-04), gives an upper bound to the number of RL vectors needed here. This is because the size of response matrix in an RPA calculation defines also the size of the Hartree potential, whose short-range components are not screened. In the present case, instead, the electron-electron interaction is screened and, for this reason, the RL vectors needed are considerably smaller than in the RPA case. So let's type (no ''-J'' option here) ...

<pre> localhost:>yambo -F Inputs/05_W </pre>
... and after waiting for some minutes the result, that is the database SAVE/ndb.em1s. Do not use the ''-J'' option here as we will need to read this database in the next examples.



=The Bethe-Salpeter equation, Excitons (yambo -o b -k sex -y h -V qp) =

The input file of BSE describes how to calculate an excitonic absorption spectrum, by solving the BS equation using the recursive Haydock Method.
The input file can be generated with the command above. Before running ''yambo'' we must have set the paramteter of the input file, in red the one we changed:

optics # [R] Linear Response optical properties
bss # [R] BSE solver
bse # [R][BSE] Bethe Salpeter Equation.
dipoles # [R] Oscillator strenghts (or dipoles)
#DegFix # Force the code to impose the energy levels to respect their degeneracy
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
BSEdiskIO= "medium" # [BSE] medium (restart possible only in corr.)|low (writes at the end)|none
BSSmod= "h" # [BSS] (h)aydock/(d)iagonalization/(s)lepc/(i)nversion/(t)ddft`
BSENGexx= 773 RL # [BSK] Exchange components
BSENGBlk= 51 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_DbGd_INTERP_mode= "NN" # [EXTQP BSK BSS] Interpolation DbGd mode
% KfnQP_E
 5.800000 | 1.000000 | 1.000000 | # [EXTQP BSK BSS] E parameters (c/v) eV|adim|adim
%
% BSEQptR
1 | 1 | # [BSK] Transferred momenta range
%
% BSEBands
 2 | 7 | # [BSK] Bands range
%
% BEnRange
 0.00000 | 20.00000 | eV # [BSS] Energy range
%
%BDmRange
0.100000 | 0.100000 | eV # [BSS] Damping range
%
BEnSteps= 100 # [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
BSHayTrs=-0.020000 # [BSS] Relative [o/o] Haydock threshold. Strict(>0)/Average(<0)

The line
<pre>% KfnQP_E
5.80000 | 1.000000 | 1.000000 | # [EXTQP BSK BSS] E parameters (c/v)
%</pre>
applies a 5.8 eV scissor in order to open the LDA gap.
<pre>BSKmod= "SEX" # [BSE] IP/Hartree/HF/ALDA/SEX</pre>
Here we are using the screened exchange (SEX) approximation for the kernel that includes both exchange and correlation terms.
The line
%BEnRange
 0.00000 | 20.00000 | eV # [BSS] Energy range
%
specify the energy range where we want to calculate the dielectric response.
The lines
<pre>% BSEBands
2 | 7 | # [BSK] Bands range
%
BSENGBlk= 51 RL # [BSK] Screened interaction block size
BSENGexx= 773 RL # [BSK] Exchange components</pre>
specify the range of bands to be used and the RL components for the screened interaction (51) (read from the ndb.em1s file) and for the Hartree term (773, the maximal available). '''Remember''' that the BS kernel is written in Bloch space and its size is given by

<pre>BS kernel size = Valence Bands × Conduction Bands × K-points in the whole BZ</pre>
In our case the size is 2304. The BS bands range must be converged with care trying to keep as a few bands as possible.
Running

<pre> localhost:>yambo -F Inputs/06_BSE -J 06_BSE</pre>
will calculate the BS kernel and store in in the ''06_BSE'' folder. Before analyzing the BS absorption spectrum we can edit the log file ''l-06_BSE_optics_bse_bss'' to observe the on-the-fly live timing feature of yambo:

<pre> [...]
<---> BSK | | [000%] --(E) --(X)
<05s> BSK |# | [007%] 05s(E) 01m-07s(X)
<11s> BSK |### | [015%] 10s(E) 01m-06s(X)
[...]
<01m-02s> BSK |################## | [093%] 01m-01s(E) 01m-06s(X)
<01m-07s> BSK |####################| [100%] 01m-06s(E) 01m-06s(X)
[...]</pre>
The two numbers at the end of each line represents the elapsed time(E) and the expected time(X). The latter estimates the time needed to complete that particular section (the BS kernel construction in this case).

Finally we can compare the resulting BS optical absorption with the experimental result. We see how strong is the excitonic effect, if compared to the RPA calculation. The BS equation is, then, able to describe the bound electron-hole state responsible for the peak observed experimentally below the QP gap.

[[File:LiF_06.png|512x384px
]]

Notice that in the input there is variable called <code>BSEdiskIO</code>to tune the BSE-IO. If you have problem with IO on your cluster you can select "low" in this case the BSE matrix will be written only at the end of the calculation, restart will not be possible.

You can also decide to store only the results: eigelvalues, eigenvectors or spectrum is you select <code>BSEdiskIO</code><code>="none"</code><code>.</code>There is also the possibility to set it to 'high' in this case restart will be possible also in the exchange part of the calculations.

=A Bethe-Salpeter based Fxc(yambo -o c -k lrc -V qp) =

As for TDDFT we might be a little bit disappointed at this stage as the only method that yields a reasonable absorption spectrum for LiF is based on MBPT. It is possible, however, to extract from the BS kernel a possible approximation for the xc-kernel. This approach is efficient, but is far too complicated, and requires the prior computation of the BS kernel. An alternative possibility is to use a simple Long Range Component model, as introduced in [[LiF#Reining|Reining]]. In this model, the yambo input file ''Inputs/07_LRC'' describes a model where

[[File:LiF_PIC_LiF_equations-14.png|190x55px]]
with F0= ''LRC_alpha''=-8.7.

[[File:LiF_07.png|512x384px
]]

We see from the previous figure that the agreement with the BS calculation is not very good. Nevertheless this simple model captures the presence of an exciton and works quite well for a large number of simple semiconductors, as show in [[LiF#silvana|Botti]]:

[[File:F_silvana.png|371x352px]]



= Additional Exercises =

<ol>
<li>Calculate the absorption spectra like in section [[LiF#02|02]]-[[LiF#03|03]] increasing the polarization RL size, the polarization bands to find the converged values.</li>
<li>Plot the EELS spectrum obtained in section [[LiF#02|02]]-[[LiF#03|03]] and check the effect of increasing the polarization bands. Do a resonant only calculation (using the -V option) to see that the EELS cannot be described using only the resonant part of the response function.</li>
<li>See the effect of using a denser BZ sampling grid running Abinit with, for example
<pre>ngkpt2 6 6 6</pre></li>
<li>Repeat the [[LiF#06|BSE]] calculation increasing the bands.</li>
<li>Calculate the [[LiF#06|BSE]] spectrum using a diagonal only interaction.</li></ol>

= Detection of excitons symmetries (requires lumen 2.1 or better) =

Here we show how to detect excitons symmetries in LiF<ref name="Nalabothula2026"/> after extracting the BSE eigenvectors, either via full diagonalization or via recursive solvers.

TO DO ...

= References =

<references>
<ref name="wikipedia"> [http://en.wikipedia.org/wiki/Lithium_fluoride| LiF on wikipedia]</ref>
<ref name="Rohlfing"> M. Rohlfing, S. Louie, Phys. Rev. Lett. '''81''', 2312 (2000)</ref>
<ref name="Marini"> A. Marini, R. Del Sole, and A. Rubio, Phys. Rev. Lett., '''91''', 256402 (2003)</ref>
<ref name="Adragna"> G. Adragna, R. Del Sole, and A. Marini, Phys. Rev. B '''68''', 165108 (2003)</ref>
<ref name="Sottile"> F. Sottile, V. Olevano, and L. Reining, Phys. Rev. Lett. '''91''', 056402 (2003)</ref>
<ref name="Reining"> Reining, V. Olevano, A. Rubio, and G. Onida, Phys. Rev. Lett. '''88''', 066404 (2002)</ref>
<ref name="Botti"> S. Botti, Phys. Rev. B '''69''', 155112 (2004)</ref>
<ref name="Sangalli"> D. Sangalli, [https://doi.org/10.1103/PhysRevMaterials.5.083803 Phys. Rev. Mat. '''5''', 083803 (2021)]</ref>
<ref name="Nalabothula2026"> M. Nalabothula, D. Sangalli, F. Paleari, S. Reichardt, L. Wirtz, [https://doi.org/10.1103/7h2w-ys9h Phys. Rev. B. '''113''', 205130 (2026)] </ref>

</references>

Lumen 2.1

2026-06-16T08:27:19Z

Davide: Left only two highlighted main features

[[File:Times_logo_red.png|right|thumb|250px|TIMES logo red|Lumen 2.1 release was supported by the EU through the [https://times.uv.es/ TIMES MSCA-DN]]]

==Release Notes==

This is the second release of the lumen code, after the reboot based on yambo 5.3. A Partial sync with branch 5.4 of the yambo code was done.

===Main new features===

* Non-linear response with the Floquet formalism: in Lumen 2.1 it is possible to calculate non-linear response functions using Floquet formalism, for the moment this implementation is limited to the independent particle approximation with or without quasi-particle correction. See tutorial: [[SHG from Floquet formulation]]

* Exciton symmetries: Lumen 2.1 is now able to detect the point group of the system based on the list of symmetry operations. It is then possible to compute exciton symmetries via ypp after solving BSE with diagonalization solver.

* Exciton radiative lifetimes: old implementation of exciton radiative lifetimes in ypp has been recovered.
* Excitonic TR-ARPES: with ypp it is possible to compute TR-ARPES from BSE post-processing
* Exciton-exciton transitions: with ypp it is possible to compute exciton-exciton transitions from BSE post-processing
* Coded symmetry average of the dielectric constant for q=0 momentum
* Added interface to the CHASE library as BSE solver
* In spin polarized systems, it is now possible to get the contribution of a given spin-channel to absorption
* Added Gaussian smearing for the phonon-assisted response functions
* Added preliminary support for system with non-symmorphic symmetry operations. At present the database conversion can be done, and later the non-symmorphic symmetry operation can be removed via ypp
* Added support to Covariant dipoles in shifted grids with multiple shifts
* Plotting of exciton wave-function: it is now possible to plot both fixed electron and fixed hole case. Moreover it is possible to plot the coefficients of the excitonic WF in k-space. File bxsf is produced.

===Code optimization===
* Memory is now distributed in the old diagonalization solver (without using the ldiago library) interfaced with scalapack
* Fixed memory duplication in Slepc solver
* I/O in BSE can now be controlled by input as none/low/medium/high

===Fixes===

* Fixed compilation on MacOsX
* Fixed reading imaginary part from ndb.QP in yambo_nl
* Different compilation issues fixed: with intel compilers, hdf5, openACC, etc....
* Collisions table: added missing n_sp_pol variable

===More===
* added support to llvm compiler
* upgraded internal version of the following libraries: libxc, slepc, petsc, devxlib
* upgraded interface with libxc to work with version 7.0, and to handle different versions (with different API)
* upgraded interface with slepc/petsc to handle different versions (with different API)
* ExcPh databases now can be loaded also if requested band ranges is smaller than available one
* support for BSE with frozen states improved

Lumen release notes

2026-06-15T14:30:50Z

Davide: Add spacing between releases

This page contains the list of the releases of the lumen code.

Starting from version 2.0, the release notes are available in a dedicated page on the wiki, while the source code is available on gitlab.

== Lumen Releases ==

* Lumen 2.1: [[Lumen 2.1|Release Notes]]
** Tag 2.1.0 (XX/06/2026)

* Lumen 2.0: [[Lumen 2.0|Release Notes]]
** [https://gitlab.com/lumen-code/lumen/-/releases/2.0.1 Tag 2.0.1] (29/01/2026); [https://zenodo.org/records/20121198 Zenodo 20121198]
** [https://gitlab.com/lumen-code/lumen/-/releases/2.0.0 Tag 2.0.0] (04/12/2025);

== Previous versions ==

* [https://www.attaccalite.com/lumen/ Lumen 1.X] (2018)

Lumen release notes

2026-06-15T09:07:12Z

Davide: Reorganizing shape

This page contains the list of the releases of the lumen code.

Starting from version 2.0, the release notes are available in a dedicated page on the wiki, while the source code is available on gitlab.

== Lumen Releases ==

* Lumen 2.1: [[Lumen 2.1|Release Notes]]
** Tag 2.1.0 (XX/06/2026)
* Lumen 2.0: [[Lumen 2.0|Release Notes]]
** [https://gitlab.com/lumen-code/lumen/-/releases/2.0.1 Tag 2.0.1] (29/01/2026); [https://zenodo.org/records/20121198 Zenodo 20121198]
** [https://gitlab.com/lumen-code/lumen/-/releases/2.0.0 Tag 2.0.0] (04/12/2025);

== Previous versions ==

* [https://www.attaccalite.com/lumen/ Lumen 1.X] (2018)

Lumen 2.1

2026-06-11T13:15:03Z

Davide: /* Release Notes */

[[File:Times_logo_red.png|right|thumb|250px|TIMES logo red|Lumen 2.1 release was supported by the EU through the [https://times.uv.es/ TIMES MSCA-DN]]]

==Release Notes==

This is the second release of the lumen code, after the reboot based on yambo 5.3. A Partial sync with branch 5.4 of the yambo code was done.

===Main new features===

* Non-linear response with the Floquet formalism: in Lumen 2.1 it is possible to calculate non-linear response functions using Floquet formalism, for the moment this implementation is limited to the independent particle approximation with or without quasi-particle correction. See tutorial: [[SHG from Floquet formulation]]
* Exciton symmetries: Lumen 2.1 is now able to detect the point group of the system based on the list of symmetry operations. It is then possible to compute exciton symmetries via ypp after solving BSE with diagonalization solver.
* Exciton radiative lifetimes: old implementation of exciton radiative lifetimes in ypp has been recovered.
* Excitonic TR-ARPES: with ypp it is possible to compute TR-ARPES from BSE post-processing
* Exciton-exciton transitions: with ypp it is possible to compute exciton-exciton transitions from BSE post-processing
* Coded symmetry average of the dielectric constant for q=0 momentum
* Added interface to the CHASE library as BSE solver
* In spin polarized systems, it is now possible to get the contribution of a given spin-channel to absorption
* Added Gaussian smearing for the phonon-assisted response functions
* Added preliminary support for system with non-symmorphic symmetry operations. At present the database conversion can be done, and later the non-symmorphic symmetry operation can be removed via ypp
* Added support to Covariant dipoles in shifted grids with multiple shifts
* Plotting of exciton wave-function: it is now possible to plot both fixed electron and fixed hole case. Moreover it is possible to plot the coefficients of the excitonic WF in k-space. File bxsf is produced.

===Code optimization===
* Memory is now distributed in the old diagonalization solver (without using the ldiago library) interfaced with scalapack
* Fixed memory duplication in Slepc solver
* I/O in BSE can now be controlled by input as none/low/medium/high

===Fixes===

* Fixed compilation on MacOsX
* Fixed reading imaginary part from ndb.QP in yambo_nl
* Different compilation issues fixed: with intel compilers, hdf5, openACC, etc....
* Collisions table: added missing n_sp_pol variable

===More===
* added support to llvm compiler
* upgraded internal version of the following libraries: libxc, slepc, petsc, devxlib
* upgraded interface with libxc to work with version 7.0, and to handle different versions (with different API)
* upgraded interface with slepc/petsc to handle different versions (with different API)
* ExcPh databases now can be loaded also if requested band ranges is smaller than available one
* support for BSE with frozen states improved

Lumen 2.1

2026-06-11T13:13:03Z

Davide: Updating list of changes

[[File:Times_logo_red.png|right|thumb|250px|TIMES logo red|Lumen 2.1 release was supported by the EU through the [https://times.uv.es/ TIMES MSCA-DN]]]

==Release Notes==

This is the second release of the lumen code, after the reboot based on yambo 5.3

===Main new features===

* Non-linear response with the Floquet formalism: in Lumen 2.1 it is possible to calculate non-linear response functions using Floquet formalism, for the moment this implementation is limited to the independent particle approximation with or without quasi-particle correction. See tutorial: [[SHG from Floquet formulation]]
* Exciton symmetries: Lumen 2.1 is now able to detect the point group of the system based on the list of symmetry operations. It is then possible to compute exciton symmetries via ypp after solving BSE with diagonalization solver.
* Exciton radiative lifetimes: old implementation of exciton radiative lifetimes in ypp has been recovered.
* Excitonic TR-ARPES: with ypp it is possible to compute TR-ARPES from BSE post-processing
* Exciton-exciton transitions: with ypp it is possible to compute exciton-exciton transitions from BSE post-processing
* Coded symmetry average of the dielectric constant for q=0 momentum
* Added interface to the CHASE library as BSE solver
* In spin polarized systems, it is now possible to get the contribution of a given spin-channel to absorption
* Added Gaussian smearing for the phonon-assisted response functions
* Added preliminary support for system with non-symmorphic symmetry operations. At present the database conversion can be done, and later the non-symmorphic symmetry operation can be removed via ypp
* Added support to Covariant dipoles in shifted grids with multiple shifts
* Plotting of exciton wave-function: it is now possible to plot both fixed electron and fixed hole case. Moreover it is possible to plot the coefficients of the excitonic WF in k-space. File bxsf is produced.

===Code optimization===
* Memory is now distributed in the old diagonalization solver (without using the ldiago library) interfaced with scalapack
* Fixed memory duplication in Slepc solver
* I/O in BSE can now be controlled by input as none/low/medium/high

===Fixes===

* Fixed compilation on MacOsX
* Fixed reading imaginary part from ndb.QP in yambo_nl
* Different compilation issues fixed: with intel compilers, hdf5, openACC, etc....
* Collisions table: added missing n_sp_pol variable

===More===
* added support to llvm compiler
* upgraded internal version of the following libraries: libxc, slepc, petsc, devxlib
* upgraded interface with libxc to work with version 7.0, and to handle different versions (with different API)
* upgraded interface with slepc/petsc to handle different versions (with different API)
* ExcPh databases now can be loaded also if requested band ranges is smaller than available one
* support for BSE with frozen states improved

Lumen 2.1

2026-06-11T12:47:37Z

Davide: Updating list of changes

Dichroism in molecules

2026-05-19T07:35:15Z

Davide: Fixed reference

[[File:Chirality with hands.svg.png|x300px|frameless|right|Chirality representation with hands, from wikipedia]]

Circular dichroism (CD) is dichroism involving circularly polarized light, i.e., the differential absorption of left- and right-handed light<ref name="wiki_CD"/>.

<math>\Delta n(\omega) = n_R(\omega) - n_L(\omega)</math>

Dichroism in absorption can originate from an asymmetry in the dielectric tensor, which can be induced by the presence of either an external magnetic field or a spontaneous magnetization (MCD). This can be computed from the off-diagonal elements of the dielectric tensor. See the tutorial on [[The magneto-optical Kerr effect]] (MOKE). Another form of dichroism can be due to the chiral nature of molecules or materials<ref name="wiki_chiral"/> and goes under the name of natural circular dichroism (NCD). This latter can be related to the beta tensor and will be discussed in the present tutorial. See these references for NCD with the Lumen code: Ref.<ref name="Molteni2022a"/>, Ref.<ref name="Molteni2022b"/>

== Few equations ==

NCD is described by the trace of the so called beta tensor

<math>\Delta n(\omega) = \frac{8\pi N\omega}{3c} Tr[\beta_{\alpha\beta}(\omega)]</math>

which related the electronic dipole to the time derivative of the applied magnetic field, and the magnetic dipole to the time derivative of the applied electric field

<math>\mu_{\alpha}(\omega) = \alpha_{\alpha\beta}(\omega) E_\beta(\omega) - \frac{i \omega}{c}\beta_{\alpha\beta}(\omega) B_\beta(\omega)</math>

<math>m_{\alpha}(\omega) = \chi_{\alpha\beta}(\omega) B_\beta(\omega) + \frac{i \omega}{c}\beta_{\alpha\beta}(\omega) E_\beta(\omega)</math>

The latter can be obtained via the so called G tensor, which is similar to the dipole-dipole response function, but with electric and magnetic dipoles

<math>\beta_{\alpha\beta}(\omega) = -\frac{i\hbar m_e}{\omega q_e^2} G_{\alpha\beta}(\omega)</math>

<math>G_{ij}(\omega) = \frac{q_e^2}{2m\hbar}\sum_{nm}
\left( \frac{r^i_{nm} m^j_{mn}}{\omega_{nm}-\omega-i\gamma}
+ \frac{m^i_{mn} r^j_{nm}}{\omega_{nm}+\omega+i\gamma} \right)
</math>

The calculation of the magnetic dipoles is what makes NCD difficult to be computed in PBC, since these are not easily defined (solutions have been recently proposed in the literature). For isolated systems instead a simple sum over states can be employed (this is not possible in PBC because it involves diagonal intra-band dipoles): 

<math>
\mathbf{m}_{nm} = \sum_{l=l_{min}}^{l_{max}} \mathbf{r}_{nl} \times \mathbf{v}_{lm}
</math>

A detailed study of the convergence of the sum over states can be found in Ref.<ref name="Molteni2022a"/>. The same approach is also used in the Siesta code.

== Prerequisites ==
=== Databases ===
[[File:R methyloxirane structure.png|x150px|frameless|right|R methyloxirane structure]]
In this tutorial we will compute NCD on Propylene oxide, also known as R-Methyloxirane<ref name="wiki_PropOxide"/>. R-Methyloxirane has often been used in the literature as a benchmark molecule for NCD calculations against experimental data, due to its rigidity. To this end we provide a <code>yambo</code> database (DB) for R_methyloxirane.

* Yambo core databases (9 Mb): [https://media.yambo-code.eu/educational/tutorials/files/R_methyloxirane_DBs.tar.gz R_methyloxirane_DBs.tar.gz]



=== Input files ===
To compute NCD a standard calculation of response function in transition space needs to be performed, also setting "dich" (short for dichroism) and eventually "trace" in the "BSE_prop" variable as below
BSEprop= "abs dich trace" # [BSS] abs/kerr/magn/dichr trace

Since the dipoles are obtained from a sum over states it is important to converge the number of states used in the equation, this can be controlled via the input variables
DipBandsALL # [DIP] Compute all bands range, not only valence and conduction
DipApproach= "G-space v" # [DIP] [G-space v/R-space x/Covariant/Shifted grids]
% DipBands
1 | 100 | # [DIP] Bands range for dipoles
%

A standard input file for response function in transition space is (IP approximation)
optics # [R] Linear Response optical properties
bse # [R][BSE] Bethe Salpeter Equation.
dipoles # [R] Oscillator strenghts (or dipoles)
BSEmod= "retarded" # [BSE] resonant/retarded/coupling
BSKmod= "IP" # [BSE] IP/Hartree/HF/ALDA/SEX/BSfxc
% BSEBands
1 | 20 | # [BSK] Bands range
%
% BEnRange
4.0000 | 7.00000 | eV # [BSS] Energy range
%
% BDmRange
0.10000 | 0.10000 | eV # [BSS] Damping range
%
BEnSteps= 1001 # [BSS] Energy steps
% BLongDir
1.000000 | 0.000000 | 0.000000 | # [BSS] [cc] Electric Field
%

A part from this, it is suggested to set the input variable which helps the convergence in isolated systems
rim_cut # [R] Coulomb potential
NonPDirs= "XYZ" # [X/BSS] Non periodic chartesian directions (X,Y,Z,XY...)
RandQpts=0 # [RIM] Number of random q-points in the BZ
RandGvec= 1 RL # [RIM] Coulomb interaction RS components
CUTGeo= "ws xyz" # [CUT] Coulomb Cutoff geometry: box/cylinder/sphere/ws X/Y/Z/XY..
% CUTBox
0.000000 | 0.000000 | 0.000000 | # [CUT] [au] Box sides
%
CUTRadius= 0.000000 # [CUT] [au] Sphere/Cylinder radius
CUTCylLen= 0.000000 # [CUT] [au] Cylinder length
CUTwsGvec= 0.700000 # [CUT] WS cutoff: number of G to be modified

== Methyloxirane ==

=== CD at the Independent Particles approximation ===
Here we use it to show how to compute NCD using the yambo code.
Here we provide a poorly converged database, just to show how a calculation of NCD works with the yambo code.
Using the input file provided above you will obtain the following plots for absorption (alpha file, on the left) and NCD (CD file, on the right):
{|style="margin: 10px auto;"
| [[File:Methylox abs tuto.png|450px|frameless|center|Yambo tutorial, absorption in methyloxirane]]
| [[File:Methylox cd tuto.png|450px|frameless|center|Yambo tutorial, circular dichroism in methyloxirane]]
|}

Below we also report the results with converged parameters from <ref name="Molteni2022a"/>.
This simulation was done with an FCC supercell of size 50 and an energy cutoff of 90 Ry, while for the tutorial we used a supercell of size 20 and an energy cutoff of 70 Ry in the QE run.
Moreover a shift of X.XX eV was applied to align the spectra with the experimental data.
{|style="margin: 10px auto;"
| [[File:Absorption spectrum of methyloxirane.png|620px|frameless|center|Absorption spectrum of methyloxirane. Yambo results compared with experimental data.]]
| [[File:Circular Dichroism spectrum of methyloxirane.png|590px|frameless|center|Circular Dichroism spectrum of methyloxirane. Yambo results compared with experimental data.]]
|}

=== Beyond the Independent Particles approximation ===

== References ==

<references>
<ref name="wiki_CD">[https://en.wikipedia.org/wiki/Circular_dichroism Circular Dichroism] on wikipedia</ref>
<ref name="wiki_chiral">[https://en.wikipedia.org/wiki/Chiral_media Chiral Media] on wikipedia</ref>
<ref name="wiki_PropOxide">[https://en.wikipedia.org/wiki/Propylene_oxide Propylene oxide] on wikipedia</ref>
<ref name="Molteni2022a">E. Molteni, G. Cappellini, and D. Sangalli, ''Ab initio circular dichroism with the Yambo code: applications to dipeptides'' [https://www.doi.org/10.1088/1757-899X/1265/1/012005 IOP Conf. Series: Materials Science and Engineering 1265, 012005 (2022)]</ref>
<ref name="Molteni2022b"> E Molteni, G Mattioli, D Sangalli, ''Ab initio Circular Dichroism with the Yambo code: beyond the Independent Particle approximation'', [https://www.sif.it/riviste/sif/ncc/econtents/2022/045/06/article/18 Il Nuovo Cimento 4 (C), 175 (2022)], [https://arxiv.org/abs/2202.12702 arXiv2202.12720]</ref>
</references>

Dichroism in molecules

2026-05-19T07:34:39Z

Davide: All equations fixed

[[File:Chirality with hands.svg.png|x300px|frameless|right|Chirality representation with hands, from wikipedia]]

Circular dichroism (CD) is dichroism involving circularly polarized light, i.e., the differential absorption of left- and right-handed light<ref name="wiki_CD"/>.

<math>\Delta n(\omega) = n_R(\omega) - n_L(\omega)</math>

Dichroism in absorption can originate from an asymmetry in the dielectric tensor, which can be induced by the presence of either an external magnetic field or a spontaneous magnetization (MCD). This can be computed from the off-diagonal elements of the dielectric tensor. See the tutorial on [[The magneto-optical Kerr effect]] (MOKE). Another form of dichroism can be due to the chiral nature of molecules or materials<ref name="wiki_chiral"/> and goes under the name of natural circular dichroism (NCD). This latter can be related to the beta tensor and will be discussed in the present tutorial. See these references for NCD with the Lumen code: Ref.<ref name="Molteni2022a"/>, Ref.<ref name="Molteni2022b"/>

== Few equations ==

NCD is described by the trace of the so called beta tensor

<math>\Delta n(\omega) = \frac{8\pi N\omega}{3c} Tr[\beta_{\alpha\beta}(\omega)]</math>

which related the electronic dipole to the time derivative of the applied magnetic field, and the magnetic dipole to the time derivative of the applied electric field

<math>\mu_{\alpha}(\omega) = \alpha_{\alpha\beta}(\omega) E_\beta(\omega) - \frac{i \omega}{c}\beta_{\alpha\beta}(\omega) B_\beta(\omega)</math>

<math>m_{\alpha}(\omega) = \chi_{\alpha\beta}(\omega) B_\beta(\omega) + \frac{i \omega}{c}\beta_{\alpha\beta}(\omega) E_\beta(\omega)</math>

The latter can be obtained via the so called G tensor, which is similar to the dipole-dipole response function, but with electric and magnetic dipoles

<math>\beta_{\alpha\beta}(\omega) = -\frac{i\hbar m_e}{\omega q_e^2} G_{\alpha\beta}(\omega)</math>

<math>G_{ij}(\omega) = \frac{q_e^2}{2m\hbar}\sum_{nm}
\left( \frac{r^i_{nm} m^j_{mn}}{\omega_{nm}-\omega-i\gamma}
+ \frac{m^i_{mn} r^j_{nm}}{\omega_{nm}+\omega+i\gamma} \right)
</math>

The calculation of the magnetic dipoles is what makes NCD difficult to be computed in PBC, since these are not easily defined (solutions have been recently proposed in the literature). For isolated systems instead a simple sum over states can be employed (this is not possible in PBC because it involves diagonal intra-band dipoles): 

<math>
\mathbf{m}_{nm} = \sum_{l=l_{min}}^{l_{max}} \mathbf{r}_{nl} \times \mathbf{v}_{lm}
</math>

A detailed study of the convergence of the sum over states can be found in Ref.<ref name="Molteni2022a"/>. The same approach is also used in the Siesta code.

== Prerequisites ==
=== Databases ===
[[File:R methyloxirane structure.png|x150px|frameless|right|R methyloxirane structure]]
In this tutorial we will compute NCD on Propylene oxide, also known as R-Methyloxirane<ref name="wiki_PropOxide"/>. R-Methyloxirane has often been used in the literature as a benchmark molecule for NCD calculations against experimental data, due to its rigidity. To this end we provide a <code>yambo</code> database (DB) for R_methyloxirane.

* Yambo core databases (9 Mb): [https://media.yambo-code.eu/educational/tutorials/files/R_methyloxirane_DBs.tar.gz R_methyloxirane_DBs.tar.gz]



=== Input files ===
To compute NCD a standard calculation of response function in transition space needs to be performed, also setting "dich" (short for dichroism) and eventually "trace" in the "BSE_prop" variable as below
BSEprop= "abs dich trace" # [BSS] abs/kerr/magn/dichr trace

Since the dipoles are obtained from a sum over states it is important to converge the number of states used in the equation, this can be controlled via the input variables
DipBandsALL # [DIP] Compute all bands range, not only valence and conduction
DipApproach= "G-space v" # [DIP] [G-space v/R-space x/Covariant/Shifted grids]
% DipBands
1 | 100 | # [DIP] Bands range for dipoles
%

A standard input file for response function in transition space is (IP approximation)
optics # [R] Linear Response optical properties
bse # [R][BSE] Bethe Salpeter Equation.
dipoles # [R] Oscillator strenghts (or dipoles)
BSEmod= "retarded" # [BSE] resonant/retarded/coupling
BSKmod= "IP" # [BSE] IP/Hartree/HF/ALDA/SEX/BSfxc
% BSEBands
1 | 20 | # [BSK] Bands range
%
% BEnRange
4.0000 | 7.00000 | eV # [BSS] Energy range
%
% BDmRange
0.10000 | 0.10000 | eV # [BSS] Damping range
%
BEnSteps= 1001 # [BSS] Energy steps
% BLongDir
1.000000 | 0.000000 | 0.000000 | # [BSS] [cc] Electric Field
%

A part from this, it is suggested to set the input variable which helps the convergence in isolated systems
rim_cut # [R] Coulomb potential
NonPDirs= "XYZ" # [X/BSS] Non periodic chartesian directions (X,Y,Z,XY...)
RandQpts=0 # [RIM] Number of random q-points in the BZ
RandGvec= 1 RL # [RIM] Coulomb interaction RS components
CUTGeo= "ws xyz" # [CUT] Coulomb Cutoff geometry: box/cylinder/sphere/ws X/Y/Z/XY..
% CUTBox
0.000000 | 0.000000 | 0.000000 | # [CUT] [au] Box sides
%
CUTRadius= 0.000000 # [CUT] [au] Sphere/Cylinder radius
CUTCylLen= 0.000000 # [CUT] [au] Cylinder length
CUTwsGvec= 0.700000 # [CUT] WS cutoff: number of G to be modified

== Methyloxirane ==

=== CD at the Independent Particles approximation ===
Here we use it to show how to compute NCD using the yambo code.
Here we provide a poorly converged database, just to show how a calculation of NCD works with the yambo code.
Using the input file provided above you will obtain the following plots for absorption (alpha file, on the left) and NCD (CD file, on the right):
{|style="margin: 10px auto;"
| [[File:Methylox abs tuto.png|450px|frameless|center|Yambo tutorial, absorption in methyloxirane]]
| [[File:Methylox cd tuto.png|450px|frameless|center|Yambo tutorial, circular dichroism in methyloxirane]]
|}

Below we also report the results with converged parameters from <ref name="Molteni2022"/>.
This simulation was done with an FCC supercell of size 50 and an energy cutoff of 90 Ry, while for the tutorial we used a supercell of size 20 and an energy cutoff of 70 Ry in the QE run.
Moreover a shift of X.XX eV was applied to align the spectra with the experimental data.
{|style="margin: 10px auto;"
| [[File:Absorption spectrum of methyloxirane.png|620px|frameless|center|Absorption spectrum of methyloxirane. Yambo results compared with experimental data.]]
| [[File:Circular Dichroism spectrum of methyloxirane.png|590px|frameless|center|Circular Dichroism spectrum of methyloxirane. Yambo results compared with experimental data.]]
|}

=== Beyond the Independent Particles approximation ===

== References ==

<references>
<ref name="wiki_CD">[https://en.wikipedia.org/wiki/Circular_dichroism Circular Dichroism] on wikipedia</ref>
<ref name="wiki_chiral">[https://en.wikipedia.org/wiki/Chiral_media Chiral Media] on wikipedia</ref>
<ref name="wiki_PropOxide">[https://en.wikipedia.org/wiki/Propylene_oxide Propylene oxide] on wikipedia</ref>
<ref name="Molteni2022a">E. Molteni, G. Cappellini, and D. Sangalli, ''Ab initio circular dichroism with the Yambo code: applications to dipeptides'' [https://www.doi.org/10.1088/1757-899X/1265/1/012005 IOP Conf. Series: Materials Science and Engineering 1265, 012005 (2022)]</ref>
<ref name="Molteni2022b"> E Molteni, G Mattioli, D Sangalli, ''Ab initio Circular Dichroism with the Yambo code: beyond the Independent Particle approximation'', [https://www.sif.it/riviste/sif/ncc/econtents/2022/045/06/article/18 Il Nuovo Cimento 4 (C), 175 (2022)], [https://arxiv.org/abs/2202.12702 arXiv2202.12720]</ref>
</references>

Dichroism in molecules

2026-05-19T07:28:13Z

Davide: Fixing some equations

[[File:Chirality with hands.svg.png|x300px|frameless|right|Chirality representation with hands, from wikipedia]]

Circular dichroism (CD) is dichroism involving circularly polarized light, i.e., the differential absorption of left- and right-handed light<ref name="wiki_CD"/>.

<math>\Delta n(\omega) = n_R(\omega) - n_L(\omega)</math>

Dichroism in absorption can originate from an asymmetry in the dielectric tensor, which can be induced by the presence of either an external magnetic field or a spontaneous magnetization (MCD). This can be computed from the off-diagonal elements of the dielectric tensor. See the tutorial on [[The magneto-optical Kerr effect]] (MOKE). Another form of dichroism can be due to the chiral nature of molecules or materials<ref name="wiki_chiral"/> and goes under the name of natural circular dichroism (NCD). This latter can be related to the beta tensor and will be discussed in the present tutorial. See these references for NCD with the Lumen code: Ref.<ref name="Molteni2022a"/>, Ref.<ref name="Molteni2022b"/>

== Few equations ==

NCD is described by the trace of the so called beta tensor

<math>\Delta n(\omega) = \frac{8\pi N\omega}{3c} Tr[\beta_{\alpha\beta}(\omega)]</math>

which related the electronic dipole to the time derivative of the applied magnetic field, and the magnetic dipole to the time derivative of the applied electric field

<math>\mu_{\alpha}(\omega) = \alpha_{\alpha\beta}(\omega) E_\beta(\omega) - \frac{i \omega}{c}\beta_{\alpha\beta}(\omega) B_\beta(\omega)</math>

<math>m_{\alpha}(\omega) = \chi_{\alpha\beta}(\omega) B_\beta(\omega) + \frac{i \omega}{c}\beta_{\alpha\beta}(\omega) E_\beta(\omega)</math>

The latter can be obtained via the so called G tensor, which is similar to the dipole-dipole response function, but with electric and magnetic dipoles

[[File:Beta G tensor relation.png|x60px|frameless|Beta G tensor relation]] 
[[File:Gtensor.png|x120px|frameless|Yambo tutorial image]] 

The calculation of the magnetic dipoles is what makes NCD difficult to be computed in PBC, since these are not easily defined (solutions have been recently proposed in the literature). For isolated systems instead a simple sum over states can be employed (this is not possible in PBC because it involves diagonal intra-band dipoles): 
[[File:Mdipoles sum over states.png|x70px|frameless|Yambo tutorial image]]

A detailed study of the convergence of the sum over states can be found in Ref.<ref name="Molteni2022a"/>. The same approach is also used in the Siesta code.

== Prerequisites ==
=== Databases ===
[[File:R methyloxirane structure.png|x150px|frameless|right|R methyloxirane structure]]
In this tutorial we will compute NCD on Propylene oxide, also known as R-Methyloxirane<ref name="wiki_PropOxide"/>. R-Methyloxirane has often been used in the literature as a benchmark molecule for NCD calculations against experimental data, due to its rigidity. To this end we provide a <code>yambo</code> database (DB) for R_methyloxirane.

* Yambo core databases (9 Mb): [https://media.yambo-code.eu/educational/tutorials/files/R_methyloxirane_DBs.tar.gz R_methyloxirane_DBs.tar.gz]



=== Input files ===
To compute NCD a standard calculation of response function in transition space needs to be performed, also setting "dich" (short for dichroism) and eventually "trace" in the "BSE_prop" variable as below
BSEprop= "abs dich trace" # [BSS] abs/kerr/magn/dichr trace

Since the dipoles are obtained from a sum over states it is important to converge the number of states used in the equation, this can be controlled via the input variables
DipBandsALL # [DIP] Compute all bands range, not only valence and conduction
DipApproach= "G-space v" # [DIP] [G-space v/R-space x/Covariant/Shifted grids]
% DipBands
1 | 100 | # [DIP] Bands range for dipoles
%

A standard input file for response function in transition space is (IP approximation)
optics # [R] Linear Response optical properties
bse # [R][BSE] Bethe Salpeter Equation.
dipoles # [R] Oscillator strenghts (or dipoles)
BSEmod= "retarded" # [BSE] resonant/retarded/coupling
BSKmod= "IP" # [BSE] IP/Hartree/HF/ALDA/SEX/BSfxc
% BSEBands
1 | 20 | # [BSK] Bands range
%
% BEnRange
4.0000 | 7.00000 | eV # [BSS] Energy range
%
% BDmRange
0.10000 | 0.10000 | eV # [BSS] Damping range
%
BEnSteps= 1001 # [BSS] Energy steps
% BLongDir
1.000000 | 0.000000 | 0.000000 | # [BSS] [cc] Electric Field
%

A part from this, it is suggested to set the input variable which helps the convergence in isolated systems
rim_cut # [R] Coulomb potential
NonPDirs= "XYZ" # [X/BSS] Non periodic chartesian directions (X,Y,Z,XY...)
RandQpts=0 # [RIM] Number of random q-points in the BZ
RandGvec= 1 RL # [RIM] Coulomb interaction RS components
CUTGeo= "ws xyz" # [CUT] Coulomb Cutoff geometry: box/cylinder/sphere/ws X/Y/Z/XY..
% CUTBox
0.000000 | 0.000000 | 0.000000 | # [CUT] [au] Box sides
%
CUTRadius= 0.000000 # [CUT] [au] Sphere/Cylinder radius
CUTCylLen= 0.000000 # [CUT] [au] Cylinder length
CUTwsGvec= 0.700000 # [CUT] WS cutoff: number of G to be modified

== Methyloxirane ==

=== CD at the Independent Particles approximation ===
Here we use it to show how to compute NCD using the yambo code.
Here we provide a poorly converged database, just to show how a calculation of NCD works with the yambo code.
Using the input file provided above you will obtain the following plots for absorption (alpha file, on the left) and NCD (CD file, on the right):
{|style="margin: 10px auto;"
| [[File:Methylox abs tuto.png|450px|frameless|center|Yambo tutorial, absorption in methyloxirane]]
| [[File:Methylox cd tuto.png|450px|frameless|center|Yambo tutorial, circular dichroism in methyloxirane]]
|}

Below we also report the results with converged parameters from <ref name="Molteni2022"/>.
This simulation was done with an FCC supercell of size 50 and an energy cutoff of 90 Ry, while for the tutorial we used a supercell of size 20 and an energy cutoff of 70 Ry in the QE run.
Moreover a shift of X.XX eV was applied to align the spectra with the experimental data.
{|style="margin: 10px auto;"
| [[File:Absorption spectrum of methyloxirane.png|620px|frameless|center|Absorption spectrum of methyloxirane. Yambo results compared with experimental data.]]
| [[File:Circular Dichroism spectrum of methyloxirane.png|590px|frameless|center|Circular Dichroism spectrum of methyloxirane. Yambo results compared with experimental data.]]
|}

=== Beyond the Independent Particles approximation ===

== References ==

<references>
<ref name="wiki_CD">[https://en.wikipedia.org/wiki/Circular_dichroism Circular Dichroism] on wikipedia</ref>
<ref name="wiki_chiral">[https://en.wikipedia.org/wiki/Chiral_media Chiral Media] on wikipedia</ref>
<ref name="wiki_PropOxide">[https://en.wikipedia.org/wiki/Propylene_oxide Propylene oxide] on wikipedia</ref>
<ref name="Molteni2022a">E. Molteni, G. Cappellini, and D. Sangalli, ''Ab initio circular dichroism with the Yambo code: applications to dipeptides'' [https://www.doi.org/10.1088/1757-899X/1265/1/012005 IOP Conf. Series: Materials Science and Engineering 1265, 012005 (2022)]</ref>
<ref name="Molteni2022b"> E Molteni, G Mattioli, D Sangalli, ''Ab initio Circular Dichroism with the Yambo code: beyond the Independent Particle approximation'', [https://www.sif.it/riviste/sif/ncc/econtents/2022/045/06/article/18 Il Nuovo Cimento 4 (C), 175 (2022)], [https://arxiv.org/abs/2202.12702 arXiv2202.12720]</ref>
</references>

Dichroism in molecules

2026-05-19T07:23:51Z

Davide: Fixed initial equation

[[File:Chirality with hands.svg.png|x300px|frameless|right|Chirality representation with hands, from wikipedia]]

Circular dichroism (CD) is dichroism involving circularly polarized light, i.e., the differential absorption of left- and right-handed light<ref name="wiki_CD"/>.

<math>\Delta n = n_R-n_L</math>

Dichroism in absorption can originate from an asymmetry in the dielectric tensor, which can be induced by the presence of either an external magnetic field or a spontaneous magnetization (MCD). This can be computed from the off-diagonal elements of the dielectric tensor. See the tutorial on [[The magneto-optical Kerr effect]] (MOKE). Another form of dichroism can be due to the chiral nature of molecules or materials<ref name="wiki_chiral"/> and goes under the name of natural circular dichroism (NCD). This latter can be related to the beta tensor and will be discussed in the present tutorial. See these references for NCD with the Lumen code: Ref.<ref name="Molteni2022a"/>, Ref.<ref name="Molteni2022b"/>

== Few equations ==

NCD is described by the trace of the so called beta tensor

[[File:Deltan Beta relation.png|x60px|frameless|Yambo tutorial image]]

which related the electronic dipole to the time derivative of the applied magnetic field, and the magnetic dipole to the time derivative of the applied electric field

[[File:Dipoles fields relations.png|x120px|frameless|Yambo tutorial image]]

The latter can be obtained via the so called G tensor, which is similar to the dipole-dipole response function, but with electric and magnetic dipoles

[[File:Beta G tensor relation.png|x60px|frameless|Beta G tensor relation]] 
[[File:Gtensor.png|x120px|frameless|Yambo tutorial image]] 

The calculation of the magnetic dipoles is what makes NCD difficult to be computed in PBC, since these are not easily defined (solutions have been recently proposed in the literature). For isolated systems instead a simple sum over states can be employed (this is not possible in PBC because it involves diagonal intra-band dipoles): 
[[File:Mdipoles sum over states.png|x70px|frameless|Yambo tutorial image]]

A detailed study of the convergence of the sum over states can be found in Ref.<ref name="Molteni2022a"/>. The same approach is also used in the Siesta code.

== Prerequisites ==
=== Databases ===
[[File:R methyloxirane structure.png|x150px|frameless|right|R methyloxirane structure]]
In this tutorial we will compute NCD on Propylene oxide, also known as R-Methyloxirane<ref name="wiki_PropOxide"/>. R-Methyloxirane has often been used in the literature as a benchmark molecule for NCD calculations against experimental data, due to its rigidity. To this end we provide a <code>yambo</code> database (DB) for R_methyloxirane.

* Yambo core databases (9 Mb): [https://media.yambo-code.eu/educational/tutorials/files/R_methyloxirane_DBs.tar.gz R_methyloxirane_DBs.tar.gz]



=== Input files ===
To compute NCD a standard calculation of response function in transition space needs to be performed, also setting "dich" (short for dichroism) and eventually "trace" in the "BSE_prop" variable as below
BSEprop= "abs dich trace" # [BSS] abs/kerr/magn/dichr trace

Since the dipoles are obtained from a sum over states it is important to converge the number of states used in the equation, this can be controlled via the input variables
DipBandsALL # [DIP] Compute all bands range, not only valence and conduction
DipApproach= "G-space v" # [DIP] [G-space v/R-space x/Covariant/Shifted grids]
% DipBands
1 | 100 | # [DIP] Bands range for dipoles
%

A standard input file for response function in transition space is (IP approximation)
optics # [R] Linear Response optical properties
bse # [R][BSE] Bethe Salpeter Equation.
dipoles # [R] Oscillator strenghts (or dipoles)
BSEmod= "retarded" # [BSE] resonant/retarded/coupling
BSKmod= "IP" # [BSE] IP/Hartree/HF/ALDA/SEX/BSfxc
% BSEBands
1 | 20 | # [BSK] Bands range
%
% BEnRange
4.0000 | 7.00000 | eV # [BSS] Energy range
%
% BDmRange
0.10000 | 0.10000 | eV # [BSS] Damping range
%
BEnSteps= 1001 # [BSS] Energy steps
% BLongDir
1.000000 | 0.000000 | 0.000000 | # [BSS] [cc] Electric Field
%

A part from this, it is suggested to set the input variable which helps the convergence in isolated systems
rim_cut # [R] Coulomb potential
NonPDirs= "XYZ" # [X/BSS] Non periodic chartesian directions (X,Y,Z,XY...)
RandQpts=0 # [RIM] Number of random q-points in the BZ
RandGvec= 1 RL # [RIM] Coulomb interaction RS components
CUTGeo= "ws xyz" # [CUT] Coulomb Cutoff geometry: box/cylinder/sphere/ws X/Y/Z/XY..
% CUTBox
0.000000 | 0.000000 | 0.000000 | # [CUT] [au] Box sides
%
CUTRadius= 0.000000 # [CUT] [au] Sphere/Cylinder radius
CUTCylLen= 0.000000 # [CUT] [au] Cylinder length
CUTwsGvec= 0.700000 # [CUT] WS cutoff: number of G to be modified

== Methyloxirane ==

=== CD at the Independent Particles approximation ===
Here we use it to show how to compute NCD using the yambo code.
Here we provide a poorly converged database, just to show how a calculation of NCD works with the yambo code.
Using the input file provided above you will obtain the following plots for absorption (alpha file, on the left) and NCD (CD file, on the right):
{|style="margin: 10px auto;"
| [[File:Methylox abs tuto.png|450px|frameless|center|Yambo tutorial, absorption in methyloxirane]]
| [[File:Methylox cd tuto.png|450px|frameless|center|Yambo tutorial, circular dichroism in methyloxirane]]
|}

Below we also report the results with converged parameters from <ref name="Molteni2022"/>.
This simulation was done with an FCC supercell of size 50 and an energy cutoff of 90 Ry, while for the tutorial we used a supercell of size 20 and an energy cutoff of 70 Ry in the QE run.
Moreover a shift of X.XX eV was applied to align the spectra with the experimental data.
{|style="margin: 10px auto;"
| [[File:Absorption spectrum of methyloxirane.png|620px|frameless|center|Absorption spectrum of methyloxirane. Yambo results compared with experimental data.]]
| [[File:Circular Dichroism spectrum of methyloxirane.png|590px|frameless|center|Circular Dichroism spectrum of methyloxirane. Yambo results compared with experimental data.]]
|}

=== Beyond the Independent Particles approximation ===

== References ==

<references>
<ref name="wiki_CD">[https://en.wikipedia.org/wiki/Circular_dichroism Circular Dichroism] on wikipedia</ref>
<ref name="wiki_chiral">[https://en.wikipedia.org/wiki/Chiral_media Chiral Media] on wikipedia</ref>
<ref name="wiki_PropOxide">[https://en.wikipedia.org/wiki/Propylene_oxide Propylene oxide] on wikipedia</ref>
<ref name="Molteni2022a">E. Molteni, G. Cappellini, and D. Sangalli, ''Ab initio circular dichroism with the Yambo code: applications to dipeptides'' [https://www.doi.org/10.1088/1757-899X/1265/1/012005 IOP Conf. Series: Materials Science and Engineering 1265, 012005 (2022)]</ref>
<ref name="Molteni2022b"> E Molteni, G Mattioli, D Sangalli, ''Ab initio Circular Dichroism with the Yambo code: beyond the Independent Particle approximation'', [https://www.sif.it/riviste/sif/ncc/econtents/2022/045/06/article/18 Il Nuovo Cimento 4 (C), 175 (2022)], [https://arxiv.org/abs/2202.12702 arXiv2202.12720]</ref>
</references>