Excited States from BSE

by Benjamin Aurich, Christian Vorwerk, & Caterina Cocchi for exciting nitrogen

Purpose: In this tutorial you will learn how to perform a basic Bethe-Salpeter-equation (BSE) calculation. As an example, the optical spectrum of lithium fluoride (LiF) will be studied.

1. Preliminary step: Ground-state 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 excited-state properties from BSE, a ground-state calculation has to be performed. In this tutorial we consider as an example LiF. Create a directory named LiF-BSE and move into it.

$ mkdir LiF-BSE
$ cd LiF-BSE

Inside the directory LiF-BSE we create a sub-directory GS where we perform the preliminary ground state calculation:

$ mkdir GS
$ cd GS

Inside the GS sub-directory we create the input file for LiF. In the structure element we include the lattice parameter and basis vectors of LiF, which has a rock-salt cubic lattice, as well as the positions of the Li and F atoms. In the groundstate element, we include a 10$\times$10$\times$10 k-point mesh (ngridk) and a value of 14.0 for gmaxvr. This value, which is larger than the default, is needed in view of the excited-state calculation.

The resulting input file is the following:

<input>
 
   <title>LiF-BSE</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="10  10  10"
      xctype="GGA_PBE_SOL"
      gmaxvr="14.0"/>
 
</input>

N. B.: Do not forget to replace in the input.xml the string "$EXCITINGROOT" by the actual value of the environment variable $EXCITINGROOT.

$ SETUP-excitingroot.sh

Start now the ground-state SCF calculation and check if it finishes gracefully. 
$ time excitingser &

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 a starting point for the BSE calculation.

2. Theoretical background: Bethe-Salpeter equation

Before starting with the calculation, we briefly review the theoretical background. Many-body perturbation theory (MBPT) offers a powerful tool such as the Bethe-Salpeter equation (BSE) formalism (see BSE-1951) to compute optical absorption spectra of materials. In solid-state physics, the BSE formalism is used to describe the excitations of the system in terms of interacting electron-hole pairs. In particular, the BSE formalism allows to resolve excitations below the band gap. These bound electron-hole pairs are quasi-particles also called excitons, which can be created when the system absorbs a photon.

The BSE Hamiltonian

The BSE can be stated as an eigenvalue problem of an effective two-particle Hamiltonian (see, e.g., SAD-2009):

(1)
\begin{align} H ^{\rm (eff)}\: | A_{\lambda} \rangle= E_{\lambda}\: |A_{\lambda}\rangle\:, \end{align}

where the effective Hamiltonian is defined as

(2)
\begin{align} H ^{\rm (eff)} = H ^{\rm (diag)}+2\:\gamma_{\rm x}\:H ^{\rm (x)} + \gamma_{\rm c}\:H ^{\rm (c)}\:. \end{align}

In this expression, three terms can be recognized, each carrying a specific physical meaning which will be clarified in the following. Moreover, the factors $\gamma_{\rm x}$ and $\gamma_{\rm c}$ allow one to choose different levels of approximation and to distinguish between spin-singlet (i.e., $\gamma_{\rm x} \equiv 1$ and $\gamma_{\rm c} \equiv 1$) and spin-triplet channels (i.e., $\gamma_{\rm x} \equiv 0$ and $\gamma_{\rm c} \equiv 1$).

Since the BSE describes interacting electron-hole pairs, it is natural to introduce a two-particle basis (see, e.g., SAN-2015) which reflects the excitation (resonant terms, labeled r) and the de-excitation (anti-resonant terms, labeled a) of independent electron-hole pairs:

(3)
\begin{align} \Upsilon^\text{r}_{\alpha,\mathbf{q}}(\mathbf{r},\mathbf{r}') = \psi_{v_\alpha \mathbf{k}^+_\alpha}\!(\mathbf{r})\, \psi^*_{c_\alpha \mathbf{k}^-_\alpha}\!(\mathbf{r}')\;\phantom{!} \\[4 mm] \Upsilon^\text{a}_{\alpha,\mathbf{q}}(\mathbf{r},\mathbf{r}') = \psi_{c_\alpha \mathbf{k}^+_\alpha}\!(\mathbf{r})\, \psi^*_{v_\alpha \mathbf{k}^-_\alpha}\!(\mathbf{r}')\;. \end{align}

Here, $\psi_{n\mathbf{k}}\!(\mathbf{r})$ are quasi-particle states. The combined index $\alpha$ incorporates $v$, indicating valence states, $c$ indicating conduction states and $\mathbf{k}$ standing for the k-points. The (crystal)-momentum difference between electron and hole is denoted by $\mathbf{q}$, such that $\mathbf{k}^\pm = \mathbf{k} \pm \mathbf{q}/2$. In the case of optical excitations, no momentum transfer from photon to electron is considered ($\mathbf{q} = 0$). In the basis presented above, the Hamiltonian has a $2 \times 2$ block structure

(4)
\begin{align} H^{\text{(eff)}} = \begin{pmatrix} H^\text{rr} & H^\text{ra} \\ H^\text{ar} & H^\text{aa} \end{pmatrix}, \end{align}

with the eigenvectors split into a resonant and an anti-resonant part:

(5)
\begin{align} H^{\text{(eff)}} \begin{pmatrix} X_\lambda \\ Y_\lambda \end{pmatrix} = E_\lambda \begin{pmatrix} X_\lambda \\ Y_\lambda \end{pmatrix}. \end{align}

We can formally write the eigenvector coefficients in the Dirac notation as

(6)
\begin{align} (X_\lambda)_\alpha &= \langle (v_\alpha \mathbf{k}^+_\alpha) (c_\alpha \mathbf{k}^-_\alpha)^* | A_\lambda \rangle \\[4 mm] (Y_\lambda)_\alpha &= \langle (c_\alpha \mathbf{k}^+_\alpha) (v_\alpha \mathbf{k}^-_\alpha)^* | A_\lambda \rangle. \end{align}

In the BSE formalism it is common to adopt the Tamm-Dancoff-approximation (TDA) which neglects the coupling between the resonant and anti-resonant components of the Hamiltonian. For $\mathbf{q} = 0$, this approximation entails that only the $H^\text{rr}$ part has to be considered. In this way, the necessary computational effort is significantly reduced. In the following, we focus only on the resonant-resonant block of the BSE and discuss its terms.

The diagonal term of the Hamiltonian, $H^{\rm (diag)}$, accounts for the contribution of independent particle transitions energies

(7)
\begin{align} H^{\text{(diag)}}_{\alpha\,\beta}\!(\mathbf{q}) = (\epsilon_{c_\alpha\,\mathbf{k}^-_\alpha} - \epsilon_{v_\alpha\,\mathbf{k}^+_\alpha}) \,\delta_{\alpha\,\beta}, \end{align}

where $\epsilon_{n\,\mathbf{k}}$ is the single-particle energy of the $n$-th band.

The exchange term, is defined as

(8)
\begin{align} H^{\text{(x)}}_{\alpha\,\beta}(\mathbf{q}) = \int \!\!d^{3}\mathbf{r}\!\!\int\!\! d^{3}\mathbf{r'} \,\Upsilon^{\text{r}*}_{\alpha,\mathbf{q}}(\mathbf{r},\mathbf{r}) \,v(\mathbf{r,r'}) \,\Upsilon^\text{r}_{\beta,\mathbf{q}}(\mathbf{r}',\mathbf{r}') \end{align}

and describes the repulsive exchange interaction between the electron-hole pairs $\alpha$ and $\beta$. The pair $\beta$ annihilates at $\mathbf{r}'$ and the pair $\alpha$ is created at $\mathbf{r}$. The interaction is mediated by the bare Coulomb potential $v(\mathbf{r,r'})$. For $\mathbf{q}=0$, only the short-range part of this term is included (see, e.g., SAD-2009, ONI-2002).

Finally, the direct term $H^{\text{(c)}}$ reads

(9)
\begin{align} H^{\text{(c)}}_{\alpha\,\beta}(\mathbf{q}) = -\int\!\! d^{3}\mathbf{r}\!\!\int\!\! d^{3}\mathbf{r'} \,\Upsilon^{\text{r}*}_{\alpha,\mathbf{q}}(\mathbf{r},\mathbf{r}') \,W(\mathbf{r,r'}) \,\Upsilon^\text{r}_{\beta,\mathbf{q}}(\mathbf{r},\mathbf{r}')\;. \end{align}

It describes the scattering of the electron-hole pair $\beta$ into the pair $\alpha$ due to the screened Coulomb potential $W(\mathbf{r,r'})$. It accounts for the attractive electron-hole interaction including the screened potential of the form:

(10)
\begin{align} W(\mathbf{r}, \mathbf{r}') = \int\!\! d^{3}\mathbf{r}'' \,v(\mathbf{r}, \mathbf{r}'')\, \varepsilon^{-1}\!(\mathbf{r}'',\mathbf{r}') \;, \end{align}

where the screening of the Coulomb potential ${v}$ is given by the microscopic dielectric function $\varepsilon^{-1}$, which is calculated here within the random-phase approximation (RPA).

In the limit of vanishing $\mathbf{q}$ (optical limit), one can obtain the macroscopic transverse dielectric tensor $\varepsilon^{ij}_{\text{M}}(\omega)$ from the BSE solutions. The imaginary part of this tensor yields the optical absorption spectrum, as expressed by

(11)
\begin{align} \operatorname{Im} \left[\varepsilon^{ij}_{\text{M}}(\omega)\right] = \dfrac{4\pi^2}{V} \sum_{\lambda}\: t^*_{\lambda,i} t_{\lambda,j} \,\frac{1}{\pi \delta} \left( \frac{\delta^2}{(\omega-E_\lambda)^2+\delta^2} - \frac{\delta^2}{(\omega+E_\lambda)^2+\delta^2} \right)\;. \end{align}

Thus, considering a diagonal element, the spectrum is constructed from scaled positive (negative) Lorentizian function of width $\delta$ centered around the positive (negative) BSE eingenenergies. The scaling factors are determined by the transition coefficients $t_{\lambda, i}$ which are a sum of dipole-transition matrix elements weighted by the BSE eigenvectors:

(12)
\begin{align} t^*_{\lambda, i}=\mathrm{i} \sum_{\alpha} \left( X^{\lambda}_{\alpha} + Y^{\lambda}_{\alpha} \right) \langle c_\alpha\,\mathbf{k}_\alpha|-\hat{r}_i|v_\alpha\,\mathbf{k}_\alpha\rangle\,, \end{align}

where, within the TDA, $Y^\lambda$ = 0.

In the actual computation, the dipole matrix elements are replaced by momentum matrix elements obtained by using the corresponding canonical commutator relation.

3. How to run a BSE calculation

i) Preparation of the input file and run

We start by performing a BSE calculation using the TDA and at $\mathbf{q}=0$. To do this, move to the parent directory LiF-BSE and create a new folder which we name BSE. Then, we move into it:

$ cd ..
$ pwd
/home/tutorials/LiF-BSE
$ mkdir BSE
$ cd BSE

As anticipated, we need the files STATE.OUT and EFERMI.OUT obtained from the ground state calculation to be inside the new directory. Copy these two files from the GS folder into the BSE one.

$ cp ../GS/{STATE.OUT,EFERMI.OUT,input.xml} ./

In this way we can skip the ground-state calculation by setting do = "skip" in the groundstate element of the input file, also copied from the GS folder.

To perform an excited-state calculation we must include the xs element in the input file:

...
   <xs 
      xstype="BSE" 
      ngridk="4 4 4" 
      vkloff="0.097 0.273 0.493"
      ngridq="4 4 4"
      nempty="30"
      gqmax="3.0"
      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>                 
...

Several attributes appear in this section of the input file. We will go through all of them in order to clarify their meaning within the BSE calculation we are about to perform:

  • The attribute xstype defines the type of excited-state calculation: In this case it is set to BSE;
  • The attribute ngridk is also set in the xs element and it is used in the (single-shot) SCF calculation which is necessary to obtain the single-particle eigenstates and eigenenergies. Usually the same mesh as in the groundstate is taken or, as in this case, a coarser grid is considered to reduce the computational costs. If not defined, the k-point mesh used in the groundstate calculation will be assumed;
  • The attribute vkloff is used to shift the k-mesh off symmetry by a small displacement, in order to break all symmetry relations among the k-points of the mesh. In this way, all the k-points in the mesh will be crystallographically inequivalent and there will be no redundant contribution to the spectrum;
  • The attribute ngridq defines the q-mesh for the calculation of the dynamical screening. It is a good practice to choose a q-mesh equivalent to the k-mesh;
  • The attribute nempty determines the number of empty states available for the construction of the BSE Hamiltonian;
  • The attribute gqmax is an energy threshold for the local field effects included in the calculation. The number of G-vectors employed in the calculation is written in the first line of the file GQPOINTS_QM001.OUT. Tuning gqmax may require to adjust the value of gmaxvr in the groundstate element;
  • The attribute 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 quasi-particle gap. This parameter should not appear when the electronic structure resulting from a GW calculation is taken as starting point for the BSE calculation;
  • The attribute tappinfo allows to append the output printed in the log file INFOXS.OUT at each run, without overwriting the existing file.
  • The attribute tevout sets the energy output in electronvolt (eV).

Additional sub-elements appear inside the xs element (refer to the Input Reference for further details):

  • The element 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 points to be used.
  • screening is a rather important element in a BSE calculation, as it defines the parameters used to compute the screened Coulomb potential. The attribute screeningtype sets the adopted 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 be always 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 triplets, where the exchange term of the Hamiltonian is neglected, RPA, where the direct term is ignored, and the independent-particle (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 q-dependent 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$).

We are now ready to run the BSE calculation:

$ time excitingser &

While the calculation is running, some strings will be printed on screen. A successfully completed BSE run will produce the following lines:

Calculating RPA Dielectric Function:                100.000 
Calculating Screened Coulomb Potential:             100.000 
Calculating Screened Coulomb Matrix Elements:       100.000 
Calculating Plane-wave matrix elements:             100.000 
Calculating Exchange Interaction Matrix Elements:   100.000 
Solving BSE Eigenvalue Problem:                     100.000
Additional information about the calculation work-flow
[Show]
ii) Analysis of the results

After running the calculation, 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. Since we are interested in plotting the optical absorption spectrum of LiF, given by the imaginary part of its macroscopic dielectric function, we consider the folder EPSILON and move into it. To do so and subsequently plot the spectrum, execute the following commands:

$ cd EPSILON
$ cp EPSILON_BSE-singlet-TDA-BAR_SCR-full_OC11.OUT.xml tmp.xml
$ xsltproc --stringparam filename dielectricfunction $EXCITINGVISUAL/xmldielectricfunction2agr.xsl tmp.xml

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). With the executions of the commands above, three xmgrace files are generated: dielectricfunction_Re.agr, dielectricfunction_Im.agr, and dielectricfunction_ReKK.agr, which contain the real part, the imaginary part, and the real part from a Kramers-Kronig transform of the macroscopic dielectric function, respectively. This information is contained inside the EPSILON_BSE-singlet-TDA-BAR_SCR-full_OC11.OUT file in the second, third, and fourth column, respectively. To visualize the plot of the imaginary part with xmgrace (see also Xmgrace: A Quickstart), type:

$ xmgrace dielectricfunction_Im.agr >/dev/null 2>&1 &
Here is the resulting output plot:
dielectricfunction_Im.png

For further analysis of the results we inspect also the output files contained in the folder EXCITON. This directory contains the solutions of the diagonalization of the BSE. Specifically, the files EXCITON_BSE-singlet-TDA-BAR_SCRfull_OCXX.OUT, where XX=11,22,33 include the relevant information about the excitation energies and intensities. Each file has six columns. After the progressive index referring to the excitation number, we find the excitation energies, exciton binding energies, oscillator strengths, and finally the real and the imaginary part of the transition coefficients from which the oscillator strengths are computed (see also SAD-2009). A relevant quantity is the exciton binding energy, which is defined as the difference between the excitation energy and the band gap, which is printed at the top of the file and labeled as E_shift. Negative values of the exciton binding energies are associated with bound excitons.

Other folders that contain useful quantities for the physical interpretation of the results are SIGMA, including the optical conductivity, and LOSS with the loss functions. The latter will be discussed more in details in the tutorial on q-dependent BSE.

iii) Scaling and convergence

BSE calculations are extremely demanding, even on state-of-the-art computer infrastructure. The scaling with respect to the size of the k-mesh 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 k-mesh. 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)
\begin{align} T_{\rm BSE} \sim \alpha_{\rm SCR}(N_v N_c N_{\bf k}) N_{\bf q}N_{\bf G}^2 + \alpha_{\rm HAM}(N_v N_c N_{\bf k})^2 N_{\bf G}^2 + \alpha_{\rm DIAG} (N_v N_c N_{\bf k})^3 \nonumber \end{align}

where $N$ stands for the "number of" and the subscripts denote the k-points, q-points, G-vectors, 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.

Description Parameter Info
k-mesh/q-mesh ngridk, ngridq inside xs crucial parameter but sometimes hard to be converged, due to the computational effort
energy cutoff (screening) nempty inside screening to be converged
local field effects G-cutoff gqmax inside xs defines the quality of the screened Coulomb potential in G-space and should be converged. High values of gqmax can be computationally very costly.
number of states included in the BSE window nstlbse inside BSE it is necessary to include more states in this window in order to explore higher energy portion of the optical spectrum. Also, adding more states to this window can induce transitions mixing, which in turn may affect the spectral features.
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 independent-particle (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 work-flow).
[Show spectra]
  • Modify the number of conduction bands considered in the BSE window by changing the third and forth 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 single-shot GW calculation. For details regarding band-structure calculations using GW we refer to Electronic bandstructure from GW.

In order to use the quasi-particle 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 Tamm-Dancoff 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 resonant-anti-resonant coupling block of the screened Coulomb interaction needs to be recalculated, since the Fourier coefficients of this block are needed at different k-points than those of the resonant-resonant block. No additional action is instead required for the exchange interaction. To do this, add to the input file 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 excitingser &

After the calculation is finished (it will take a few minutes) new files will appear in the working directory as well as in the sub-folders. 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 -TDA-BAR. To plot the imaginary part of the dielectric function, follow the same procedure illustrated above. To simply overlay the new spectrum on top of the TDA one produced before, you may simply plot the first and third column of the file EPSILON_BSE-singlet_SCR-full_OC11.OUT. The resulting graph will look like this:

dielectricfunction_Im_tdafull.png

It can be noticed that 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

  • BSE-1951: The BSE has been introduced for the first time in E.E. Salpeter and H.A. Bethe, Phys. Rev. 84, 1232 (1951).
  • SAD-2009: Details on implementation and application of the BSE formalism in exciting: Stephan Sagmeister and Claudia Ambrosch-Draxl, Time-dependent density functional theory versus Bethe-Salpeter equation: An all-electron 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 the excerpt from the Master's thesis of B. Aurich here.
  • STR-1988: "Application of the Green’s functions method to the study of the optical properties of semiconductors", G. Strinati, Riv. Nuovo Cim. (1988) (PDF).
  • SAN-2015: "Beyond the Tamm-Dancoff approximation for extended systems using exact diagonalization", T. Sander, et al (2015) (PDF).
  • ONI-2002: "Electronic excitations: density-functional versus many-body Green's-function approaches" G. Onida, et al (2015) (PDF).
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License