Category:Bethe-Salpeter equations

The formalism of the Bethe-Salpeter equation (BSE) allows for calculating the polarizability with the electron-hole interaction and constitutes the state of the art for calculating absorption spectra in solids.

Theory

Bethe-Salpeter equation

In the BSE, the excitation energies correspond to the eigenvalues ${\displaystyle \omega_\lambda}$ of the following linear problem^[1]

${\displaystyle \left(\begin{array}{cc} \mathbf{A} & \mathbf{B} \\ \mathbf{B}^* & \mathbf{A}^* \end{array}\right)\left(\begin{array}{l} \mathbf{X}_\lambda \\ \mathbf{Y}_\lambda \end{array}\right)=\omega_\lambda\left(\begin{array}{cc} \mathbf{1} & \mathbf{0} \\ \mathbf{0} & -\mathbf{1} \end{array}\right)\left(\begin{array}{l} \mathbf{X}_\lambda \\ \mathbf{Y}_\lambda \end{array}\right)~. }$

The matrices ${\displaystyle A}$ and ${\displaystyle A^*}$ describe the resonant and anti-resonant transitions between the occupied ${\displaystyle v,v'}$ and unoccupied ${\displaystyle c,c'}$ states

${\displaystyle A_{vc}^{v'c'} = (\varepsilon_v-\varepsilon_c)\delta_{vv'}\delta_{cc'} + \langle cv'|V|vc'\rangle - \langle cv'|W|c'v\rangle. }$

The energies and orbitals of these states are usually obtained in a ${\displaystyle G_0W_0}$ calculation, but DFT and Hybrid functional calculations can be used as well. The electron-electron interaction and electron-hole interaction are described via the bare Coulomb ${\displaystyle V}$ and the screened potential ${\displaystyle W}$.

The coupling between resonant and anti-resonant terms is described via terms ${\displaystyle B}$ and ${\displaystyle B^*}$

${\displaystyle B_{vc}^{v'c'} = \langle vv'|V|cc'\rangle - \langle vv'|W|c'c\rangle. }$

Due to the presence of this coupling, the Bethe-Salpeter Hamiltonian is non-Hermitian.

Tamm-Dancoff approximation

A common approximation to the BSE is the Tamm-Dancoff approximation (TDA), which neglects the coupling between resonant and anti-resonant terms, i.e., ${\displaystyle B}$ and ${\displaystyle B^*}$. Hence, the TDA reduces the BSE to a Hermitian problem

${\displaystyle AX_\lambda=\omega_\lambda X_\lambda~. }$

In reciprocal space, the matrix ${\displaystyle A}$ is written as

${\displaystyle A_{vc\mathbf{k}}^{v'c'\mathbf{k}'} = (\varepsilon_v-\varepsilon_c)\delta_{vv'}\delta_{cc'}\delta_{\mathbf{kk}'}+ \frac{2}{\Omega}\sum_{\mathbf{G}\neq0}\bar{V}_\mathbf{G}(\mathbf{q})\langle c\mathbf{k}|e^{i\mathbf{Gr}}|v\mathbf{k}\rangle\langle v'\mathbf{k}'|e^{-i\mathbf{Gr}}|c'\mathbf{k}'\rangle -\frac{2}{\Omega}\sum_{\mathbf{G,G}'}W_{\mathbf{G,G}'}(\mathbf{q},\omega)\delta_{\mathbf{q,k-k}'} \langle c\mathbf{k}|e^{i(\mathbf{q+G})}|c'\mathbf{k}'\rangle \langle v'\mathbf{k}'|e^{-i(\mathbf{q+G})}|v\mathbf{k}\rangle, }$

where ${\displaystyle \Omega}$ is the cell volume, ${\displaystyle \bar{V}}$ is the bare Coulomb potential without the long-range part

${\displaystyle \bar{V}_{\mathbf{G}}(\mathbf{q})=\begin{cases} 0 & \text { if } G=0 \\ V_{\mathbf{G}}(\mathbf{q})=\frac{4 \pi}{|\mathbf{q}+\mathbf{G}|^2} & \text { else } \end{cases}~, }$

and the screened Coulomb potential ${\displaystyle W_{\mathbf{G}, \mathbf{G}^{\prime}}(\mathbf{q}, \omega)=\frac{4 \pi \epsilon_{\mathbf{G}, \mathbf{G}^{\prime}}^{-1}(\mathbf{q}, \omega)}{|\mathbf{q}+\mathbf{G}|\left|\mathbf{q}+\mathbf{G}^{\prime}\right|}. }$

Here, the dielectric function ${\displaystyle \epsilon_\mathbf{G,G'}(\mathbf{q})}$ describes the screening in ${\displaystyle W}$ within the random-phase approximation (RPA)

${\displaystyle \epsilon_{\mathbf{G}, \mathbf{G}^{\prime}}^{-1}(\mathbf{q}, \omega)=\delta_{\mathbf{G}, \mathbf{G}^{\prime}}+\frac{4 \pi}{|\mathbf{q}+\mathbf{G}|^2} \chi_{\mathbf{G}, \mathbf{G}^{\prime}}^{\mathrm{RPA}}(\mathbf{q}, \omega). }$

Although the dielectric function is frequency-dependent, the static approximation ${\displaystyle W_{\mathbf{G}, \mathbf{G}^{\prime}}(\mathbf{q}, \omega=0)}$ is considered a standard for practical BSE calculations.

Scaling

The scaling of the BSE equation strongly limits its application for large systems. The main limiting factor is the diagonalization of the BSE Hamiltonian. The rank of the Hamiltonian is

${\displaystyle N_{\rm rank} = N_k\times N_c\times N_v}$,

where ${\displaystyle N_k}$ is the number of k-points in the Brillouin zone and ${\displaystyle N_c}$ and ${\displaystyle N_v}$ are the number of conduction and valence bands, respectively. The diagonalization of the matrix scales cubically with the matrix rank, i.e., ${\displaystyle N_{\rm rank}^3}$.

Despite the fact that this matrix diagonalization is usually the bottleneck for bigger systems, the construction of the BSE Hamiltonian also scales unfavorably and can play a dominant role in big systems, i.e.,

${\displaystyle N_k\times N_q\times (N_v\times N_v\times N_G\times N_c\times N_c)}$,

where ${\displaystyle N_q}$ is the number of q-points and ${\displaystyle N_G}$ number of G-vectors.

Exact diagonalization

The diagonalization of the BSE Hamiltonian can be perform using various eigensolvers provided in ScaLAPACK, ELPA, and cuSolver libraries. The advantage of this approach is that the eigenvectors can be directly obtained and used for the analysis of the excitons. Using the eigenvalues ${\displaystyle \omega_\lambda}$ and eigenvectors ${\displaystyle X_\lambda}$ of the BSE Hamiltonian, the macroscopic dielectric which accounts for the excitonic effects can be found

${\displaystyle \epsilon_M(\mathbf{q},\omega)= 1+2\lim_{\mathbf{q}\rightarrow 0}v(q)\sum_{\lambda} \left|\sum_{c,v,k}\langle c\mathbf{k}|e^{i\mathbf{qr}}|v\mathbf{k}\rangle X_\lambda^{cv\mathbf{k}}\right|^2 \times \left(\frac{1}{\omega_\lambda - \omega - i\delta}\right)~. }$

The following features are currently supported:

• Calculating the dielectric function and eigenvectors
• Calculations beyond Tamm-Dancoff approximation
• Calculations of ${\displaystyle \varepsilon(\mathbf{q},\omega)}$ for ${\displaystyle \mathbf{q}\neq0}$
• Fatband plot

Time evolution

Alternatively, it is possible to use the time-evolution algorithm which applies a short Dirac delta pulse of electric field and then follows the evolution of the dipole moments. The dielectric function is found via a Fourier transform ^[2]

${\displaystyle \epsilon_M(\omega)=1-\frac{4\pi}{\Omega}\int_0^{\infty} \mathrm{d} t \sum_{c,v,\mathbf{k}}\left(\langle\mu_{cv\mathbf{k}}| \xi_{cv\mathbf{k}}(t)\rangle+c . c .\right) e^{-i(\omega-i \delta) t} }$,

where ${\displaystyle \mu}$ and ${\displaystyle \xi(t)}$ are the dipole moments.

The solution found this way is strictly equivalent to the same solution as the exact diagonalization and can be used for obtaining the absorption spectrum, but does not yield the eigenvectors, which can be limiting for the analysis of the excitons. The advantage of this approach is the quadratic scaling with the size of the BSE Hamiltonian ${\displaystyle N_{rank}^2}$.

The time-evolution algorithm can be selected by setting IBSE = 1 in a BSE calculation. The required number of steps in the time-evolution calculation depends on the broadening CSHIFT and the maximum energy OMEGAMAX. The precision can be selected via tag BSEPREC.

Mind: The required number of steps does not depend on the size of the Hamiltonian

The following features are currently supported:

• Calculating the dielectric function
• Calculations beyond the Tamm-Dancoff approximation

Performing BSE calculations on GPU

As of VASP 6.5, the BSE and TDDFT (Casida) calculations can be fully run on NVIDIA GPUs. To be able to offload the BSE calculations to GPUs one has to compile VASP with the cuSOLVERMp and cuBLASMp libraries provided with NVHPC-SDK 24.7 or newer. The easiest way to add these libraries in your environment is loading the NVHPC-SDK module:

module load nvhpc-hpcx-cuda12/24.9

To enable these libraries in VASP, make sure to include the following lines in your makefile.include

CPP_OPTIONS+= -DCUSOLVERMP -DCUBLASMP
LLIBS      += -cudalib=cusolvermp,cublasmp -lnvhpcwrapcal

To be able to perform the BSE calculation on GPUs, VASP needs to store the full BSE Hamiltonian in the GPU memory, which is often the limiting factor. The memory required to store the BSE Hamiltonian can be estimated as ${\displaystyle N_{\rm rank}^2\times 16\cdot 10^{-9}}$ in Gb for ANTIRES = 0. In the case of exact diagonalization IBSE = 2, the eigensolver requires an additional scratch space.

Mind: In NVHPC-SDEK 24.7 one has to set the following variable to avoid hangups in the eigensolver export CAL_UCC_TLS=^nccl

How to

• Practical guide for solving the Bethe-Salpeter equation via diagonalization BSE calculations
• Practical guide for solving the Casida equation via diagonalization TDDFT calculations

References