# BSE for Molecules

In this tutorial we will show the reader how to perform the Bethe-Salpeter equation (BSE) calculation on molecules using FHI-aims to study the neutral excitations. The implementation of the BSE in FHI-aims based on NAO are described by Liu, et al. 1

## Background

The physical process described by the BSE equation is the charge neutral optical excitation as shown in this figure.

Figure from Golze et al. 2

As shown, the optical excitation involves the interaction between a hole and a excited electron. To effectively describe such phenomenon, the electron-hole vector space, which is the combination of a occupied orbital and an unoccupied orbital, is a natural choice to expand the BSE equation. In our implementation, the BSE equation is written as

$\left[\begin{array}{cc} A & B \\ -B^{\dagger} & -A^{\dagger} \end{array}\right]\left[\begin{array}{c} X_{s} \\ Y_{s} \end{array}\right]=E_{s}\left[\begin{array}{c} X_{s} \\ Y_{s} \end{array}\right]$

where, $$X_s$$ and $$Y_s$$ are the eigenvectors within the electron-hole vector space. With the Tamm-Dancoff approximation3, we ignore the coupling matrix $$B$$, thus a reduced matrix can be solved. Let us look at the elements in $$A$$ matrix.

$A_{i a}^{j b}=\left(\varepsilon_{a}^{G W}-\varepsilon_{i}^{G W}\right) \delta_{i j} \delta_{a b}-\alpha^{S / T}(i a|V| j b)+(i j|W(\omega=0)| a b),$

where $$|i a)$$ and $$|j b)$$ are the electron-hole vector space prepared from the GW quasi-particle wavefunctions. $$\varepsilon_{i}$$ is the GW quasiparticle energies. $$V$$ and $$W$$ are the Coulomb and screened Coulomb interactions.

The eigenvalues of the BSE matrix are the calculated optical excitation energies.

## Exercise

In this BSE exercise, you will perform a BSE calculation based on $$G_0W_0$$ calculation which is started from PBE exchange-correlation functional DFT calculation (BSE@$$G_0W_0$$@PBE). Tier 2 + aug2 basis set for ethylene C$$_2$$H$$_4$$ is used. The quantities of interest in this case are the BSE eigenvalues. You can proceed as follows:

• Generate the geometry.in file for ethylene from the experimental data (extracted from here):
atom    0.0000  0.0000  0.6695 C
atom    0.0000  0.0000  -0.6695 C
atom    0.0000  0.9289  1.2321 H
atom    0.0000  -0.9289 1.2321 H
atom    0.0000  0.9289  -1.2321 H
atom    0.0000  -0.9289 -1.2321 H

• Here are the commands for the BSE calculations in the control.in file. As can be seen, in general, we separate the commands into three parts: The DFT part, the GW part, and the BSE part. There is also the basis set information for different species. As recommended in our paper,1 we use Tier2 + aug2 as a relatively converged basis sets for molecular BSE calculation. We can find these basis set information in the folder AIMS_DIR/species_defaults/non-standard/Tier2_aug2.
  #########################################
xc                 pbe
spin none
relativistic   atomic_zora scalar
calculate_all_eigenstates # Make sure all the unoccupied orbitals are calculated
#########################################

#########################################
qpe_calc           gw     # G0W0
n_anacon_par 16
#########################################
neutral_excitation    bse
bse_s_t singlet                # singlet or triplet excitation
#########################################

• Run the calculation with the FHI-aims executable.
mpirun -n 4 aims.x | tee output

• Read the first excitation energy from the output. Given the best theoretical estimation of the first excitation energy to be 7.8 eV.4, how good is the result from BSE? Here is an example output, where the 6.1698 eV is the first excitation state.
  |        1. Singlet:      6.1698 eV     [  0.2267 Ha]  -  Oscillator Strength:       0.0527
|            Transition Moments    X:   -0.2811   Y:    0.0000   Z:    0.0000
|   Dominant single-particle level transitions in this excitation
|   eigenvector: [occupied ==> virtual, eigenvector element]
|     8==>    10   -0.6261
|     8==>    15    0.2002
|     8==>    19    0.2304

• Optional: Start the calculation with the DFT calculation with PBE0 functional. i.e. BSE@$$G_0W_0$$@PBE0. Play around with the hybrid mixing parameter in PBE0. We expect the hybrid mixing parameter to have an influence on the BSE first excitation energy. Here is a reference on this topic 5. To perform a BSE@$$G_0W_0$$@PBE0 calculation, we need to modify the xc pbe line in control.in to these two lines.
    xc                 pbe0
hybrid_xc_coeff    0.25

• For convenience, here we provide the full control.in file where xc part can be modified accordingly to perform above exercises.
  #########################################
xc                 pbe
spin none
relativistic   atomic_zora scalar
calculate_all_eigenstates
#########################################

#########################################
qpe_calc           gw
anacon_type  1
n_anacon_par 16
#########################################
neutral_excitation    bse
bse_s_t singlet
#########################################
################################################################################
#
#  FHI-aims code project
#  VB, Fritz-Haber Institut, 2007
#
#  Suggested "safe" defaults for H atom (to be pasted into control.in file)
#
################################################################################
species        H
#     global species definitions
nucleus             1
mass                1.00794
#
l_hartree           8
#
cut_pot             4.0  2.0  1.0
basis_dep_cutoff    0.d0
#
angular_grids       specified
division   0.1930   50
division   0.3175  110
division   0.4293  194
division   0.5066  302
division   0.5626  434
division   0.5922  590
division   0.6227  974
division   0.6868 1202
outer_grid  1202

################################################################################
#
#  Definition of "minimal" basis
#
################################################################################
#     valence basis states
valence      1  s   1.
#     ion occupancy
ion_occ      1  s   0.5
################################################################################
#
#  Suggested additional basis functions. For production calculations,
#  uncomment them one after another (the most important basis functions are
#  listed first).
#
#  Basis constructed for dimers: 0.5 A, 0.7 A, 1.0 A, 1.5 A, 2.5 A
#
################################################################################
#  "First tier" - improvements: -1014.90 meV to -62.69 meV
hydro 2 s 2.1
hydro 2 p 3.5
#  "Second tier" - improvements: -12.89 meV to -1.83 meV
hydro 1 s 0.85
hydro 2 p 3.7
hydro 2 s 1.2
hydro 3 d 7
#  "Third tier" - improvements: -0.25 meV to -0.12 meV
#     hydro 4 f 11.2
#     hydro 3 p 4.8
#     hydro 4 d 9
#     hydro 3 s 3.2

###############################################################################
#
# Gaussian augmentation functions.
# Obtained by diff cc-pV5Z aug-cc-pV5Z
#
###############################################################################

gaussian 0 1 0.0207000
gaussian 1 1 0.0744000

################################################################################
#
#  FHI-aims code project
#  VB, Fritz-Haber Institut, 2007
#
#  Suggested "safe" defaults for C atom (to be pasted into control.in file)
#
################################################################################
species        C
#     global species definitions
nucleus             6
mass                12.0107
#
l_hartree           8
#
cut_pot             4.0  2.0  1.0
basis_dep_cutoff    0.d0
#
angular_grids specified
division   0.2187   50
division   0.4416  110
division   0.6335  194
division   0.7727  302
division   0.8772  434
division   0.9334  590
division   0.9924  770
division   1.0230  974
division   1.5020 1202
outer_grid  1202

################################################################################
#
#  Definition of "minimal" basis
#
################################################################################
#     valence basis states
valence      2  s   2.
valence      2  p   2.
#     ion occupancy
ion_occ      2  s   1.
ion_occ      2  p   1.
################################################################################
#
#  Suggested additional basis functions. For production calculations,
#  uncomment them one after another (the most important basis functions are
#  listed first).
#
#  Constructed for dimers: 1.0 A, 1.25 A, 1.5 A, 2.0 A, 3.0 A
#
################################################################################
#  "First tier" - improvements: -1214.57 meV to -155.61 meV
hydro 2 p 1.7
hydro 3 d 6
hydro 2 s 4.9
#  "Second tier" - improvements: -67.75 meV to -5.23 meV
hydro 4 f 9.8
hydro 3 p 5.2
hydro 3 s 4.3
hydro 5 g 14.4
hydro 3 d 6.2
#  "Third tier" - improvements: -2.43 meV to -0.60 meV
#     hydro 2 p 5.6
#     hydro 2 s 1.4
#     hydro 3 d 4.9
#     hydro 4 f 11.2
#  "Fourth tier" - improvements: -0.39 meV to -0.18 meV
#     hydro 2 p 2.1
#     hydro 5 g 16.4
#     hydro 4 d 13.2
#     hydro 3 s 13.6
#     hydro 4 f 17.6
#  Further basis functions - improvements: -0.08 meV and below
#     hydro 3 s 2
#     hydro 3 p 6
#     hydro 4 d 20

###############################################################################
#
# Gaussian augmentation functions.
# Obtained by diff cc-pV5Z aug-cc-pV5Z
#
###############################################################################

gaussian 0 1 0.0394000
gaussian 1 1 0.0272000