14.2 Excitation energies from the Bethe–Salpeter equation

14.2.1 Theoretical background

The Bethe–Salpeter (BS) excitation energies $\omega_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 general references are given below. The orbitals are assumed to be complex. \[\begin{equation} \left(\begin{array}{cc} {\bf A} & {\bf B} \\ {\bf B*} & {\bf A*} \end{array}\right) \left(\begin{array}{c}{\bf X}_n\\{\bf Y}_n\end{array}\right) = \omega_n \left(\begin{array}{cc} {\bf 1} & \phantom{-}{\bf 0} \\ {\bf 0} & -{\bf 1} \end{array}\right) \left(\begin{array}{c}{\bf X}_n\\{\bf Y}_n\end{array}\right) \; . \end{equation}\](14.7) In the following, the indices $a,b,\ldots$ ($i,j,\ldots$) refer to unoccupied (occupied) spinors, the indices $p,q,r,\ldots$ include occupied and unoccupied spinors. The definitions of the matrices ${\bf A}$ and ${\bf B}$ differ in the one-particle energies, which are the orbital energies in case of TD-HF and the $GW$ quasi-particle energies ${\bf \varepsilon}^{\text{QP}}$ in case of BS. Also, the exchange interaction is screened in case of BS, \[\begin{align} A_{ia,jb} &= \Delta \varepsilon_{ia,jb}^{\text{QP}} + (ai|jb) - (ji\bar{|}ab) \; ,\\ B_{ia,jb} &= (ai|bj) - (ja\bar{|}ib) \; . \end{align}\](14.8–14.9) The exchange interaction $(pq\bar{|}rs)$, which is defined as \[\begin{equation} (pq\bar{|}rs) = \sum_{tu} ({\bf \epsilon}^{-1})_{pq,tu} (tu|rs) \Leftrightarrow {\bf W} = {\bf \epsilon}^{-1} {\bf v} \; , \end{equation}\](14.10) is screened by the inverse of the dielectric function, \[\begin{align} \epsilon_{pq,rs} = \delta_{pr}\delta_{qs} - (pq|rs) \left[\chi_0(\omega=0)\right]_{rs,rs} \; , \\ \left[\chi_0(\omega=0)\right]_{rs,rs} = \sum_{ia} \frac{\delta_{ra}\delta_{si} + \delta_{ri}\delta_{sa}}{\varepsilon_i^{\text{QP}} - \varepsilon_a^{\text{QP}}} \; , \end{align}\](14.11–14.12) which depends on the $GW$ quasi-particle energies. It can be expressed in terms of an infinite-order perturbation expansion in the Coulomb interaction, \[\begin{equation} \begin{split} {\bf W} =& {\bf \epsilon}^{-1} {\bf v} \\ =& {\bf v} + {\bf v} {\bf \chi}_0(\omega=0) {\bf v} + {\bf v} {\bf \chi}_0(\omega=0) {\bf v} {\bf \chi}_0(\omega=0) {\bf v} + \ldots \; . \end{split} \end{equation}\](14.13) To obtain ${\bf W}$, the RI-approximation is introduced for the two-electron integrals. As a direct result, the inversion ${\bf \epsilon}$ is replaced by the inversion of a matrix in the auxiliary basis: \[\begin{equation} \begin{split} {\bf W} \approx {\bf b}^\top {\bf \tilde{i}} \; {\bf b} =& {\bf b}^\top \big({\bf 1}+{\bf i}+{\bf i}^2+{\bf i}^3+\ldots\big) {\bf b} \\ =& {\bf b}^\top \big({\bf 1}-{\bf i}\big)^{-1}{\bf b} \; , \end{split} \end{equation}\](14.14) \[\begin{equation} {\bf i} = 2 {\bf b} {\bf \chi}_0(\omega=0) {\bf b}^\top \Leftrightarrow i_{PQ} = 2 \sum_{ia} \frac{b_{P,ai}b_{Q,ai}}{\varepsilon_i^{\text{QP}} - \varepsilon_a^{\text{QP}}} \; . \end{equation}\](14.15) For further details, we refer to Ref. [34].

14.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 D$_{\text{2h}}$ 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 D$_{\text{2h}}$. Two-component references ($soghf) including spin-orbit interactions can be used.

Excitations

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\times 2$ “super-matrix” can be computed, or of ${\bf 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. For two-component calculations only the TDA variant is implemented. Note that 2c-BSE is, compared to TD-DFT, not limited to closed shell cases but also implemented for open shell cases.

Quasiparticle energies

The BSE method reads in quasiparticle energies obtained from a preceding $GW$ calculation. It supports all different $GW$ methods available in escf. Since Turbomole 7.3 the number of frozen orbitals can be different in the $GW$ and $BSE$ step. However qpenergies.dat files from Turbomole 7.2 and earlier are not directly compatible with Turbomole 7.3 and higher, and it is recommended to redo the $GW$ calculation. The reason for this is that the quasiparticle states for frozen orbitals are also projected since Turbomole 7.3, and always all orbitals are therefore updated and reported in the qpenergies.dat file.

Quasiparticle energies in ${\bf A}$ and ${\bf W}$

By default, the matrices ${\bf A}$ and ${\bf W}$ are set up using quasiparticle energies. To set up ${\bf A}$ using KS/HF orbital energies instead, i.e. \[\begin{equation} A_{ia,jb} = \Delta \varepsilon_{ia,jb} + (ai|jb) - (ji\bar{|}ab) \; ,\\ \end{equation}\](14.16) the option noqpa has to be added to the $bse block in control. Similarly, if the denominator in ${\bf \chi_0}$ and, consequently, the static screened interaction should be constructed using KS/HF orbital energies instead of $GW$ quasiparticle energies, the option noqpw has to be set. The quasiparticle 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 ${\bf \tilde{i}}$. By default, this is computed by Cholesky decomposition and inversion of $({\bf 1} - {\bf i})^{-1}$. An alternative iterative procedure, \[\begin{equation} {\bf \tilde{i}}^{(i+1)} = {\bf 1} + {\bf i} \; {\bf \bar{i}}^{(i)} \; , \end{equation}\](14.17) 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.

14.2.3 General recipe for a BSE calculation

dscf or ridft calculation to converge reference orbitals
Provide $gw or $rigw control flags and perform a $GW$ calculation
Choose reasonable cbas auxiliary basis sets if not done in the $GW$ step
escf calculation to obtain qpenergies.dat
Remove $gw and $rigw control flags
Set BSE control flags
Add the keyword $rik or $rick if not already present
escf calculation to obtain excitation energies