by Benjamin Aurich, Christian Vorwerk, Caterina Cocchi, & Keith Gilmore, for exciting neon
Purpose: In this tutorial you will learn how to perform a basic BetheSalpeterequation (BSE) calculation. As an example, the optical spectrum of lithium fluoride (LiF) will be studied.
Table of Contents

1. Theoretical background: BetheSalpeter equation
2. Preliminary step: Groundstate calculation
Before starting, be sure that relevant environment variables are already defined as specified in How to set environment variables for tutorials scripts.
Important note: All input parameters are given in atomic units!
As a preliminary step to calculate excitedstate properties from BSE, a groundstate calculation has to be performed. In this tutorial we consider as an example LiF. Create a directory named LiFBSE and move into it.
$ mkdir LiFBSE
$ cd LiFBSE
Inside the directory LiFBSE we create a subdirectory GS where we perform the preliminary ground state calculation:
$ mkdir GS
$ cd GS
Inside the GS subdirectory we create the input file for LiF. In the structure element we include the lattice parameter and basis vectors of LiF, which has a rocksalt cubic lattice, as well as the positions of the Li and F atoms. In the groundstate element, we include a 8$\times$8$\times$8 kpoint mesh (ngridk) and a value of 14.0 for gmaxvr. This value, which is larger than the default, is needed in view of the excitedstate calculation.
Copy and paste the following into the new file input.xml:
<input> <title>LiFBSE</title> <structure speciespath="$EXCITINGROOT/species/"> <crystal scale="7.608"> <basevect>0.5 0.5 0.0</basevect> <basevect>0.5 0.0 0.5</basevect> <basevect>0.0 0.5 0.5</basevect> </crystal> <species speciesfile="Li.xml"> <atom coord="0.0000 0.0000 0.0000" /> </species> <species speciesfile="F.xml"> <atom coord="0.5000 0.5000 0.5000" /> </species> </structure> <groundstate do="fromscratch" ngridk="8 8 8" xctype="GGA_PBE_SOL" gmaxvr="14.0"/> </input>
N. B.: Do not forget to replace the string "$EXCITINGROOT" in input.xml by the actual value of the environment variable $EXCITINGROOT by running
$ SETUPexcitingroot.sh
Now start the groundstate SCF calculation and check if it finishes gracefully.
$ time exciting_smp
In case of a successful run the files STATE.OUT and EFERMI.OUT should be present in the directory. These two files are needed as the starting point for the BSE calculation.
3. How to run a BSE calculation
i) Preparation of the input file
We start by performing a BSE calculation using the TDA. To do this, move to the parent directory LiFBSE and create a new folder which we name BSE. Then, we move into it:
$ cd ..
$ mkdir BSE
$ cd BSE
Copy the necessary groundstate files from the GS folder into the BSE path with
$ cp ../GS/{STATE.OUT,EFERMI.OUT,input.xml} ./
We can then skip the groundstate calculation by setting do = "skip" in the groundstate element of the input file.
To perform an excitedstate calculation we must include the xs element in the input file:
... <xs xstype="BSE" ngridk="3 3 3" vkloff="0.097 0.273 0.493" ngridq="3 3 3" nempty="30" gqmax="2.5" broad="0.007" scissor="0.20947" tappinfo="true" tevout="true"> <energywindow intv="0.0 1.0" points="1200"/> <screening screentype="full" nempty="100"/> <BSE bsetype="singlet" nstlbse="1 5 1 4" /> <qpointset> <qpoint>0.0 0.0 0.0</qpoint> </qpointset> </xs> ...
Now, we consider the attributes that appear in this section of the input file. We will go through all of them in order to clarify their meaning:
 xstype defines the type of excitedstate calculation: here we set it to BSE;
 ngridk within the xs element is used in the nonselfconsistentfield calculation required to obtain the singleparticle eigenstates and eigenenergies. If not defined, it defaults to the kpoint mesh used in the groundstate calculation. Calculations should be converged with respect to the kmesh;
 vkloff is used to shift the kmesh off symmetry by a small displacement in order to break all symmetry relations among the kpoints of the mesh. In this way, all the kpoints in the mesh will be crystallographically inequivalent and there will be no redundant contribution to the spectrum;
 ngridq defines the qmesh for the calculation of the screening. It is a good practice to choose a qmesh equivalent to the kmesh;
 nempty determines the number of empty states used in the construction of the BSE Hamiltonian;
 gqmax is an energy threshold for the local field effects included in the calculation. The number of Gvectors employed in the calculation is written in the first line of the file GQPOINTS_QM001.OUT. Tuning gqmax may require adjusting the value of gmaxvr in the groundstate element and rerunning the groundstate calculation;
 broad defines the Lorentzian broadening $\delta$ for the calculated spectra;
 When the attribute scissor is set different from zero, a scissors operator is applied to correct the band gap obtained from the ground state calculation in order to mimic the quasiparticle gap. This parameter should not appear when the electronic structure resulting from a GW calculation is taken as starting point for the BSE calculation;
 tappinfo causes output to be appended to the log file INFOXS.OUT without overwriting the existing file.
 tevout sets the energy output to electronvolt (eV).
Additional subelements appear inside the xs element (refer to the Input Reference for further details):
 energywindow defines the energy window and the number of points for the calculation of the optical spectrum. The attribute intv indicates the lower and upper energy limits for the calculation of the spectrum, while points defines the number of energy points to be used.
 screening defines the parameters used to compute the screened Coulomb potential. The attribute screeningtype selects the approximation for the screening: in this case we include the full screening in the calculation. For other options refer to the Input Reference. The attribute nempty defines the number of empty states to be included in the calculation of the screening matrix. Note that this is a different attribute than the one within the xs element and must always be specified when performing a BSE calculation.
 In the BSE element the actual parameters for a BSE calculation are set. The attribute bsetype defines the level of approximation in the solution of the BSE Hamiltonian. Here we are calculating the "singlet" states, meaning that the full BSE Hamiltonian is diagonalized. Other options include the calculations of "triplet", where the exchange term of the Hamiltonian is neglected, RPA, where the direct term is ignored, and the independentparticle (IP) approximation, where only the diagonal part of the Hamiltonian is considered. For further details see Input Reference. The attribute nstlbse defines the range of occupied and empty states included in the BSE calculation. The first two integers define the lowest and highest occupied bands, while the last two refer to the unoccupied bands. To define the band range we refer to the file EIGVAL.OUT obtained from the ground state calculation. Please notice that the occupied bands are numbered starting from the lowest occupied state, while unoccupied bands are numbered starting from the lowest unoccupied band. In this example, we include in the BSE calculation all the 5 occupied bands of LiF and only the first 4 empty ones.
 The element qpointset is related to the dependence of the response function upon the momentum transfer along different directions (see qdependent TDDFT). Here, the qpoint is set to "0.0 0.0 0.0" since we are neglecting momentum transfer in this calculation ($\mathbf{q}=0$).
ii) Additional information about the calculation workflow
iii) Running exciting and analysing the results
We can now run the BSE calculation:
$ time exciting_smp
After the calculation has finished, we proceed with the analysis of the results. A number of files and folders are present in the working directory. Most of them contain technical information about the calculation that is not strictly related to the physical interpretation of the results. We are interested in plotting the optical absorption spectrum of LiF, given by the imaginary part of its macroscopic dielectric function. Therefore, we consider the results contained in the folder EPSILON. In order to plot the spectrum, execute the following commands:
$ cp EPSILON/EPSILON_BSEsingletTDABAR_SCRfull_OC11.OUT singletTDA
$ PLOTfiles.py f singletTDA lx 'Energy [eV]' ly 'Im $\varepsilon_M$' t 'Macroscopic dielectric function' g rc cy 3 x 0 27 nl
The result is stored in PLOT.png and PLOT.eps and should resemble this plot:
Notice that the energy scale is given in units of eV.
Since we are dealing with a cubic system, all optical components are equivalent. For this reason, we consider only the first one along the xx direction (OC11). The above command plots the imaginary part (indicated by the i after the filename) of the macroscopic dielectric function, contained in the third column, versus energy in eV, in the first column. The file EPSILON_BSEsingletTDABAR_SCRfull_OC11.OUT in the directory EPSILON also stores the real part of the macroscopic dielectric function calculated directly and via a KramersKronig transformation of the imaginary part in the second and fourth columns, respectively. The real part may be plotted similarly to above by replacing the cy 3 after the filename with cy 2 (for the real part).
iv) Scaling and convergence
BSE calculations are extremely demanding, even on stateoftheart computer infrastructure. The scaling with respect to the size of the kmesh is quadratic with respect to the setup of the Hamlitonian (direct term). Moreover, for a full diagonalization to the eigenvalue problem the scalability goes like the third power of the kmesh. An estimate for the scaling of the screened Coulomb interaction, which enters the direct term of the BSE Hamiltonian, and the full BSE eigenvalue problem is the following
(13)where $N$ stands for the "number of" and the subscripts denote the kpoints, qpoints, Gvectors, valence and conduction states.
The most important convergence parameters are listed here as a supplement to those explained in the tutorial Excited states from TDDFT.

Exercises
 For the given input file, decrease the local fields cutoff gqmax and check how the results change.
 Modify the attribute bsetype to calculate the independentparticle (IP), RPA (exchange only) and triplet absorption spectrum for LiF. How is the spectrum modified? What are the main differences between singlet and triplet spectra? What happens when the direct term of the BSE Hamiltonian is switched off (RPA and IP)? If you have already computed the singlet spectra, you only need to execute the last task of the BSE program flow, i.e., "bse", when switching the bsetype (see workflow).
 Modify the number of conduction bands considered in the BSE window by changing the third and fourth indices in the attributes nstlbse. What is the effect on the spectrum?
4. BSE calculations on top of GW
The procedure to run BSE calculations in exciting illustrated above can be applied after computing the underlying electronic properties from a singleshot GW calculation. For details regarding bandstructure calculations using GW we refer to Electronic bandstructure from GW.
In order to use the quasiparticle energies computed from GW as starting point for the BSE calculation, it is sufficient to include in the input file the gw element. With this option the electronic eigenstates and eigenenergies are read from the GW output. The scissors operator specified by the attribute scissor is automatically ignored. Important: Make sure that the working directory of the BSE calculation contains the GW output file EVALQP.OUT.
5. BSE beyond the TDA
To go beyond the TammDancoff approximation we need to add coupling = "true" in the BSE element of the input file. Compared to the previous calculations on LiF performed within the TDA as described above, now also the resonantantiresonant coupling block of the screened Coulomb interaction needs to be recalculated, since the Fourier coefficients of this block are needed at different kpoints than those of the resonantresonant block. No additional action is required for the exchange interaction. To do this, add to the input file the plan element with related attributes shown below:
... <BSE bsetype="singlet" coupling="true" nstlbse="1 5 1 4" /> ... <plan> <doonly task="screen"/> <doonly task="scrcoulint"/> <doonly task="bse"/> </plan> </xs> ...
To run this calculation, once again type
$ time exciting_smp
After the calculation is finished, new files will appear in the working directory as well as in the subfolders. Again, we are interested in analyzing the new spectra and the new excitation energies. Please note that the new files DO NOT contain the additional string TDABAR. To plot the imaginary part of the dielectric function, follow the following procedure:
 Execute the following commands
$ cp EPSILON/EPSILON_BSEsinglet_SCRfull_OC11.OUT singletFull
 Then, run the script PLOTfiles.py (for a detailed description of the script arguments see The python script "PLOTfiles.py") with the following options:
$ PLOTfiles.py f singletTDA singletFull lx 'Energy [eV]' ly 'Im $\varepsilon_M$' t 'Macroscopic dielectric function' g rc cy 3 x 0 27 lp 2
The resulting graph (PLOT.png) will look like this:
In the case of LiF there is no significant difference between the optical absorption spectrum computed within the TDA and by solving the full BSE.
Literature
 BSE1951: The BSE has been introduced for the first time in E.E. Salpeter and H.A. Bethe, Phys. Rev. 84, 1232 (1951).
 SAD2009: Details on implementation and application of the BSE formalism in exciting: Stephan Sagmeister and Claudia AmbroschDraxl, Timedependent density functional theory versus BetheSalpeter equation: An allelectron study, Phys. Chem. Chem. Phys. 11, 4451 (2009) (PDF).
 More details on the implementation of the BSE formalism within the LAPW method and applications to organic semiconductors: S. Sagmeister, PhD thesis, University of Graz, August 2009 (PDF).
 Additional information on the BSE implementation beyond the TDA and zero momentum transfer can be found in B. Aurich's Master thesis (PDF).
 STR1988: "Application of the Green’s functions method to the study of the optical properties of semiconductors", G. Strinati, Riv. Nuovo Cim. (1988) (PDF).
 SAN2015: "Beyond the TammDancoff approximation for extended systems using exact diagonalization", T. Sander, et al (2015) (PDF).
 ONI2002: "Electronic excitations: densityfunctional versus manybody Green'sfunction approaches" G. Onida, et al (2015) (PDF).