13.2 Excitation energies from the Bethe–Salpeter equation

13.2.1 Theoretical background

The Bethe–Salpeter (BS) excitation energies ωn are obtained as solutions of a pseudo-Hermitian eigenvalue equation, which has the same form as in time-dependent Hartree–Fock (TD-HF) theory. For the sake of simplicity, only the equations for singlet excitation energies based on a closed-shell reference are given below. The orbitals are assumed to be real.

(       ) (     )      (       ) (     )
  A   B      Xn          1   0      Xn
  B   A      Yn   = ωn   0  - 1     Yn   .

In the following, the indices a,b, (i,j,) refer to unoccupied (doubly occupied) spatial orbitals, the indices p,q,r, include both doubly and unoccupied orbitals. The definitions of the matrices A and B differ in the one-particle energies, which are the orbital energies in case of TD-HF and the GW quasi-particle energies εQP in case of BS. Also, the exchange interaction is screened in case of BS,

Aia,jb  =  Δ εia,jb + 2(ai|jb)- (ji|ab),                (13.8)
Bia,jb  =  2(ai|bj)- (ja|ib).                         (13.9)
The exchange interaction (pq|rs), which is defined as
(pq|rs) = ∑  (ϵ-1)   (tu|rs) ⇔ W  = ϵ- 1v,
        tu     pq,tu

is screened by the inverse of the dielectric function,

ϵ    = δ  δ - 2(pq|rs)[χ (ω = 0)]  ,                 (13.11)
 pq,rs   prqs        ∑  0       rs,rs
    [χ0(ω = 0)]   =     δraδsi +-δriδsa,               (13.12)
             rs,rs   ia  εQiP - εQaP
which depends on the GW quasi-particle energies. It can be expressed in terms of an infinite-order perturbation expansion in the Coulomb interaction,
W  =ϵ-1v

   =v + 2vχ0(ω = 0)v+ 4vχ0(ω = 0)vχ0 (ω = 0)v+ ....

To obtain W, the RI-approximation is introduced for the two-electron integrals. As a direct result, the inversion ϵ is replaced by the inversion of a matrix in the auxiliary basis:

     ⊤ ˜   ⊤(       2   3     )
W  ≈ b ib =b  1 +i +i + i + ... b
         =b⊤(1 - i)-1b,

                         ∑   b   b
i = 2bχ0(ω = 0)b⊤ ⇔ iPQ = 4  -QPP,aiQ,aQiP.
                          ia εi  - εa

For further details, we refer to Ref.  [178].

13.2.2 BSE features

The BSE method is implemented in the escf module. It is activated by the keyword $bse. Its features and different variants are discussed in the following.

Reference determinant

The BSE method supports both Hartree–Fock and Kohn–Sham references, that is, orbitals obtained from a dscf or ridft computation. Note that symmetry is only supported up to D2h and its subgroups. In case of wave functions of higher symmetry, the orbitals have to be converted using define’s option susy to a subgroup of D2h.


Singlet, triplet and spin-unrestricted excitation energies based on the Bethe–Salpeter equation are available in escf. Similar to the RPA and TDA variants in TD-HF, the eigenvalues of the full 2 × 2 “super-matrix” can be computed, or of A only. To define these options, the same input blocks as for a TD-HF computation are read in: $scfinstab and $soes. See Section 8 for details.

Quasi-particle energies

The BSE method reads in quasi-particle energies obtained from a preceding GW calculation. It supports all different GW methods available in escf. Note that if orbitals are excluded from excitations by the flag $freeze, the same number of occupied orbitals have to be frozen in the GW and BSE calculations.

Quasi-particle energies in A and W

By default, the matrices A and W are set up using quasi-particle energies. To set up A using KS/HF orbital energies instead, i.e.

Aia,jb = Δεia,jb + 2(ai|jb)- (ji|ab),

the option noqpa has to be added to the $bse block in control. Similarly, if the denominator in χ0 and, consequently, the static screened interaction should be constructed using KS/HF orbital energies instead of GW quasi-particle energies, the option noqpw has to be set. The quasi-particle energies are read in from the file qpenergies.dat, which has to be created in a preceding GW calculation using escf. A different file name can be defined using the option file in the $bse block.

Computation of the static screened interaction

The static screened interaction is obtained from the auxiliary matrix ĩ. By default, this is computed by Cholesky decomposition and inversion of (1 - i)-1. An alternative iterative procedure,

˜i(i+1) = 1+ ii(i),

is activated by the option iterative. The convergence threshold (default: 1.0d-12) and maximum number of iterations (default: 100) can be adjusted by the options thrconv and iterlim, respectively.

Auxiliary basis sets

The static screened interaction is computed using the RI approximation, therefore, an auxiliary basis set is required. Since the integrals are similar to those in MP2, the cbas auxiliary basis sets are recommended. They can be conveniently defined using the cc sub-menu of define.

13.2.3 General recipe for a BSE calculation

dscf or ridft calculation to converge reference orbitals
Provide GW control flags
escf calculation to obtain qpenergies.dat
Remove GW control flags
Set BSE control flags
Choose reasonable cbas auxiliary basis sets
Add the keyword $rik.
escf calculation to obtain excitation energies