Excited States from BSE

for helium 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. Preparations…

Check them out at Excited States from TDDFT.

1. Ground state

First, a ground-state calculation has to be performed. Make yourself familiar again with the corresponding part of the TDDFT tutorial. For convenience the structure input is provided. 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>

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

<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 the TDDFT tutorial. 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. BSE

We will calculate the macroscopic dielectric function of LiF now via the BSE formalism. 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 ${\bf q}$-dependent (static) screening is calculated, where the ${\bf 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.

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 - within 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 to save computation time).

$ excitingser >& outputXS &

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

A few additional remarks

• 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.7eV = 0.20947H is defined to mimic the quasi-particle energy shifts (assumed to be constant).
• The 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 the 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.

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

$ xsltproc --stringparam filename dielectricfunction $EXCITINGROOT/xml/visualizationtemplates/xmldielectricfunction2agr.xsl EPSILON_NAR_BSEsinglet_SCRfull_OC11.OUT.xml

This template will create three xmgrace files, namely

dielectricfunction_Re.agr
dielectricfunction_Im.agr
dielectricfunction_ReKK.agr

These files contain the real part, the imaginary part and the real part from a Kramers-Kronig transform of the dielectric function.

They are viewed by invoking the xmgrace command (see xmgrace-quickstart for help on xmgrace). We are mostly interested in the imaginary part.

$ xmgrace dielectricfunction_Im.agr

For further analysis it can be useful to look directly at the excitation energies, i.e. the eigenenergies from the BSE Hamiltonian. They are printed out in the exciton files

EXCITON_NAR_BSEsinglet_SCRfull_OC11.OUT
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 the 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 k-point set 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	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 6x6x6 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 the tutorial 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 the 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 xmgrace graphics files:

$ xsltproc --stringparam filename dielectricfunction-tddft $EXCITINGROOT/xml/visualizationtemplates/xmldielectricfunction2agr.xsl EPSILON_NAR_FXCMB1_OC11_QMT001.OUT.xml

At the end of this tutorial we give you results which are partly converged with respect to the k-point set for comparison: Converged 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 <strucutre> definition. 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 used directly as scaling parameter within crystal.

• do the groundstate first
• re-think the choice of parameters for the excited state compared to LiF
• start with a 2x2x2 k-mesh and increase to 4x4x4 for the BSE
• try a larger k-mesh 6x6x6 for a long calculation
• use the RPA and singlet type of BSE Hamiltonian for the spectra
• compare your RPA result obtained from the BSE framework to the RPA result obtained from the TDDFT framework

Literature

Tutorial Talk 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)

4. Outlook

• Soon the whole TDDFT/BSE part will be extended to the spin-polarized case.