Excited States from BSE

for lithium version of Exciting

---

Purpose: In this tutorial you will learn how to perform a basic Bethe-Salpeter-equation (BSE) calculation. As an example, the optical spectrum of LiF will be studied. Moreover, the BSE-derived xc kernel will be used for TDDFT.

---

---

0. Preparation

Proceed in the same way as for Excited States from TDDFT.

---

1. Groundstate calculation

First, a groundstate calculation has to be performed. Make yourself familiar again with the corresponding part of Excited States from TDDFT. For convenience, the part of the input file concerning the element structure is provided in the following. Since in LiF the lattice has a basis composed of two atoms, one additional species has to be added (see below).

  <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" bfcmt="0.0 0.0 0.0" />
    </species>
    <species speciesfile="F.xml">
      <atom coord="0.5000  0.5000  0.5000" bfcmt="0.0 0.0 0.0" />
    </species>
  </structure>

N. B.: Do not forget to replace in the input.xml the string "$EXCITINGROOT" by the actual value of the environment variable $EXCITINGROOT. In order to do this, use the command

$ SETUP-excitingroot.sh

We further recommend a finer radial grid, higher G-vector cutoff for the density and the potential, as well as a higher lmaxvr value for density and potential (check Input Reference to get familiar with the parameters used here).

...
<groundstate ... lradstep="2" gmaxvr="14.0" lmaxvr="8" ...>
...

For the k-point mesh you can choose the same as for the TDDFT example for the ground state.

...
<groundstate ... ngridk="..." ...>
...

Complete these pieces of input generating a valid input file out of it. Start the ground state SCF calculation in a directory of your choice as described in Excited States from TDDFT. Check if the calculation finishes gracefully and if the STATE.OUT and EFERMI.OUT files are present. These two files are (as in the case of TDDFT) the starting point of the BSE calculation.

---

2. Bethe-Salpeter equation

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

Work-flow

Similarly to the TDDFT, the work-flow is as follows: The eigenvalues and wavefunctions are determined for a given k-mesh. The momentum matrix elements are then calculated. In a next step, the q-dependent (static) screening is calculated, where the q-dependent matrix elements are calculated on the fly. After that, the BSE Hamiltonian is set up, consisting of two parts which have to be calculated explicitly: The (attractive) direct term which depends on the screening, and the (repulsive) exchange term. In a last step, the eigenvalue problem for the BSE Hamiltonian is solved and the spectral summation is performed for the macroscopic dielectric function. See a summary in the following figure.

Constructing the input file and running

The actual setup of the calculation could look as follows (solely the xs element is given here, you have to add it to the input file at the proper place inside the input element).

...
<xs xstype="BSE" 
    ngridq="4 4 4"
    ngridk="4 4 4" vkloff="0.05 0.15 0.25" 
    nempty="3"
    lmaxapwwf="3"
    gqmax="3.0"
    broad="0.0073499"
    scissor="0.20947"
    tevout="true" nosym="true">
 
    <energywindow intv="0.0 1.0" 
                  points="1200" />
 
    <screening screentype="full"
               nempty="115" />
 
    <BSE bsetype="singlet"
         nstlbsemat="1 5 1 4"
         nstlbse="1 5 1 4" 
         aresbse="false"/>
 
    <qpointset>
      <qpoint>0.0 0.0 0.0</qpoint>
    </qpointset>
 
</xs>                        
...

Now, execute the program binary again for the modified input (do not forget to skip the groundstate calculation to save computation time).

$ excitingser >& outputXS &

While the program is running make yourself familiar with the input parameters appearing in above using Input Reference and check if and how the chosen values differ from the default values.

A few additional remarks

[show more]

• The k-point mesh is crucial for a good resolution of the spectra, though choosing a dense mesh can result into a time-consuming calculation. This is especially true for the BSE (see the scaling information below).
• The q-point mesh is chosen to be the same as the k-point mesh. It must be a commensurate super-mesh with respect to the k-point mesh in the current implementation. Keep it in sync with the k-point mesh!
• The symmetry relations between the k-points are not exploited in the code, but for the case of the screening which is calculated only for a reduced set of q-points. In any case, it is a common practice - also for the BSE - to shift the k-mesh off symmetry as in the case of TDDFT.
• The lmaxapwwf parameter is reduced since the full lmaxapw value (the default) is rarely needed.
• A scissors correction of 5.7 eV = 0.20947 Ha is defined to mimic the quasi-particle energy shifts (assumed to be constant).
• The attribute bsetype is set to "singlet". It can be changed to "RPA" at a later stage to compare with the TDDFT calculation.
• We choose 5 valence and 4 conduction bands. Check the definition of the nstlbse and nstlbsemat parameters in Input Reference. The nstlbse parameter vector (a list of 4 integers) must be less or equal to nstlbsemat for each of the entries. This is because with the nstlbse parameter one can select a subset of valence/conduction states for the diagonalization, after the parts of the Hamiltonian have already been calculated.
• Once the screening is calculated, use the do="skip" attribute for the screening element to avoid calculating it again.

Analyzing results

Once the calculation has finished, a xmgrace plot can be generated as follows

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

In this way, three xmgrace files, namely, dielectricfunction_Re.agr, dielectricfunction_Im.agr, and dielectricfunction_ReKK.agr, are generated. These files contain the real part, the imaginary part, and the real part from a Kramers-Kronig transform of the dielectric function, respectively.

They are viewed by invoking the xmgrace command (see Xmgrace: A Quickstart for help on xmgrace). We are mostly interested in the imaginary part.

$ xmgrace dielectricfunction_Im.agr

Here is the output plot:

For further analysis it can be useful to look directly at excitation energies, i.e., at eigenenergies of the BSE Hamiltonian. They are printed out in the files EXCITON_NAR_BSEsinglet_SCRfull_OC11.OUT and EXCITON_SORTED_NAR_BSEsinglet_SCRfull_OC11.OUT. These two files contain the eigenvalues and the oscillator strengths of the individual excitation. In the third column, the energy of the bandgap is subtracted and negative values indicate bound excitons. The last column contains the magnitue of the oscillator strength. The first file is sorted with respect to the eigenenergies, whereas the _SORTED file is sorted with respect to the oscillator strengths and allows for a quick identification of the eigenvalues with the strongest oscillator strengths. Check if you can identify a bound exciton in the case of LiF.

Converging results

The spectra obtained from the BSE 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 in the setup of the Hamlitonian (direct term), and, if a full diagonalization is applied to the eigenvalue problem, even with the third power.

The most important convergence parameters are listed here as a supplement to those given already for the TDDFT case.

Description Parameter Info k-mesh/q-mesh input/xs/@ngridk, input/xs/@ngridq very crucial but cannot always be converged energy cutoff (screening) input/xs/screening/@nempty has to be converged local field effects G-cutoff input/xs/@gqmax defines the quality of the screened Coulomb potential in G-space states below/above the Fermi level input/xs/BSE/@nstlbse, input/xs/BSE/@nstlbsemat affects the range of the spectrum towards higher energies and possible mixing of transitions

A 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

(1)

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

Excercise

Once you successfully run the previous BSE calculation you can play around with the input parameters to achieve a better quality for the dielectric function.

• Change the bsetype to calculate the RPA spectrum with and without (IP) local field effects (two different calculations, since the BSE Hamiltonian is different).
• Try to choose different values for the nstlbse parameter to narrow the range of states and check how the spectra change.
• Decrease the local fields cutoff gqmax, how do the results change.
• If you have time you can try to increase the k-mesh/q-mesh (use the same values for both), though, doing a 6$\times$6$\times$6 mesh might already take a couple of hours on one processor.

---

3. TDDFT using the BSE-derived xc kernel

We will calculate the dielectric function of LiF now via TDDFT but utilizing the BSE-derived kernel, which requires many parts of a standard BSE calculation.

The work-flow of the algorithm is a combination of the TDDFT linear response calculation (see Excited States from TDDFT) and the calculation of the direct term of the BSE Hamiltonian, which is then used to set up a MBPT-derived kernel in first order. This kernel then enters the Dyson equation for the response function in the last stage of the TDDFT formalism.

Since we do not need to calculate the groundstate again, we add the do attribute setting it to "skip".

...
<groundstate do="skip" ... >
...

For the TDDFT calculation we change the xstype attribute from "BSE" to "TDDFT".

...
<xs xstype="TDDFT"  ... >
...

and add the tddft element inside the xs element:

...
<tddft fxctype="MB1" aresdf="false" aresfxc="false" />
...

Again, make yourself familiar with input parameters you do not know so far. We skip the anti-resonant parts of the kernel and also for the BSE in the previous calculation for a simpler comparison (ares* attributes above). A hint for impatient users: The calculation of the screening can be skipped since it has already been performed for the BSE. Simply add the do attribute and set it to "skip" inside the screening element.

...
<screening do="skip" ... >
...

Having modified the input file accordingly we are now in the position to start the calculation for the TDDFT spectrum using the BSE-derived xc kernel. Simply execute

$ excitingser >& outputXS &

We again are in the position to generate a xmgrace graphics file.

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

The resulting plot is

Partly converged results with respect to the size of the k-point basis set are given for comparison in Converged TDDFT and BSE Results for LiF.

Excercise

• Additionally, you can calculate the excitonic spectrum of bulk Silicon, a standard example for the BSE. You can start from the LiF input and change the structure element. The basis vectors within crystal can remain the same. Change both atomic species to Si and check out the lattice constant of Si from WebElements. The lattice constant (in Bohr) can be given directly as scaling parameter in the attribute scale of the element crystal. Then you can perform the following actions.
  1. Do the groundstate first.
  2. Re-think the choice of parameters for the excited state compared to LiF.
  3. Start with a 2$\times$2$\times$2 k-mesh and increase it to 4$\times$4$\times$4 for the BSE.
  4. Try a larger k-mesh 4$\times$4$\times$4 for a long calculation.
  5. In the BSE approach, use the RPA and singlet type of Hamiltonian for the spectra.
  6. Compare your RPA result obtained within the BSE framework to the RPA result obtained using TDDFT.

---

Literature

• Tutorial talk by Stephan Sagmeister (PDF)
• Details on implementation and application of the BSE formalism to GaAs: 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)

---