View source for LiF


= 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 f<sub>xc</sub> (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]. <br>
'''Important''': <span style="color:red"> Increase the number of bands in the NSCF DFT up to 50</span>

== 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] &amp; Occupations
 <---> [03] Transferred momenta grid: Indexes
 <---> [04] Game Over &amp; 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] &amp; 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.

<span id="02"></span><span id="03"></span>

=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 &quot;yambo -o c&quot; 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 51<sup>st</sup> 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]]

<span id="04"></span>
= 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= <span style="color:red">51 </span>                RL    # [TDDFT] XC-kernel RL size
 NGsBlkXd= <span style="color:red">51 </span>               RL    # [Xd] Response block size
 % QpntsRXd
 <span style="color:red"> 1 | 1 | </span>                    # [Xd] Transferred momenta
 %
 % BndsRnXd
 <span style="color:red">   1 | 10 | </span>                  # [Xd] Polarization function bands
 %
 % EnRngeXd
 <span style="color:red">  7.50000 | 25.00000 | </span>  eV  # [Xd] Energy range
 %
 % DmRngeXd
  0.100000 | 0.100000 |   eV  # [Xd] Damping range
 %
 ETStpsXd= <span style="color:red">300 </span>               # [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).

<span id="05"></span>

=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.

<span id="06"></span>

=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=  <span style="color:red"> 773</span>            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
   <span style="color:red"> 5.800000 </span> | 1.000000 | 1.000000 |        # [EXTQP BSK BSS] E parameters  (c/v) eV|adim|adim
 %
 % BSEQptR
  1 | 1 |                             # [BSK] Transferred momenta range
 %
 % BSEBands
  <span style="color:red">   2 |  7 |  </span>                         # [BSK] Bands range
 %
 % BEnRange
   <span style="color:red">  0.00000 | 20.00000 |   </span>      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= &quot;SEX&quot;                # [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
   <span style="color:red">  0.00000 | 20.00000 |   </span>      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 F<sub>xc</sub>(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 F<sub>0</sub>= ''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]]

<span id="exercises"></span>

= 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><p>See the effect of using a denser BZ sampling grid running Abinit with, for example</p>
<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>