GW + Bethe-Salpeter equation

The Bethe-Salpeter equation (BSE) is a method for computing electronic excitation energies and optical absorption spectra. We repeat the theory and implementation of BSE from [Graml2026] in Sec. 1, in Sec. 2 the BSE input keywords and in Sec. 3 a full CP2K input file of a BSE calculation and the corresponding output. Sec. 4 describes the BSE section under PROPERTIES for large molecules. For further references on BSE, see [Blase2018, Blase2020, Bruneval2015, Sander2015].

1. Theory and implementation of BSE

1.1 Electronic excitation energies and the excitation wave function

A central goal of a BSE calculation is to compute excited state properties such as the electronic excitation energies \(\Omega^{(n)} (n=1,2,\ldots)\).

The following ingredients are necessary for such a computation:

Occupied Kohn-Sham (KS) orbitals \(\varphi_i(\mathbf{r})\) and empty KS orbitals \(\varphi_a(\mathbf{r})\) from a DFT calculation, where \(i=1,\ldots,N_\mathrm{occ}\) and \(a=N_\mathrm{occ}+1,\ldots,N_\mathrm{occ}+N_\mathrm{empty}\),
GW eigenvalues \(\varepsilon_i^{GW}\) and \(\varepsilon_a^{GW}\) of corresponding KS orbitals.

It is possible to use G₀W₀, evGW₀ or evGW eigenvalues, see details in GW and in Ref. [Golze2019], i.e. we perform BSE@G ₀W₀/evGW₀/evGW@DFT. Thus, also input parameters for a DFT and GW calculation influence the BSE calculation (see full input in Sec. 3.1). We obtain the properties of the excited states from the BSE by solving the following generalized eigenvalue problem that involves the block matrix \(ABBA\):

\[\begin{split} \left( \begin{array}{cc}A & B\\B & A\end{array} \right)\left( \begin{array}{cc}\mathbf{X}^{(n)}\\\mathbf{Y}^{(n)}\end{array} \right) = \Omega^{(n)}\left(\begin{array}{cc}1&0\\0&-1\end{array}\right)\left(\begin{array}{cc}\mathbf{X}^{(n)}\\\mathbf{Y}^{(n)}\end{array}\right) \quad . \end{split}\]

We abbreviate \(A\) and \(B\) as matrices with index \(A_{ia,jb}\), i.e. they have \(N_\mathrm{occ}N_\mathrm{empty}\) rows and \(N_\mathrm{occ}N_\mathrm{empty}\) columns. The entries of \(A\) and \(B\) are

\[\begin{split} A_{ia,jb} &= (\varepsilon_a^{GW}-\varepsilon_i^{GW})\delta_{ij}\delta_{ab} + \alpha^\mathrm{S/T} v_{ia,jb} - W_{ij,ab}(\omega=0) \quad ,\\ B_{ia,jb} &= \alpha^\mathrm{(S/T)} v_{ia,bj} - W_{ib,aj}(\omega=0) \quad . \end{split}\]

where \(\delta_{ij}\) is the Kronecker delta. The user sets \(\alpha^S=2\) for computing singlet excitations and \(\alpha^T=0\) for computing triplet excitations. \(v_{pq,rs}\) is the bare Coulomb interaction and \(W_{pq,rs}(\omega=0)\) the statically (\(\omega=0\)) screened Coulomb interaction, where \(p,q,r,s \in [ 1, N_\mathrm{occ}+N_\mathrm{empty}]\) are KS orbital indices. Note here, that the screened Coulomb interaction is always computed from DFT quantitites, i.e. \(W_{0}(\omega=0)\) enters the BSE. \((\mathbf{X}^{(n)},\mathbf{Y}^{(n)})\) with elements \(X_{ia}^{(n)}\) and \(Y_{ia}^{(n)}\) are the eigenvectors of the excitation \(n\), which relate to the wave function of the electronic excitation,

\[ \Psi_\text{excitation}^{(n)}(\mathbf{r}_e,\mathbf{r}_h) = \sum_{ia} X_{ia}^{(n)} \varphi_i(\mathbf{r}_h) \varphi_a(\mathbf{r}_e) + Y_{ia}^{(n)} \varphi_i(\mathbf{r}_e) \varphi_a(\mathbf{r}_h) \quad , \]

i.e. \(X_{ia}^{(n)}\) and \(Y_{ia}^{(n)}\) describe the transition amplitude between occupied orbital \(\varphi_i\) and empty orbital \(\varphi_a\) of the \(n\)-th excitation.

The Tamm-Dancoff approximation (TDA) is also implemented, which constrains \(B=0\) and \(\mathbf{Y}^{(n)}=0\). In TDA, \(\Omega^{(n)}_\mathrm{TDA}\) and \(\mathbf{X}^{(n)}_\mathrm{TDA}\) can be computed from the Hermitian eigenvalue problem

\[ A \mathbf{X}^{(n)}_\mathrm{TDA} = \Omega^{(n)}_\mathrm{TDA} \mathbf{X}^{(n)}_\mathrm{TDA} \quad . \]

Diagonalizing \(A\) in TDA, or the full block-matrix \(ABBA\), takes in the order of \((N_\mathrm{occ} N_\mathrm{empty})^3\) floating point operations. This translates to a computational scaling of \(O(N^6)\) in the system size \(N\). Alternatively, the lowest excitations can be obtained iteratively with a block Davidson solver that applies \(A\) and \(B\) to trial vectors without forming them (BSE_DIAG_METHOD ITERDIAG, settings in BSE_ITERAT).

1.2 Optical absorption spectrum

The BSE further allows the investigation of optical properties. For example, the optical absorption spectrum can be computed as the imaginary part of the dynamical dipole polarizability tensor \(\alpha_{\mu,\mu'}(\omega) \) with \((\mu,\mu'\in\{x,y,z\})\):

\[ \alpha_{\mu,\mu'}(\omega) = - \sum_n \frac{2 \Omega^{(n)} d^{(n)}_{\mu} d^{(n)}_{\mu'}}{(\omega+i\eta)^2-\left(\Omega^{(n)}\right)^2} \quad , \]

where we have introduced an artificial broadening \(\eta\). The transition moments \(d^{(n)}_{\mu}\) are computed in the length gauge \((\mu\in\{x,y,z\})\) as

\[ d^{(n)}_{\mu} = \sqrt{2} \sum_{i,a} \langle \varphi_i|\hat{\mu}| \varphi_a \rangle (X_{ia}^{(n)} + Y_{ia}^{(n)}) \quad . \]

When the molecules are not aligned, e.g. for gas phase and liquids, the spatial average is sufficient, i.e. the optical absorption spectrum can be computed as

\[ \mathrm{Im}\left[\bar{\alpha}(\omega)\right] = \frac{1}{3} \sum_{\mu\in\{x,y,z\}} \mathrm{Im}\left[\alpha_{\mu,\mu}(\omega)\right] \quad . \]

We can rewrite the last equation as

\[ \mathrm{Im}\left[\bar{\alpha}(\omega)\right] = - \mathrm{Im}\left[ \sum_n \frac{f^{(n)}}{(\omega+i\eta)^2-\left(\Omega^{(n)}\right)^2} \right] \quad . \]

where we introduced the oscillator strengths \(f^{(n)}\), which are defined by

\[ f^{(n)} = \frac{2}{3} \Omega^{(n)} \sum_{\mu\in\{x,y,z\}} | d^{(n)}_{\mu} |^2 \quad . \]

Additionally, the photoabsorption cross section tensor

\[ \sigma_{\mu,\mu'}(\omega) = \frac{4 \pi \omega}{c} \mathrm{Im}\left[\alpha_{\mu,\mu'}(\omega) \right] \quad . \]

is printed, where \(c\) denotes the speed of light.

1.3 Visualizing excited states using Natural Transition Orbitals (NTOs)

In order to analyse the excitation wave function independent of a specific choice of the molecular orbitals \(\varphi_p(\mathbf{r})\), we can rewrite it as

\[ \Psi_\text{excitation}^{(n)}(\mathbf{r}_e,\mathbf{r}_h) = \sum_I {\lambda_I^{(n)}} \phi_I^{(n)}(\mathbf{r}_e) \chi_I^{(n)}(\mathbf{r}_h) \quad . \]

in terms of the natural transitions orbitals (NTOs) \(\phi_I^{(n)}(\mathbf{r}_e) \) and \(\chi_I^{(n)}(\mathbf{r}_h)\). Here, we introduce the idea of electron and holes: In the excitation process, an electron leaves a hole behind in the (formerly) occupied orbitals \(\chi_I^{(n)}(\mathbf{r}_h)\) and is afterwards located in the empty orbitals \(\phi_I^{(n)}(\mathbf{r}_e) \). We call \(\phi_I^{(n)}(\mathbf{r}_e) \) and \(\chi_I^{(n)}(\mathbf{r}_h)\) a NTO pair of electron and hole with a weight \(\lambda_I^{(n)}\), where \(I=1,\cdots,N_b\) with \(N_b=N_\mathrm{occ}+N_\mathrm{empty}\).

The weights \(\lambda_I^{(n)}\) are sorted in decreasing order, i.e.:

\[\lambda_1^{(n)} \ge \lambda_2^{(n)} \ge \cdots \ge 0\]

For many cases, we observe \(\lambda_1^{(n)} \gg \lambda_2^{(n)}\), which indicates a dominant single-particle excitation character for this \(n\) and allows a visual analysis of a small number of NTO pairs.

Assuming \(\lambda_1^{(n)} = 1\) and \(\lambda_{I\neq 1}=0\), the excitation wave function simply is given as a product

\[ \Psi_\text{excitation}^{(n)}(\mathbf{r}_e,\mathbf{r}_h) = \phi_1^{(n)}(\mathbf{r}_e) \chi_1^{(n)}(\mathbf{r}_h) \quad . \]

In this case, the electron is excited from the occupied NTO \(\chi_1^{(n)}(\mathbf{r}_h)\) to the empty state \(\phi_1^{(n)}(\mathbf{r}_e)\). This process leaves a hole at \(\mathbf{r}_h\) and creates an excited electron at \(\mathbf{r}_e\).

Computationally, the NTO pairs \(\phi_I^{(n)}(\mathbf{r}_e) \) and \(\chi_I^{(n)}(\mathbf{r}_h)\) as well as \(\lambda_I^{(n)}\) are obtained by a singular value decomposition of the transition density matrix

\[\begin{split} T = \begin{pmatrix} 0 & {X}^{(n)}\\ \left({Y}^{(n)} \right)^T & 0 \end{pmatrix} \quad , \end{split}\]

i.e.:

\[\begin{split} {T}^{(n)} &= {U}^{(n)} {\Lambda^{(n)}} \left({V}^{(n)}\right)^T \\ \phi_I^{(n)}(\mathbf{r}_e) &= \sum_{p=1}^{N_b} \varphi_p(\mathbf{r}_e) V_{p,I}^{(n)} \quad , \\ \chi_I^{(n)}(\mathbf{r}_h) &= \sum_{q=1}^{N_b} \varphi_q(\mathbf{r}_h) U_{q,I}^{(n)} \quad . \end{split}\]

1.4 Measures for the size of an excited state

In this subsection, we introduce several measures, which allow us to analyse the wave function of an excited state \(\Psi_\text{excitation}^{(n)}(\mathbf{r}_e,\mathbf{r}_h)\) and its spatial properties, following Ref. [Mewes2018].

To that end, we define the exciton expectation value with respect to a generic operator \(\hat{O}\) as

\[ {\langle \hat{O} \rangle}_\text{exc}^{(n)} = \frac{ \langle \Psi_\text{excitation}^{(n)} | \hat{O} | \Psi_\text{excitation}^{(n)}\rangle }{ \langle \Psi_\text{excitation}^{(n)} | \Psi_\text{excitation}^{(n)}\rangle } \quad , \]

where we drop the excitation index \(n\) from now on for better readability.

For each excitation level \(n\), we have then several quantities, which allow us to quantify the spatial extent of the electron, the hole and their combined two-particle character as combined electron-hole pair, i.e. as an exciton. By that, we can often determine the type of the excited state, i.e. distinguish between, e.g., valence, Rydberg or charge-transfer states.

First, we define the distance between electron and hole as

\[ d_{h \rightarrow e} = | {\langle \mathbf{r}_h - \mathbf{r}_e \rangle}_\mathrm{exc} | \quad , \]

which can be used to distinguish different classes of excitations: For example in a charge-transfer state, electron and hole sit on different parts of the molecule and therefore have a non-vanishing electron-hole distance \(d_{h \rightarrow e}\).

Further, we can measure the size of electron and hole, respectively, as

\[ \sigma_{e/h} = \sqrt{ {\langle \mathbf{r}_{e/h}^2 \rangle}_\mathrm{exc} - {\langle \mathbf{r}_{e/h} \rangle}_\mathrm{exc} ^2 } \quad , \]

which allow us to distinguish between Rydberg states, where \(\sigma_h \ll \sigma_e\), and valence states, where \(\sigma_h \approx \sigma_e\).

Closely related to these quantities, we can also define the exciton size

\[ d_\mathrm{exc} = \sqrt{ {\langle |\mathbf{r}_h - \mathbf{r}_e|^2 \rangle}_\mathrm{exc} } \quad . \]

which quantifies the spatial extent of the combined electron-hole pair. As one would expect, the exciton size \(d_\mathrm{exc}\) increases when \(d_{h \rightarrow e}\), \(\sigma_{e}\) or \(\sigma_{h}\) increase.

Finally, we quantify the correlation of electron and hole by the electron-hole correlation coefficient

\[ R_{eh} = \frac{1}{\sigma_e \sigma_h} \left( {\langle \mathbf{r}_h \cdot \mathbf{r}_e \rangle}_\mathrm{exc} - {\langle \mathbf{r}_h \rangle}_\mathrm{exc} \cdot {\langle \mathbf{r}_e \rangle}_\mathrm{exc} \right) \quad , \]

which allows us to distinguish between correlated (\(R_{eh}>0\)) motion, i.e. bound excitons, and anticorrelated (\(R_{eh}<0\)) motion, where electron and hole try to avoid each other.

2. BSE input

For starting a BSE calculation one needs to set the RUN_TYPE to ENERGY and the following sections for GW and BSE:

&GW
  SELF_CONSISTENCY      evGW0         ! We strongly recommend to use evGW0 (and PBE) for BSE runs
  &BSE  
    TDA                 TDA+ABBA      ! Diagonalizing ABBA and A
    SPIN_CONFIG         SINGLET       ! or TRIPLET
    NUM_PRINT_EXC       15            ! Number of printed excitations
    ENERGY_CUTOFF_OCC   -1            ! Set to positive numbers (eV) to
    ENERGY_CUTOFF_EMPTY -1            ! truncate matrices A_ia,jb and B_ia,jb
    NUM_PRINT_EXC_DESCR -1            ! Number of printed exciton descriptors, -1 to default to NUM_PRINT_EXC
    &BSE_SPECTRUM                     ! Activates computation and output of optical absorption spectrum
      ETA_LIST 0.01 0.02              ! Multiple broadenings can be specified within one run
    &END BSE_SPECTRUM
    &NTO_ANALYSIS
      STATE_LIST 1 2                  ! List of states for which NTOs are computed
      CUBE_FILES T                    ! Write cube files for NTOs
      STRIDE 6 6 6                    ! Coarse grid for smaller example cube files
    &END NTO_ANALYSIS
  &END BSE
&END GW

In the upper GW/BSE section, the following keywords have been used:

SELF_CONSISTENCY: Determines which GW self-consistency (G₀W₀, evGW₀ or evGW) is used to calculate the single-particle GW energies \(\varepsilon_p^{GW}\) needed in the BSE calculation. We recommend to always use evGW₀ for BSE runs. The screened coulomb interaction is always computed from the DFT level, i.e. \(W_0(\omega=0)\) is used for the \(ABBA\)-matrix, including BSE@evGW@DFT.
TDA: Three options available:
- ON diagonalize \(A\),
- OFF generalized diagonalization of \(ABBA\),
- TDA+ABBA CP2K diagonalizes \(ABBA\) as well as \(A\).
BSE_DIAG_METHOD: FULLDIAG (default) diagonalizes the matrices; ITERDIAG runs a block Davidson solver for the lowest NUM_EXC_EN excitations without forming \(A\) and \(B\), with the settings of BSE_ITERAT.
SPIN_CONFIG: Two options available: Choose between SINGLET for computing singlet excitation energies \((\alpha^S = 2)\) and TRIPLET for computing triplet excitation energies \((\alpha^T=0)\). Standard is SINGLET as an electronic excitation directly after photoexcitation is a singlet due to angular momentum conservation; triplet excited states can form by intersystem crossing.
NUM_PRINT_EXC: Number of excitations \(N_\text{exc}^\text{print}\) to be printed.
ENERGY_CUTOFF_OCC \(E_\text{cut}^\text{occ}\): Restrict occupied molecular orbital (MO) indices \(i\) and only use occupied MOs with \(\varepsilon_i\in[\varepsilon_{i=\text{HOMO}}^{DFT}-E_\text{cut}^\text{occ},\varepsilon_{i=\text{HOMO}}^{DFT}]\). The cutoff is applied to the DFT eigenvalues, also when the BSE takes quasiparticle energies from the GW. Setting a small ENERGY_CUTOFF_OCC drastically reduces the computation time and the memory consumption, but also might affect the computed excitation energies \(\Omega^{(n)}\). Recommended to use for large systems with more than 30 atoms, but we recommend a careful convergence test by increasing ENERGY_CUTOFF_OCC and observing the effect on \(\Omega^{(n)}\) [Graml2026].
ENERGY_CUTOFF_EMPTY \(E_\text{cut}^\text{empty}\): Analogous to ENERGY_CUTOFF_OCC, but for the empty states, i.e. only empty states in the interval \(\varepsilon_a\in[\varepsilon_{a=\text{LUMO}}^{DFT},\varepsilon_{a=\text{LUMO}}^{DFT}+E_\text{cut}^\text{empty}]\).
NUM_PRINT_EXC_DESCR: Number of excitations, for which the exciton descriptors are printed.
BSE_SPECTRUM: Activates computation and printing of the optical absorption spectrum. For each chosen option in TDA, one file is printed with columns frequency \(\omega\), \(\mathrm{Im}\left[\bar{\alpha}(\omega)\right]\) and the imaginary part of the elements of the dynamical dipole polarizability tensor \(\mathrm{Im}\left[{\alpha_{\mu,\mu'}}(\omega)\right]\) as well as another file with the respective entries for the photoabsorption cross section tensor \({\sigma_{\mu,\mu'}}(\omega)\). The frequency range, step size and the broadening \(\eta\) can be specified by the user (cf. keywords in BSE_SPECTRUM). Further, multiple broadenings \(\eta\) can be given for one cp2k run (ETA_LIST), resulting in one file per specified \(\eta\) and TDA (see Sec. 3.2).
NTO_ANALYSIS: Activates computation and printing of the natural transition orbitals. Within this section, the user can specify the NTO pairs, which should be printed, e.g. by the excitation indices of the desired levels \(n\), their minimum oscillator strength \(f^{(n)}\) or by the minimum weight of their NTO coefficient \(\lambda_I\). For each NTO pair \(\phi_I^{(n)}(\mathbf{r}_e) \) and \(\chi_I^{(n)}(\mathbf{r}_h)\) of a certain excitation \(n\), which fulfills these requirements, and the chosen options in TDA, the NTO coefficients \(\lambda_I\) are computed and the corresponding NTOs are printed to separate files if CUBE_FILES is activated.

The following settings from DFT will also have an influence on the BSE excitation energies and optical properties:

XC_FUNCTIONAL: Choose between one of the available xc-functionals. The starting point can have a profound influence on the excitation energies [Graml2026]. Motivated by the discussion in [Graml2026], we strongly recommend to use BSE@evGW₀@PBE, i.e. the PBE functional as DFT starting point (see also SELF_CONSISTENCY).
BASIS_SET: Specify the basis set, which affects \(N_\mathrm{empty}\) and thus the size of the matrices \(A_{ia,jb}\) and \(B_{ia,jb}\). The aug-cc-pVDZ basis set should be sufficient for most calculations, but needs to be checked regarding convergence, e.g. using aug-cc-pVTZ.

The memory consumption of the BSE algorithm is large, it is approximately \(100 \cdot N_\mathrm{occ}^2 N_\mathrm{empty}^2\) Bytes. You can see \(N_\mathrm{occ}\), \(N_\mathrm{empty}\) and the estimated memory consumption from the BSE output. The BSE implementation is well parallelized, i.e. you can use several nodes that can provide the memory. With BSE_DIAG_METHOD ITERDIAG these matrices are never allocated; the memory is then set by the RI three-center tensors, and the solver prints its own estimate.

We have benchmarked the numerical precision of our BSE implementation in [Graml2026] and compared its results to the BSE implementation in FHI aims [Liu2020]. For our recommended settings, i.e. BSE@evGW₀@PBE with the aug-cc-pVDZ basis set, we have found excellent agreement with less than 5 meV mean absolute deviation averaged over the first 10 excitation levels and the 28 molecules in Thiel’s set for Singlet excitations.

The current BSE implementation in CP2K works for molecules. The inclusion of periodic boundary conditions in a Γ-only approach and with full k-point sampling is work in progress.

3. Minimal example for a BSE calculation

3.1 Input file

In this section, we provide a minimal example of a BSE calculation on H₂. For the calculation you need the input file BSE_H2.inp and the aug-cc-pVDZ basis (Download).

Please copy both files into your working directory and run CP2K by

mpirun -n 1 cp2k.psmp BSE_H2.inp

which requires 5 GB RAM and takes roughly 45 seconds on 1 core. You can find the input and output file here. We use the basis sets aug-cc-pVDZ and aug-cc-pVDZ-RIFIT from the file BASIS-aug. These basis sets can be obtained from the Basis Set Exchange Library as aug-cc-pVDZ, aug-cc-pVDZ-RIFIT. The geometry for H₂ was taken from [vanSetten2015].

3.2 Output

The BSE calculation outputs the excitation energies \(\Omega^{(n)}\) for n = 1, …, \(N_\text{exc}^\text{print}\):

BSE| Excitation energies from solving the BSE without the TDA:
BSE|
BSE|     Excitation n   Spin Config       TDA/ABBA   Excitation energy Ω^n (eV)
BSE|                1       Singlet         -ABBA-                      12.0361
BSE|                2       Singlet         -ABBA-                      13.0327
BSE|                3       Singlet         -ABBA-                      15.8729
BSE|                4       Singlet         -ABBA-                      15.8729

Afterwards, the single-particle transitions, i.e. the eigenvector elements \(X_{ia}^{(n)}\) and \(Y_{ia}^{(n)}\), with \(|X_{ia}^{(n)}|\) > EPS_X or \(|Y_{ia}^{(n)}|\) > EPS_X, are printed. The third column is either

\(\Rightarrow\): for excitations, i.e. entries of X_ia^n, or
\(\Leftarrow\): for deexcitations, i.e. entries of Y_ia^n.

BSE| The following single-particle transitions have significant contributions,
BSE| i.e. |X_ia^n| > 0.100 or |Y_ia^n| > 0.100, respectively :
BSE|         -- Quick reminder: HOMO i =    1 and LUMO a =    2 --
BSE|
BSE| Excitation n           i =>/<=     a      TDA/ABBA       |X_ia^n|/|Y_ia^n|
BSE|
BSE|            1           1    =>     2        -ABBA-                  0.9457
BSE|            1           1    =>     4        -ABBA-                  0.3437
BSE|
BSE|            2           1    =>     3        -ABBA-                  0.9985
BSE|
BSE|            3           1    =>     5        -ABBA-                  1.0010
BSE|
BSE|            4           1    =>     6        -ABBA-                  1.0010

In the case of H₂, the lowest excitation n = 1 is mainly built up by a transition from the HOMO (i=1) to the LUMO (a=2), what is apparent from \(X_{i=\text{HOMO},a=\text{LUMO}}^{(n=1)}= 0.9457\), and also contains a considerable contribution from the 1=>4 (HOMO=>LUMO+2) transition (\(X_{i=\text{HOMO},a=\text{LUMO+2}}^{(n=1)}=0.3437\) ). The contributions larger than 1, e.g. \(X_{i=1,a=5}^{(n=3)}= 1.0010\) are due to the unusual normalization condition of the generalized eigenvalue problem for the \(ABBA\)-matrix:

\[ \sum_{i,a} X_{ia}^{(m)} X_{ia}^{(n)} - Y_{ia}^{(m)} Y_{ia}^{(n)} = \pm \delta_{mn} \quad . \]

In the next section, the optical properties, i.e. the transition moments \(d^{(n)}_{r}\) and the oscillator strengths \(f^{(n)}\), are printed:

BSE| Optical properties from solving the BSE without the TDA:
BSE|
BSE|  Excitation n  TDA/ABBA        d_x^n     d_y^n     d_z^n Osc. strength f^n
BSE|             1    -ABBA-        0.000    -0.000    -1.191             0.418
BSE|             2    -ABBA-       -0.000    -0.000     0.000             0.000
BSE|             3    -ABBA-       -0.829     0.721    -0.000             0.469
BSE|             4    -ABBA-       -0.721    -0.829     0.000             0.469

We can see that \(n=1,3,4\) are optically active since they have a nonzero oscillator strength \(f^{(n)}\).

In case BSE_SPECTRUM is invoked, the frequency-resolved optical absorption spectrum is printed in addition. For each specified \(\eta\) in ETA_LIST, a separate file, e.g. BSE-ABBA-eta=0.010.spectrum, is created, which starts with

# Imaginary part of polarizability α_{µ µ'}(ω) = -\sum_n [2 Ω^n d_µ^n d_µ'^n] / [(ω+iη)²- (Ω^n)²] from Bethe Salpeter equation for method -ABBA-
#     Frequency (eV)       Im{α_{avg}(ω)}          Im{α_xx(ω)}          Im{α_xy(ω)}          Im{α_xz(ω)}          Im{α_yx(ω)}          Im{α_yy(ω)}          Im{α_yz(ω)}          Im{α_zx(ω)}          Im{α_zy(ω)}          Im{α_zz(ω)}
          0.00000000           0.00000000           0.00000000           0.00000000           0.00000000           0.00000000           0.00000000           0.00000000           0.00000000           0.00000000           0.00000000
          0.10000000           0.00005283           0.00003298          -0.00000000          -0.00000000          -0.00000000           0.00003298           0.00000000          -0.00000000           0.00000000           0.00009251

The next table summarizes some of the exciton descriptors of each excitation level \(n\):

BSE| Exciton descriptors from solving the BSE   without the TDA:
BSE|
BSE|    n      c_n     d_eh [Å]     σ_e [Å]     σ_h [Å]   d_exc [Å]        R_eh
BSE|    1    1.026       0.0000      2.0542      0.8576      2.2521     -0.0332
BSE|    2    1.004       0.0000      2.9922      0.8578      3.1127     -0.0000
BSE|    3    1.005       0.0000      1.5859      0.8568      1.8134     -0.0145
BSE|    4    1.005       0.0000      1.5859      0.8568      1.8134     -0.0145

where we obtain the normalization factor \(c_n = \langle \Psi_\text{excitation}^{(n)} | \Psi_\text{excitation}^{(n)}\rangle \), the vectorial distance between electron and hole \(d_{e\rightarrow h}\), sizes of electron \(\sigma_e\) and hole \(\sigma_h\), the exciton size \(d_\mathrm{exc}\) and the electron-hole correlation coefficient \(R_{eh}\).

Further, the expectation values \(\langle \mathbf{r}_e \rangle_\mathrm{exc}\) and \(\langle \mathbf{r}_h \rangle_\mathrm{exc}\) as well as the reference point of origin and the atomic geometry are printed.

The last section summarizes the characteristics of the NTO analysis:

BSE| Number of excitations, for which NTOs are computed                       2
BSE| Threshold for oscillator strength f^n                                  ---
BSE| Threshold for NTO weights (λ_I)^2                                    0.010
BSE|
BSE| NTOs from solving the BSE without the TDA:
BSE|
BSE|  Excitation n  TDA/ABBA   Index of NTO I               NTO weights (λ_I)^2
BSE|
BSE|             1    -ABBA-                1                           1.01320
BSE|             1    -ABBA-                2                           0.01320
BSE|
BSE|             2    -ABBA-                1                           1.00210

In the input file, we have specified a list of three excitation levels, no constraint on the oscillator strengths \(f^{(n)}\) and used the default threshold on the weights, given by \(\lambda_I^2 \leq 0.010\), which is repeated at the beginning of the section. Therefore, we obtain two NTO pairs for \(n=1\) with significant weight and one NTO pair for \(n=2\). For \(n=1\), we can verify that the NTO weights satisfy \(1.01320 + 0.01320 \approx \sum_I {\left( \lambda_I^{(n)}\right)}^2 = c_n \approx 1.026\) when comparing to the output of the exciton descriptor section.

3.3 Large scale calculations

Going to larger systems is a challenge for a GW+BSE-calculation, since the memory consumption increases with \(N_\mathrm{occ}^2 N_\mathrm{empty}^2\). The used \(N_\mathrm{occ}\), \(N_\mathrm{empty}\) and the required memory of a calculation are printed in the output file to estimate the memory consumption. Here, you can find a sample output of a BSE@evGW0@PBE calculation on a nanographene with 206 atoms, which has a peak memory requirement of 2.5 TB RAM:

BSE| Cutoff occupied orbitals [eV]                                       80.000
BSE| Cutoff empty orbitals [eV]                                          10.000
BSE| First occupied index                                                   155
BSE| Last empty index (not MO index!)                                       324
BSE| Energy of first occupied index [eV]                                -25.080
BSE| Energy of last empty index [eV]                                      9.632
BSE| Energy difference of first occupied index to HOMO [eV]              20.289
BSE| Energy difference of last empty index to LUMO [eV]                  10.549
BSE| Number of GW-corrected occupied MOs                                    334
BSE| Number of GW-corrected empty MOs                                       324
BSE|
BSE| Total peak memory estimate from BSE [GB]                         8.725E+02
BSE| Peak memory estimate per MPI rank from BSE [GB]                      5.453

To enable BSE calculations on large molecules, we recommend to use large clusters with increased RAM and explicitly setting the keywords ENERGY_CUTOFF_OCC and ENERGY_CUTOFF_EMPTY, see details given above.

4. BSE on the GW of `PROPERTIES%BANDSTRUCTURE%GW` (`PROPERTIES%BSE`)

The BSE section exists in two places. Each runs on the GW next to it, the keywords of both are identical, and both treat non-periodic systems (molecules) only:

RI_RPA%GW%BSE, in full FORCE_EVAL%DFT%XC%WF_CORRELATION%RI_RPA%GW%BSE, runs on the \(O(N^4)\)-scaling GW of the RI_RPA%GW section, as in the sections above.
PROPERTIES%BSE, in full FORCE_EVAL%PROPERTIES%BSE, runs on the low-scaling GW of PROPERTIES%BANDSTRUCTURE%GW, which must be present in the same input and is described in Band structure from GW.

The GW of PROPERTIES%BANDSTRUCTURE%GW is designed for large systems [Graml2024], so PROPERTIES%BSE is the BSE section to use for large molecules. PROPERTIES%BSE has the following requirements of its own:

A non-periodic cell, i.e. PERIODIC NONE.
No RI_RS subsection in PROPERTIES%BANDSTRUCTURE%GW. By default, PROPERTIES%BANDSTRUCTURE%GW computes the three-center integrals \((\mu\nu|P)\) of the resolution of the identity (RI) between atomic orbitals \(\mu, \nu\) and RI basis functions \(P\) explicitly; RI_RS replaces them by a representation on a real-space grid. PROPERTIES%BSE builds its kernel from the explicit three-center integrals.
The Γ-point only, i.e. no BANDSTRUCTURE_PATH subsection and no DOS%KPOINTS in PROPERTIES%BANDSTRUCTURE.
No RTBSE section in the same run.
With NTO_ANALYSIS, CUBE_FILES F: PROPERTIES%BSE writes no cube files of the natural transition orbitals.

PROPERTIES%BSE takes the quasiparticle energies, the three-center integrals and the static screened interaction \(W(i\omega = 0)\) from PROPERTIES%BANDSTRUCTURE%GW; no second RPA calculation is run for the screening. \(W(i\omega = 0)\) is summed from the screened interaction on the imaginary-time grid of PROPERTIES%BANDSTRUCTURE%GW. The numerical parameters of PROPERTIES%BANDSTRUCTURE%GW are therefore convergence parameters of PROPERTIES%BSE as well, beyond their effect on the quasiparticle energies:

NUM_TIME_FREQ_POINTS sets the imaginary-time grid from which \(W(i\omega = 0)\) is summed.
CUTOFF_RADIUS_RI is the radius of the truncated Coulomb metric of the RI. PROPERTIES%BSE builds its kernel in the same metric.
EPS_FILTER filters the three-center integrals, which enter the screening and the kernel of PROPERTIES%BSE.

A minimal input:

&PROPERTIES
  &BANDSTRUCTURE
    &GW
      NUM_TIME_FREQ_POINTS 30
      CUTOFF_RADIUS_RI 20.0
      MEMORY_PER_PROC 2
    &END GW
  &END BANDSTRUCTURE
  &BSE
    BSE_DIAG_METHOD FULLDIAG
    ENERGY_CUTOFF_EMPTY 50.0
    ENERGY_CUTOFF_OCC 60.0
    TDA OFF
  &END BSE
&END PROPERTIES

A CUTOFF_RADIUS_RI beyond the extent of the molecule turns the truncated Coulomb metric of PROPERTIES%BANDSTRUCTURE%GW into the Coulomb metric of RI_RPA%GW. With such a radius and enough time and frequency points, the two GW implementations, and with them PROPERTIES%BSE and RI_RPA%GW%BSE, can be converged to each other. With looser settings, a residual deviation between the excitation energies of PROPERTIES%BSE and of RI_RPA%GW%BSE remains, which stems from the GW step.

PROPERTIES%BSE transforms the three-center integrals to molecular orbitals in batches of RI atoms that fit the memory MEMORY_PER_PROC leaves per MPI rank after what the rank holds; by default the BSE detects that memory as the GW does for its own MEMORY_PER_PROC. The same remainder bounds the subspace of the iterative solver, and with FULLDIAG the memory of the diagonalization is estimated before the GW starts and the run stops when it does not fit. The output reports the memory per rank, what is left for each step, the largest batch and the number of batches. The three-center integrals of the BSE stay until the BSE ends and take about \(8 N_\mathrm{RI} (N_\mathrm{occ}^2 + N_\mathrm{occ} N_\mathrm{virt} + N_\mathrm{virt}^2)\) bytes over all ranks and per spin, with \(N_\mathrm{RI}\) RI functions and \(N_\mathrm{occ}\), \(N_\mathrm{virt}\) levels in the active window; a value set by hand covers them and what SCF and GW hold.

The budget holds per step and is no total for the run. The transformation of the integrals and the solver run one after the other, and each may take up to the budget on top of what the run holds when the step starts. By default, each step takes half of the memory that is free at that moment, so no value has to be estimated. A budget set by hand should not exceed the memory per MPI rank minus what the run holds already: the memory that SCF and GW have not released, and the transformed three-center integrals, which stay until the BSE ends and take about \(8\,N_\text{RI}\,(N_\text{occ}^2 + N_\text{occ} N_\text{virt} + N_\text{virt}^2)\) bytes over all ranks and per spin, with \(N_\text{RI}\) RI functions and \(N_\text{occ}\), \(N_\text{virt}\) levels in the active window.