Excited States from BSE

by Benjamin Aurich, Christian Vorwerk, Caterina Cocchi, & Keith Gilmore for exciting sodium

(Jupyter notebook by Maximilian Schebek and Mara Voiculescu)

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.

Table of Contents

0. Before Starting

1. Theoretical Background: Bethe-Salpeter Equation 2. Preliminary Step: Ground-State Calculation

3. How to Run a BSE Calculation

4. BSE Calculations on Top of GW

5. BSE Beyond the TDA

Literature

0. Before Starting

Read the following paragraphs before starting with the rest of this tutorial!

Before running any Jupyter tutorials, please refer to the 00_before_starting.md document on how to correctly set up the environment. This only needs to be done once. After which, the venv can be (re)activated from exciting's root directory:

source $EXCITINGROOT/tools/excitingjupyter/venv/excitingvenv/bin/activate

As a first step, you may create a running directory for the notebook.

In [1]:
%%bash
mkdir -p run_LiF_BSE

1. Theoretical Background: Bethe-Salpeter Equation

$\Rightarrow$ Click to show theoretical background!

Before starting with the calculation, we briefly review the theoretical background. Many-body perturbation theory (MBPT) offers powerful tools such as the Bethe-Salpeter equation (BSE) formalism (see BSE-1951) to compute bound two-particle states. In the context of optical absorption, the BSE formalism is used to describe electron-hole pair excitations. Of particular interest are bound electron-hole pairs with energies below the band gap, which are quasi-particles referred to as excitons.

The BSE Hamiltonian

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

$$ H ^{\rm (eff)}\: | A_{\lambda} \rangle= E_{\lambda}\: |A_{\lambda}\rangle\:, $$

where the effective Hamiltonian is defined as

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

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:

$$ \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}')\;. $$

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

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

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

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

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

$$ \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

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

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

The exchange term, is defined as

$$ 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}') $$

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

$$ 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}')\;. $$

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:

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

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

$$ \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)\;. $$

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:

$$ 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\,, $$

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.

2. Preliminary Step: Ground-State Calculation

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.

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

In [2]:
%%bash
cd run_LiF_BSE
mkdir -p GS
cd ..

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 8×8×8 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.

<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="8 8 8"
      xctype="GGA_PBE_SOL"
      radialgridtype="cubic-2"
      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

In [4]:
%%bash
cd run_LiF_BSE/GS
python3 -m excitingscripts.setup.excitingroot
cd ../..

Now start the ground-state SCF calculation.

In [ ]:
%%bash
cd run_LiF_BSE/GS
time $EXCITINGROOT/install/bin/exciting_smp input.xml
cd ../..

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, create a new folder in the running directory which we name BSE. Then, copy the necessary ground-state files from the GS folder into the BSE directory.

In [6]:
%%bash
cd run_LiF_BSE
mkdir -p BSE && cd BSE
cp ../GS/{STATE.OUT,EFERMI.OUT,input.xml} ./
cd ../..

We can then skip the ground-state calculation by setting do="skip" in the groundstate element of the input file.

To perform an excited-state 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 excited-state calculation: here we set it to BSE;
  • ngridk within the xs element is used in the non-self-consistent-field calculation required to obtain the single-particle eigenstates and eigenenergies. If not defined, it defaults to the k-point mesh used in the groundstate calculation. Calculations should be converged with respect to the k-mesh;
  • 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;
  • ngridq defines the q-mesh for the calculation of the screening. It is a good practice to choose a q-mesh equivalent to the k-mesh;
  • 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 G-vectors 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 re-running the ground-state 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 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;
  • tappinfo causes output to be appended to the log file INFOXS.OUT without overwriting the existing file.
  • dtevouto sets the energy output to electronvolt (eV).

Additional sub-elements 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 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 tutorial 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$).

ii) Additional Information about the Calculation Workflow

$\Rightarrow$ Click to show workflow!

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

No description has been provided for this image

  1. As a first step, the xs element triggers a one-step groundstate 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 (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_{i, n m \mathbf{k}}=\langle \psi_{n\mathbf{k}}|\hat p_{i}|\psi_{m\mathbf{k}} \rangle$ are calculated. The size of this matrix is determined by the number of KS orbitals obtained in the previous step. The matrix elements are used later in the work-flow for calculating the coefficients $t_{\lambda,i}$ defined in Section 1.

  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_{i, 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 Fourier transform of the screened Coulomb interaction matrix $W_{\mathbf{G,G'}}(\mathbf{q})$, and finally building the matrix $H^{\rm(c)}_{\alpha\,\beta}$. 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 nstlbse. 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 nstlbse.

  5. Next the exchange term of the BSE Hamiltonian, $H^{\rm (x)}_{\alpha\,\beta}$, 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 nstlbse, as in the previous step.

  6. Finally, BSE effective Hamiltonian $H^{\rm(eff)}_{\alpha\,\beta}$ is built and diagonalized. The band ranges are defined by the attribute nstlbse. Once the BSE effective Hamiltonian is diagonalized, the macroscopic dielectric function (and other derived quantities) is calculated using the expression given in Section 1. The momentum matrix elements used for the evaluation of the coefficients $t_{\lambda,i}$ are those calculated in step 2.

The execution of these steps can be controlled directly from the input file by adding the plan element to the xs block and specifying which task should be run. The standard BSE execution follows implicitly the plan reported below:

...
   <xs...>
 
      <plan>
         <doonly task="xsgeneigvec"/>
         <doonly task="writepmatxs"/>
         <doonly task="scrgeneigvec"/>
         <doonly task="scrwritepmat"/>
         <doonly task="screen"/>
         <doonly task="scrcoulint"/>
         <doonly task="exccoulint"/>
         <doonly task="bse"/>
      </plan>

   </xs>                 
...

iii) Running exciting and analysing the results

We can now run the BSE calculation:

In [ ]:
%%bash
cd run_LiF_BSE/BSE
time $EXCITINGROOT/install/bin/exciting_smp input.xml
cd ../..

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:

In [9]:
%%bash
cd run_LiF_BSE/BSE
cp EPSILON/EPSILON_BSE-singlet-TDA-BAR_SCR-full_OC11.OUT singlet-TDA
python3 -m excitingscripts.plot.files -f singlet-TDA  -lx 'Energy [eV]'  -ly 'Im $\varepsilon_M$'  -t 'Macroscopic dielectric function'  -g  -rc  -cy 3  -x 0 27  -nl
cd ../..

The result is stored in PLOT.png and should resemble this plot:

No description has been provided for this image

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_BSE-singlet-TDA-BAR_SCR-full_OC11.OUT in the directory EPSILON also stores the real part of the macroscopic dielectric function calculated directly and via a Kramers-Kronig 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 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

$$ 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 $$

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 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 transition 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 workflow in Section 3).

No description has been provided for this image
  • 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 single-shot GW calculation. For details regarding band-structure calculations using GW we refer to the tutorial 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 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

In [ ]:
%%bash
cd run_LiF_BSE/BSE
time $EXCITINGROOT/install/bin/exciting_smp input.xml
cd ../..

After the calculation is finished, 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 following procedure:

In [12]:
%%bash
cd run_LiF_BSE/BSE
cp EPSILON/EPSILON_BSE-singlet_SCR-full_OC11.OUT singlet-Full
python3 -m excitingscripts.plot.files -f singlet-TDA singlet-Full  -lx 'Energy [eV]'  -ly 'Im $\varepsilon_M$'  -t 'Macroscopic dielectric function'  -g  -rc  -cy 3  -x 0 27  -lp 2
cd ../..

The resulting graph (PLOT.png) will look like this:

No description has been provided for this image

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 B. Aurich's Master thesis PDF.
  • 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.
excited_states_from_bse

Excited States from BSE

by Benjamin Aurich, Christian Vorwerk, Caterina Cocchi, Keith Gilmore, & Benedikt Maurer for exciting sodium

(Jupyter notebook by Maximilian Schebek and Mara Voiculescu)

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.

Table of Contents

0. Before Starting

1. Theoretical Background: Bethe-Salpeter Equation 2. Preliminary Step: Ground-State Calculation

3. How to Run a BSE Calculation

4. BSE Calculations on Top of GW

5. BSE Beyond the TDA

6. Execute Parallel BSE Calculations

Literature

0. Before Starting

Read the following paragraphs before starting with the rest of this tutorial!

Before running any Jupyter tutorials, please refer to the 00_before_starting.md document on how to correctly set up the environment. This only needs to be done once. After which, the venv can be (re)activated from exciting's root directory:

source $EXCITINGROOT/tools/excitingjupyter/venv/excitingvenv/bin/activate

As a first step, you may create a running directory for the notebook.

In [ ]:
%%bash
mkdir -p run_LiF_BSE

1. Theoretical Background: Bethe-Salpeter Equation

$\Rightarrow$ Click to show theoretical background!

Before starting with the calculation, we briefly review the theoretical background. Many-body perturbation theory (MBPT) offers powerful tools such as the Bethe-Salpeter equation (BSE) formalism (see BSE-1951) to compute bound two-particle states. In the context of optical absorption, the BSE formalism is used to describe electron-hole pair excitations. Of particular interest are bound electron-hole pairs with energies below the band gap, which are quasi-particles referred to as excitons.

The BSE Hamiltonian

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

$$ H ^{\rm (eff)}\: | A_{\lambda} \rangle= E_{\lambda}\: |A_{\lambda}\rangle\:, $$

where the effective Hamiltonian is defined as

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

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:

$$ \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}')\;. $$

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

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

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

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

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

$$ \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

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

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

The exchange term, is defined as

$$ 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}') $$

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

$$ 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}')\;. $$

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:

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

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

$$ \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)\;. $$

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:

$$ 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\,, $$

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.

2. Preliminary Step: Ground-State Calculation

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.

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

In [ ]:
%%bash
cd run_LiF_BSE
mkdir -p GS
cd ..

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 8×8×8 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.

<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="8 8 8"
      xctype="GGA_PBE_SOL"
      radialgridtype="cubic-2"
      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

In [ ]:
%%bash
cd run_LiF_BSE/GS
python3 -m excitingscripts.setup.excitingroot
cd ../..

Now start the ground-state SCF calculation.

In [ ]:
%%bash
cd run_LiF_BSE/GS
time $EXCITINGROOT/install/bin/exciting_smp input.xml
cd ../..

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, create a new folder in the running directory which we name BSE. Then, copy the necessary ground-state files from the GS folder into the BSE directory.

In [ ]:
%%bash
cd run_LiF_BSE
mkdir -p BSE && cd BSE
cp ../GS/{STATE.OUT,EFERMI.OUT,input.xml} ./
cd ../..

We can then skip the ground-state calculation by setting do="skip" in the groundstate element of the input file.

To perform an excited-state 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 excited-state calculation: here we set it to BSE;
  • ngridk within the xs element is used in the non-self-consistent-field calculation required to obtain the single-particle eigenstates and eigenenergies. If not defined, it defaults to the k-point mesh used in the groundstate calculation. Calculations should be converged with respect to the k-mesh;
  • 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;
  • ngridq defines the q-mesh for the calculation of the screening. It is a good practice to choose a q-mesh equivalent to the k-mesh;
  • 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 G-vectors 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 re-running the ground-state 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 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;
  • tappinfo causes output to be appended to the log file INFOXS.OUT without overwriting the existing file.
  • dtevouto sets the energy output to electronvolt (eV).

Additional sub-elements 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 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 tutorial 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$).

ii) Additional Information about the Calculation Workflow

$\Rightarrow$ Click to show workflow!

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

No description has been provided for this image

  1. As a first step, the xs element triggers a one-step groundstate 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 (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_{i, n m \mathbf{k}}=\langle \psi_{n\mathbf{k}}|\hat p_{i}|\psi_{m\mathbf{k}} \rangle$ are calculated. The size of this matrix is determined by the number of KS orbitals obtained in the previous step. The matrix elements are used later in the work-flow for calculating the coefficients $t_{\lambda,i}$ defined in Section 1.

  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_{i, 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 Fourier transform of the screened Coulomb interaction matrix $W_{\mathbf{G,G'}}(\mathbf{q})$, and finally building the matrix $H^{\rm(c)}_{\alpha\,\beta}$. 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 nstlbse. 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 nstlbse.

  5. Next the exchange term of the BSE Hamiltonian, $H^{\rm (x)}_{\alpha\,\beta}$, 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 nstlbse, as in the previous step.

  6. Finally, BSE effective Hamiltonian $H^{\rm(eff)}_{\alpha\,\beta}$ is built and diagonalized. The band ranges are defined by the attribute nstlbse. Once the BSE effective Hamiltonian is diagonalized, the macroscopic dielectric function (and other derived quantities) is calculated using the expression given in Section 1. The momentum matrix elements used for the evaluation of the coefficients $t_{\lambda,i}$ are those calculated in step 2.

The execution of these steps can be controlled directly from the input file by adding the plan element to the xs block and specifying which task should be run. The standard BSE execution follows implicitly the plan reported below:

...
   <xs...>
 
      <plan>
         <doonly task="xsgeneigvec"/>
         <doonly task="writepmatxs"/>
         <doonly task="scrgeneigvec"/>
         <doonly task="scrwritepmat"/>
         <doonly task="screen"/>
         <doonly task="scrcoulint"/>
         <doonly task="exccoulint"/>
         <doonly task="bse"/>
      </plan>

   </xs>                 
...

iii) Running exciting and analysing the results

We can now run the BSE calculation:

In [ ]:
%%bash
cd run_LiF_BSE/BSE
time $EXCITINGROOT/install/bin/exciting_smp input.xml
cd ../..

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:

In [ ]:
%%bash
cd run_LiF_BSE/BSE
cp EPSILON/EPSILON_BSE-singlet-TDA-BAR_SCR-full_OC11.OUT singlet-TDA
python3 -m excitingscripts.plot.files -f singlet-TDA  -lx 'Energy [eV]'  -ly 'Im $\varepsilon_M$'  -t 'Macroscopic dielectric function'  -g  -rc  -cy 3  -x 0 27  -nl
cd ../..

The result is stored in PLOT.png and should resemble this plot:

No description has been provided for this image

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_BSE-singlet-TDA-BAR_SCR-full_OC11.OUT in the directory EPSILON also stores the real part of the macroscopic dielectric function calculated directly and via a Kramers-Kronig 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 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

$$ 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 $$

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 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 transition 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 workflow in Section 3).

No description has been provided for this image
  • 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 single-shot GW calculation. For details regarding band-structure calculations using GW we refer to the tutorial 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 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

In [ ]:
%%bash
cd run_LiF_BSE/BSE
time $EXCITINGROOT/install/bin/exciting_smp input.xml
cd ../..

After the calculation is finished, 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 following procedure:

In [ ]:
%%bash
cd run_LiF_BSE/BSE
cp EPSILON/EPSILON_BSE-singlet_SCR-full_OC11.OUT singlet-Full
python3 -m excitingscripts.plot.files -f singlet-TDA singlet-Full  -lx 'Energy [eV]'  -ly 'Im $\varepsilon_M$'  -t 'Macroscopic dielectric function'  -g  -rc  -cy 3  -x 0 27  -lp 2
cd ../..

The resulting graph (PLOT.png) will look like this:

No description has been provided for this image

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.

6. Execute Parallel BSE Calculations

Before starting

In the following, we explain how to execute parallel BSE calculations with $\texttt{exciting}$. This requires building $\texttt{exciting}$ with MPI and ScaLAPACK on a machine with multiple cores. Also you should be familiar with the basic concepts of Message Passing Interface (MPI) parallelization and how to start an application on multiple processes.

Instead of giving a step-by-step tutorial that you can follow as before, we now provide some details on the implementation. Based on these we explain how to choose the number of processes for the different parts of the workflow. The following steps can be applied to any $\texttt{exciting}$ BSE calculation. However, to benefit from parallel execution, the system under study must be large enough. The easiest way to increase the problem size is by increasing the number of $\mathbf{k}$ and $\mathbf{q}$ points via the ngridk and ngridq attributes in the xs element. Please note that for BSE calculations, ngridk and ngridq must be the same.

Background on the Implementation

Solving the BSE in $\texttt{exciting}$ consists of three parts. First, one performs a ground-state one-shot runs to set up the dielectric matrix and the BSE Hamiltonian (BSH), along with the momentum matrix elements. Second, one computes the dielectric matrix, which is required for constructing the screened Coulomb kernel that describes the attractive electron-hole interaction. Third, one solves the BSE by setting up the BSH and computing its eigenvalues, and optionally eigenvectors. These can then be used to calculate, for instance, the absorption spectrum.

Calculating the ground-state one-shot runs involves solving the Kohn–Sham equations once using the ground-state density. Therefore, this step — as well as the calculation of momentum matrix elements — can be done independently for each $\mathbf{k}$-point. extbf{Note that this step is not computationally intensive compared to the others and does not represent a bottleneck.}

Indeed, the most computationally intensive steps are computing the dielectric matrix and solving the BSE, the latter of which becomes the main bottleneck with increasing $\mathbf{k}$-grid size.

These steps are parallelized differently across MPI tasks. To optimize performance on an HPC system, it is necessary to understand this distribution. For that, we examine how the matrix elements of the screened Coulomb kernel are computed:

$$ 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},\mathbf{r'})\, \Upsilon^\text{r}_{\beta,\mathbf{q}}(\mathbf{r},\mathbf{r'}) \;, $$

where the exciton wavefunctions $\Upsilon^{\text{r}*}_{\alpha, \mathbf{q}}(\mathbf{r}, \mathbf{r'})$ are defined as above. By defining the indices as $\alpha = (v, \mathbf{k} - \mathbf{q}/2, c, \mathbf{k} + \mathbf{q}/2)$ and $\beta = (v', \mathbf{k}' - \mathbf{q}/2, c', \mathbf{k}' + \mathbf{q}/2)$, and regrouping, the product of valence and conduction wavefunctions at $\mathbf{r}$ and $\mathbf{r'}$ becomes separated by a reciprocal vector $\mathbf{p} = \mathbf{k} - \mathbf{k'}$.

In practical calculations, the screened Coulomb interaction (which includes the dielectric matrix) is expanded in plane waves:

$$ W(\mathbf{r}, \mathbf{r}') = \sum_{\mathbf{G}, \mathbf{G}'} \sum_{\mathbf{p}} e^{i(\mathbf{G} + \mathbf{p}) \mathbf{r}}\, w_{\mathbf{G}, \mathbf{G}'}(\mathbf{p})\, e^{-i(\mathbf{G}' + \mathbf{p}) \mathbf{r}'} \;. $$

Since symmetry is not exploited in constructing the BSH matrix elements, the $\mathbf{k}$ grid should be shifted by an arbitrary vector (see the vkloff attribute in the xs element) to ensure all points are symmetry-inequivalent. This improves convergence with respect to the $\mathbf{k}$-grid.

How to choose the number of Processes

The distribution for the ground-state one-shots and momentum matrix elements is limited by the number of irreducible $\mathbf{k}$ points $N_\mathbf{k}$. The number of processes $N_p$ should be chosen such that:

$$ N_\mathbf{k} \bmod N_p = 0 \;. $$

As for the computation of the screened Coulomb interaction, this depends only on differences between $\mathbf{k}$ points. Therefore, as the offset is subtracted, symmetries can be exploited, and the computation can be splitted into $N_\mathbf{k}^{\text{red}}$ independent parts. The number of processes should be chosen such that:

$$ N_\mathbf{k}^{\text{red}} \bmod N_p = 0 \;. $$

The BSH consists of $N_\mathbf{k}^2$ blocks. For MPI distribution, the matrix is split on a square grid, assigning one process to each block. The upper limit of parallelization is thus the number of $\mathbf{k}, \mathbf{k}'$ blocks. A good practice is to choose the number of processes $N_p$ as a perfect square, such that:

$$ N_\mathbf{k}^2 \bmod N_p = 0 \;. $$

Finding $N_\mathbf{k}$ and $N_\mathbf{k}^{red}$

To get the number of non-reduced $\mathbf{k}$ points, simply take the product of the three numbers defining ngridkq attribute in the xs element:

$$ N_\mathbf{k} = \text{ngridk(1)}\cdot\text{ngridk(2)}\cdot\text{ngridk(3)} $$

To get the number of reduced $\mathbf{p}$ points, we need to test first how the points reduce with symmetry operations. Doing this by hand would be a bit tedious. Therefore we implement the dryrun attribute in the screening element. If set to "true" $\texttt{exciting}$ terminates the calculation after setting up the grids in the "screen" task and prints out the number of reduced \mathbf{p}$ points as an error message:

Error(screenlauncher): input%xs%screening%dryrun is set to "true". Set to "false" to continue calculation. Number of q points: 3

Alternatively, you find the number of reduced $\mathbf{p}$ points in the end of the INFOXS.OUT file:

================================================================================
| exciting NEON started for task screen ( 430)                                 =
| version hash id: badc2e8d4742d1e7cb7c0db0b00cfe3fff5edac                     =
| MPI version using      1 processor(s)                                        =
| |  using MPI-2 features                                                      =
|                                                                              =
| Date (DD-MM-YYYY) : 03-07-2025                                               =
================================================================================

Info(xsinit): mapping screening-specific parameters

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Info(screenlauncher): Calculating screening for unshifted q-grid.
Info(screenlauncher): Number of q points:       3
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Since $\texttt{exciting}$ terminates after printing this information, there is no end tag of the "screen" task section.

Note that the $\mathbf{p}$ points are called q points in $\texttt{exciting}$ internally.

How to execute the three Parts independently

To run the three parts independently, you can use the plan element in the xs section and select only the required tasks. For ground-state one-shots and momentum matrix elements, use:

   ...
   <plan>
      <doonly task="scrgeneigvec"/>
      <doonly task="scrwritepmat"/>
      <doonly task="xsgeneigvec"/>
      <doonly task="writepmatxs"/>
   </plan>
   ...

To calculate the dielectric matrix:

   ...
   <plan>
      <doonly task="screen"/>
   </plan>
   ...

To solve the BSE:

   ...
   <plan>
      <doonly task="scrcoulint"/>
      <doonly task="exccoulint"/>
      <doonly task="bse"/>
   </plan>      
   ...

Hints

  • For computing the dielectric matrix on the reduced $\mathbf{k}$ grid, ensure that the reduceq attribute in the xs element is set to true.
  • Set skipdoneq in the screening element to true to skip already computed $\mathbf{q}$ points.
  • To check the number of reduced $\mathbf{k}$ points, start the "screen" task, stop it immediately, and inspect the INFOXS.OUT file.
  • Disable distribution of the BSH by setting the distribute attribute to "false" in the BSE element.
  • To verify that the solution of the BSE was indeed distributed, check the "bse" task in the INFOXS.OUT file. If the BSH was distributed you will find the following line:
Info(bse): Invoking scalapack routine PZHEEVX

else

Info(bse): Invoking Lapack routine ZHEEVR
  • If your system allows for a hybrid OMP and MPI parallelization, setting NUM_OMP_THREADS=4 is usually a good choice. For optimal settings you need to do a scaling test on the mashine you plan to use.

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 B. Aurich's Master thesis PDF.
  • 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.