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

Background

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

optical_absorption_figure

Figure from Golze et al. ²

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 approximation³, 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,¹ 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
    anacon_type  1            # PADE
    n_anacon_par 16
  #########################################
    neutral_excitation    bse
    read_write_qpe rw
    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.⁴, 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 ⁵. 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
    read_write_qpe rw
    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
#
    radial_base         24 7.0
    radial_multiplier   6
    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
#
    radial_base         34 7.0
    radial_multiplier   6
    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