Excited States from BSE

by Weine Olovsson, Stephan Sagmeister, Juan Pablo Echeverry, & Caterina Cocchi for exciting boron

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 that will appear will be 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
      ngridk="10  10  10"
      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.

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. In order to compute the optical absorption spectra of materials, many-body perturbation theory (MBPT) offers a powerful tool such as the Bethe-Salpeter equation (BSE) formalism (see BSE-1951). In solid-state physics, the BSE formalism is used to describe the electron-hole pair (i.e., the exciton), which can be created when the system absorbs a photon. In this case an effective two-particle Hamiltonian is introduced. In the following we use the notation of SAD-2009.

In the BSE formalism, one solves the eigenvalue problem for the effective Hamiltonian $H ^{\rm (eff)}$,

(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$).

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

(3)

\begin{align} H^{\rm (diag)}_{v\,c\,\mathbf{k},\,v'c'\mathbf{k'}} = (\epsilon_{c\,\mathbf{k}} - \epsilon_{v\,\mathbf{k}}) \:\delta_{v\,v'}\:\delta_{c\,c'}\:\delta_{\mathbf{k\,k'}}\:, \end{align}

where $\epsilon_{n\,\mathbf{k}}$ is the single-particle energy of the $n$-th band. Here, the label $v$ ($c$) is used for valence (conduction) states.

The exchange term, which is caused by the repulsive interaction between the electron and the hole, is defined as

(4)

\begin{align} H^{\rm (x)}_{v\,c\,\mathbf{k},\,v'c'\mathbf{k'}} = \int \!\!d^{3}\mathbf{r}\!\!\int\!\! d^{3}\mathbf{r'} \:\varphi_{v\,\mathbf{k}}(\mathbf{r})\: \varphi_{c\,\mathbf{k}}^* (\mathbf{r})\: \overline{v}(\mathbf{r,r'}) \:\varphi^*_{v'\mathbf{k'}}(\mathbf{r'})\: \varphi_{c'\mathbf{k'}} (\mathbf{r'})\:, \end{align}

where the $\varphi_{n\mathbf{k}}(\mathbf{r})$ are the single-particle states and $\overline{v}(\mathbf{r,r'})$ is the bare (unscreened) Coulomb potential which is defined without the long-range term.

Finally, the correlation term (also known as direct term), $H^{\rm (c)}$, which reads

(5)

\begin{align} H^{\rm (c)}_{v\,c\,\mathbf{k},\,v'c'\mathbf{k'}} = -\int\!\! d^{3}\mathbf{r}\!\!\int\!\! d^{3}\mathbf{r'}\: \varphi_{v\,\mathbf{k}}(\mathbf{r})\: \varphi_{c\,\mathbf{k}}^* (\mathbf{r'}) \:W(\mathbf{r,r'}) \:\varphi_{v'\mathbf{k'}}^*(\mathbf{r}) \:\varphi_{c'\mathbf{k'}} (\mathbf{r'})\:, \end{align}

contains the attractive electron-hole correlation, due to the screened Coulomb potential $W(\mathbf{r,r'})$. This screened potential is expressed in reciprocal space as:

(6)

\begin{align} W_{\mathbf{G,G'}}(\mathbf{q}) = \varepsilon^{-1}_{\mathbf{G,G'}}(\mathbf{q})\: {v}_{\mathbf{G'}}(\mathbf{q})\:, \end{align}

where the screening of the Coulomb potential ${v}$ is given by the microscopic dielectric function $\varepsilon^{-1}$.

The optical absorption, which is the main physical quantity obtained from a BSE calculation, is given by the imaginary part of the frequency dependent macroscopic dielectric function,

(7)

\begin{align} \operatorname{Im} \left[\varepsilon_{\!M}(\omega)\right] = \dfrac{8\pi^2}{V}\sum_{\lambda}\: |t_{\lambda}|^2\:\delta(\omega - E_{\lambda}) \end{align}

where $t_{\lambda}$ is related to the momentum matrix elements as well as to the components $A^{\lambda}_{v\,c\,\mathbf{k}} = \langle v\,c\,\mathbf{k}| A_\lambda\rangle$ of the BSE eigenvectors:

(8)

\begin{align} t_{\lambda}=\sum_{v\,c\,\mathbf{k}}\:A^{\lambda}_{v\,c\,\mathbf{k}}\: \dfrac{\langle v\,\mathbf{k}|\hat{\mathbf{p}}|c\,\mathbf{k}\rangle}{\epsilon_{c\,\mathbf{k}} - \epsilon_{v\,\mathbf{k}}} \end{align}

Now, we calculate the macroscopic dielectric function of LiF applying the BSE formalism.

---

3. How to run a BSE calculation

i) Preparation of the input file and run

We now perform the BSE calculation. Move to the parent directory LiF-BSE and create here 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. In this way we can skip the ground state by adding 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"
      tevout="true">
 
      <energywindow 
         intv="0.0 1.0" 
         points="1200"/>
 
      <screening 
         screentype="full"
         nempty="100"/>
 
      <BSE 
         bsetype="singlet"
         nstlbse="1 5 1 4"
         nstlbsemat="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 can be set also in the xs element. It is used in the SCF calculation necessary to obtain 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 to calculate the response function, as explained in more details in Excited states from TDDFT. In the BSE calculation the response function is needed to compute the screening;
• 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 for the calculate spectra;
• When the attribute scissor is set different from zero, a scissor operator is applied to correct the band gap obtained from the ground state calculation (often understimated in LDAor GGA) in order to mimic the quasi-particle gap. This parameter should not appear when the electronic structure resulting from a GW calculation are taken as a starting point for a BSE calculation;
• 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 see 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 the triplets, where the exchange term of the Hamiltonian is neglected, and RPA, where the correlation 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 attributes nstlbse and nstlbsemat define 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 numerated starting from the lowest occupied state, while unoccupied bands are numerated 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 four empty bands.
• The element qpointset is related to the dependence of the response 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.

We are now ready to run the BSE calculation:

$ time excitingser &

Additional information about the calculation work-flow

[Show]

The work-flow for a BSE calculation is shown and described in the following:

1. As a first step, the xs element triggers a one-step ground state calculation to generate a set of Kohn-Sham (KS) eigenvalues $\varepsilon_{n\mathbf{k}}$ and eigenfunctions $\psi_{n\mathbf{k}}$. This ground-state calculation uses a previously obtained KS potential (e.g., from a previous groundstate run). It is important to note that the parameters for this one-step ground-state run are completely independent of those in the groundstate element. In particular, the ngridk and nempty parameters that determine the range of $n\mathbf{k}$ indices are xs%ngridk and xs%nempty , respectively. These parameters are independent, and do not conflict with groundstate%ngridk and groundstate%nempty.
2. Next, the momentum matrix elements $p_{\alpha, n m \mathbf{k}}=\langle \psi_{n\mathbf{k}}|\hat p_{\alpha}|\psi_{m\mathbf{k}+\mathbf{q}} \rangle$ are calculated. The size of this matrix is determined by the number of KS orbitals obtained in the previous step. These matrix elements are used later in the work-flow for calculating the coefficients $t_{\lambda}$ defined in Section 2.
3. Then, the static dielectric matrix in the RPA approximation, $\varepsilon^\textrm{RPA}_{\mathbf{G,G'}}(\mathbf{q})$, is calculated. This calculation involves three important intermediate steps: (i) A new one-step ground state calculation is performed, using the same KS potential as in the first step. The parameters for this additional ground state calculation are independent of those in the xs and groundstate elements. There is however a requirement: screening%nempty $\ge$ xs%nempty. The reason for this is that the calculation of the screening must involve typically many more empty KS orbitals than those determining the electron-hole pairs in the definition of the BSE Hamiltonian, which usually involve a small number of empty bands. (ii) For this newly generated set of KS orbitals, the momentum matrix elements $p_{\alpha, n m \mathbf{k}}$ are also calculated and written into a separate file. (iii) The independent particle density response function $\chi_0$ is assembled using the momentum matrix elements obtained in (ii) and plane-wave matrix elements calculated on-the-fly and not subsequently written into files. Finally, the dielectric matrix is computed as $\varepsilon^{RPA}=\mathbf{1}-v\chi_0$ and written into separate files for each $\mathbf{q}$.
4. Next the direct term of the BSE Hamiltonian is calculated. This is done by first reading $\varepsilon^{RPA}$ from the files generated in the previous step, then inverting it to build the screened Coulomb interaction matrix $W_{\mathbf{G,G'}}$ as defined in Section 2, and finally building the matrix $H^{\rm(c)}_{v\,c\,\mathbf{k},\,v'c'\mathbf{k'}}$. The KS orbitals used to build this matrix are those obtained in the first step, and the band ranges for valence ($v$) and conduction ($c$) band indices are determined by the attribute nstlbsemat. This means that the number of empty states generated in the first step (i.e., xs%nempty) must be equal or larger than the highest index for conduction bands set in the attribute nstlbsemat.
5. Next the exchange term of the BSE Hamiltonian, $H^{\rm (x)}_{v\,c\,\mathbf{k},\,v'c'\mathbf{k'}}$, is obtained. For its calculation only plane-wave matrix elements between KS orbitals calculated in step 1 are used. These matrix elements are calculated on the fly and not stored in any file. The valence and conduction band ranges are defined by the attribute nstlbsemat, as in the previous step.
6. Finally, the BSE effective Hamiltonian $H^{\rm(eff)}_{v\,c\,\mathbf{k},\,v'c'\mathbf{k'}}$ is built and diagonalized. The band ranges are defined by the attribute nstlbse, therefore the ranges defined in nstlbse must be contained in those defined by nstlbsemat, otherwise the program terminates with an error message. This independent definition of the band limits for building the BSE Hamiltonian and for the actual calculation of the matrix elements offers the flexibility of studying the behavior of the dielectric function for different selections of the electron-hole pairs without recalculating the screening (to skip the calculation of the screening set the screening%do attribute equal to skip). Once the BSE effective Hamiltonian is diagonalized, the macroscopic dielectric function is calculated using the expression given in Section 2. The momentum matrix elements used for the evaluation of the coefficients $t_{\lambda}$ are those calculated in step 2.

ii) Analysis of the results

Once the calculation finishes (it will take a few minutes), a xmgrace plot can be generated as follows:

$ cp EPSILON_BSEsinglet_SCRfull_OC11.OUT.xml tmp.xml
$ xsltproc --stringparam filename dielectricfunction $EXCITINGVISUAL/xmldielectricfunction2agr.xsl tmp.xml

Three xmgrace files are generated: dielectricfunction_Re.agr, dielectricfunction_Im.agr, and dielectricfunction_ReKK.agr, which contain, respectively, the real part, the imaginary part, and the real part from a Kramers-Kronig transform of the dielectric function.
The optical absorption spectrum is given by the imaginary part of the dielectric function. To visualize the plot with xmgrace (see also Xmgrace: A Quickstart), type:

$ xmgrace dielectricfunction_Im.agr &

Here is the output plot:

For further analysis of the results we inspect some of the output files (YY being either of the three diagonal components 11, 22, and 33).

• The file EXCITON_BSEsinglet_SCRfull_OCYY.OUT, which contains the relevant information about the excitation energies and intensities. These files contain six columns. After the progressive index referring to the excitation number, we find the excitation energies, exciton binding energies, oscillator strengths, and finally real and imaginary part of the excitation amplitudes, from which the oscillator strengths themselves are computed as the absolute values (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 end of the file and labeled as E_gap. Negative values of the exciton binding energies are associated with bound excitons.
• Inside EXCITON_SORTED_BSEsinglet_SCRfull_OCYY.OUT, the list of excitations is sorted with respect to the corresponding oscillator strengths.
• In the files EPSILON_BSEsinglet_SCRfull_OCYY.OUT, the real and imaginary part of the macroscopic dielectric function are stored in the second and third columns, respectively. In the first column, energies are indicated. In the fourth column, one can find the real part of the macroscopic dielectric function computed from the imaginary part through the Kramers-Kronig relations.

iii) Scaling and convergence

BSE calculations are extremely demanding, even on our state-of-the-art computer infrastructure. The scaling with respect to the size of the k-point mesh is very bad and enters quadratically the setup of the Hamlitonian (direct term), and, if a full diagonalization to the eigenvalue problem is applied, even with the third power. An estimate for the scaling of the screening, the direct term of the BSE Hamiltonian, and the BSE eigenvalue problem (which are the most time-consuming parts of the calculation in general) is the following

(9)

\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 states, 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, nstlbsemat 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)?

[Show spectra]

• Modify the number of conduction bands considered in the BSE window by changing the third and forth indices in the attributes nstlbse and nstlbsemat. What is the effect on the spectrum?

---

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).
• Tutorial talk by Stephan Sagmeister (PDF).
• 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).