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

`element).`

**<input>**<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
and**nstlbse**parameters in the Input Reference. The**nstlbsemat**parameter vector (a list of 4 integers) must be less or equal to**nstlbse**for each of the entries. This is because with the**nstlbsemat**parameter one can select a subset of valence/conduction states for the diagonalization, after the parts of the Hamiltonian have already been calculated.**nstlbse** - Once the screening is calculated, use the
attribute for the**do="skip"**element to avoid calculating it again.**screening**

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)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
to calculate the**bsetype***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
parameter to narrow the range of states and check how the spectra change.**nstlbse** - Decrease the local fields cutoff
, how do the results change.**gqmax** - 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

`element:`

**<xs>**<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

`can remain the same. Change both atomic species to`

**crystal****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
k-mesh and increase to**2x2x2**for the BSE**4x4x4** - try a larger k-mesh
for a long calculation**6x6x6** - 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.