14.1 Single particle spectra based on the $GW$ approximation

14.1.1 Theoretical background

A method to systematically improve upon DFT estimates of single-particle excitation spectra, that is, ionization potentials and electron affinities, is the $GW$ method. Its central object is the single-particle Green’s function $G$; its poles describe single-particle excitation energies and lifetimes. In particular, the poles up to the Fermi-level correspond to the primary vertical ionization energies. The $GW$-approach is based on an exact representation of $G$ in terms of a power series of the screened Coulomb interaction $W$, which is called the Hedin equations. The $GW$-equations are obtained as an approximation to the Hedin-equations, in which the screened Coulomb interaction $W$ is calculated neglecting so called vertex corrections. In this approximation the self–energy $\Sigma$, which connects the fully interacting Green’s function $G$ to a reference non-interacting Green’s function $G_0$, is given by $\Sigma=GW$.

This approach can be used to perturbatively calculate corrections to the Kohn–Sham spectrum. To this end, the Green’s function is expressed in a spectral representation as a sum of quasi particle states. \[\begin{equation} G(\mathbf{r},\mathbf{r}';z)=\sum_n \frac{\Psi_{\text{r},n}(\mathbf{r},z)\Psi_{\text{l},n}^\dagger(\mathbf{r}',z)}{z - \varepsilon_n(z)+i\eta\,\mathrm{sgn}(\varepsilon_n-\mu)}. \end{equation}\](14.1) Under the approximation that the KS states are already a good approximation to these quasi–particle states $\Psi_{\text{l},n}$ the leading order correction can be calculated by solving the zeroth order quasi–particle equation: \[\begin{equation} \varepsilon_n = \epsilon_n +\langle n|\Sigma[G_{\text{KS}}](\varepsilon_n) -V_{\mathrm{xc}}|n\rangle \end{equation}\](14.2) An approximation to the solution of this equation can be obtained by linearizing it: \[\begin{equation} \varepsilon_n = \epsilon_n + Z_n\left<n\right\vert\Sigma(\epsilon_n)-V_{\mathrm{xc}}\left\vert n\right> \end{equation}\](14.3) here, $Z_n$ is given by: \[\begin{equation} Z_n = \left[1 - \left<n\right\vert\left.\frac{\partial\Sigma(E)}{\partial E}\right|_{E=\epsilon_n}\left\vert n\right>\right]^{-1} \end{equation}\](14.4) reducing the computational effort to a single iteration.

The self–energy $\Sigma$ appearing in Eqn. (14.2) is calculated in the $GW$ approximation from the KS Green’s function and screening. This is the so-called $G_0W_0$ approximation. The Self–energy splits in an energy independent exchange part $\Sigma^{\mathrm{x}}$ and a correlation part $\Sigma^{\mathrm{c}}(E)$ that does depend on energy. Their matrix elements are given by: \[\begin{align} \left<n\right\vert \Sigma^{\mathrm{x}}\left\vert n'\right>&= =-\sum_{i}\left(n i\vert i n'\right), \end{align}\](14.5) and \[\begin{align} \left<n\right\vert\Sigma^{\mathrm{c}}(\epsilon_n)\left\vert n\right>&= \sum_m\sum_{\underline{n}} \frac{\left|\left(\underline{n}n\vert\rho_{m}\right)\right|^2}{\epsilon_{n}-\epsilon_{\underline{n}}-Z_m\mathrm{sgn}(\epsilon_{\underline{n}}-\mu)}. \end{align}\](14.6) Where $Z_m = \Omega_m - i\eta$ are the excitation energies shifted infinitesimally into the complex plane. The $\rho_m$ are the corresponding excitation densities. More details, tests and benchmark calculations are can be found in Ref. .

14.1.2 $GW$ features

The $GW$ method is implemented in TURBOMOLE in the escf module supporting the following features:

LDA, GGA, meta-GGA and their Hybrid functionals can be used for the underlying DFT calculation.
Single-shot G$_0$W$_0$
Quasiparticle selfconsistent GW (qsGW)
Eigenvalue-only selfconsistent GW (evGW)
In $G_0W_0$, the linearized, Eqn. (14.3), and solved, Eqn. (14.2), quasiparticle equation.
Both RPA and TDDFT response functions can be used to screen the Coulomb interaction in constructing $W$, although we only recommend to use RPA response.
Closed shell and open shell references are supported in all calculations. Additionally, Kramers-restricted closed shell and also Kramers-unrestricted closed and open shell systems within the two-component relativistic framework (inclusion of spin-orbit coupling) can be treated using $soghf. For open shell two-component calculations collinear and non-collinear approaches are implemented.

14.1.3 General recipe for $G_0W_0$ calculations

Since TURBOMOLE 7.5 two fast $GW$ options are available: Just add $g0w0 OR $evgw to the control file and run escf with the -gw flag: escf -gw. This automatically sets all parameters to curated sensible values and will provide very good results and starting points for a BSE calculation in the majority of all cases.

The general recipe for a $G_0W_0$ and ev$GW$ calculations with dRPA response using RI is as follows:

define session
Choose reasonable cbas auxiliary basis sets (define)
dscf or ridft calculation
Provide $gw OR $rigw flags and keywords
To trigger the fast RI algorithm add $rick
escf calculation

Ad 1) The def2-TZVPP basis seems to be the most useful, it comes for all tested systems within 0.1 eV of the def2-QZVP result with about half the number of basis functions. A cbas must be set, the use of jbas-type fitting bases is discontinued for $GW$ since Turbomole 7.3. In the final define menu select the gw menu to set options for the actual $GW$ calculation.

The gw menu in define will set all needed variables for a gw/rigw calculation. The according $scfinstab and $soes flags will be set, and also the $gw or $rigw flags are written to the control file with the chosen options. Also the $rick flag is set. The control file provided by definecan usually be used for a calculation without further modifications; which are nevertheless described below of one wishes to modify it manually.

Ad 3) Symmetry up to D$_{\text{2h}}$ is available for both closed shell (rpas) and open shell system (urpa) for $GW$, for $gw and $rigw. Especially in the gw module the use of symmetry leads to large speedups and the user is encouraged to exploit symmetry if possible. Two-component calculations can only be performed in C$_1$. The $gw module can also handle open-shell 2c calculations, while $rigw is formally limited to Kramers-restricted closed shell molecules in the 2c case. Moreover, the simplified methods x$GW$ and s$GW$, which neglect contributions from excitation vectors are available for all point groups implemented in Turbomole.

Ad 4) In the last define step the gw menu can be selected to set up a $G_0W_0$ (and also $evGW$ and $qsGW$) calculation. There are three distinct versions available for $G_0W_0$ and ev$GW$:

$gw uses spectral representations just as the previous version, which scales as $N^6$ with system size. This version is to be preferred when full quasiparticle spectra are desired (i.e.: QP energies are also needed for non-valence orbitals, or generally for all). It is compatible with nearly all available starting points:
1. scalar/non-relativistic closed-shell systems (up to D$_{2h}$ symmetry)
2. scalar/non-relativistic open-shell systems (up to D$_{2h}$ symmetry)
3. two-component (2c) Kramers-restricted systems (C$_1$ only)
4. two-component (2c) open-shell systems (C$_1$ only)
$rigw; RI-AC-$G_0W_0$ and RI-AC-ev$GW$ variants construct the self-energy ${\Sigma}_C$ on the imaginary axis using numerical integration. The obtained result is then analytically continued to the real axis using Pade approximants. This ansatz scales roughly as $N^4$, and yields reliable results for valence orbitals, especially for HOMO and LUMO quasiparticle energies. RI-AC-$GW$ is the default variant for $rigw calculations. It is available for the following starting points:
1. scalar/non-relativistic closed-shell systems (up to D$_{2h}$ symmetry)
2. scalar/non-relativistic open-shell systems (up to D$_{2h}$ symmetry)
3. two-component (2c) Kramers-restricted systems (C$_1$ only)
4. two-component (2c) open-shell systems (C$_1$ only)
$rigw; RI-CD-$G_0W_0$ and RI-CD-ev$GW$ variants construct the self-energy ${\Sigma}_C$ using contour deformation. This ansatz scales roughly as $N^4$ (valence orbitals) - $N^5$ (core orbitals), and yields reliable results for most orbitals. RI-CD-$GW$ is invoked by the $contour$ keyword within the $rigw datagroup. It is available for the following starting points:
1. scalar/non-relativistic closed-shell systems (up to D$_{2h}$ symmetry))
2. scalar/non-relativistic open-shell systems (up to D$_{2h}$ symmetry))
3. two-component (2c) Kramers-restricted systems (C$_1$ only)
4. two-component (2c) open-shell systems (C$_1$ only)

Especially for 2c calculations the prefactor of 256 for a $gw calculation is reduced to 4-8 in a $rigw calculation, making calculations feasible also for large systems. This prefactor reduction is valid for both $rigw variants, RI-AC-$GW$ and RI-CD-$GW$. Magnetic fields and other Kramers-unristricted references can be used throughout all $GW$ and BSE calculations since Turbomole 7.7.

In cases where the number of states between a given quasiparticle energy and the gap becomes large or if many quasiparticle energies are to be optimized, the RI-CD-$GW$ approach may be approximated, by using a sampling strategy for the frequencies required in the calculation of the dielectric function and residues entering the contour. This approach follows the ideas described in Ref. [314], which employs the analytic continuation approach for the evaluation of most residues. The frequency grid is build by sampling the exact set of frequencies required for RI-CD-$GW$, where a frequency is removed if it is close to an already included grid point. Note that calculations within the generalized two-component framework are limited to the Kramers-symmetric systems using this approximation.

Performance & Accuracy: Since Turbomole 7.3 specialized $GW$ algorithms are implemented in escf, leading to large speedups. The analytic continuation ($rigw) variant usually yields self-energies with $meV$ accuracy compared to standard $gw for HOMO/LUMO orbitals, and is much more conservative with memory and CPU requirements. For two-component calculations (including spin-orbit interactions) also the prefactor of $G_0W_0$ and $evGW$ is largely reduced in $rigw. The standard $gw variant is however accurate and reliable throughout all quasiparticle states, also describing non-valence orbitals well and therefore recommended whenever possible.

Possible source of errors: The DFT options used in dscf or ridft should not be altered before starting escf. Otherwise erratic quasiparticle energies may be obtained. Also any files containing excitations (sing, unrs) from other escf runs should be removed prior to the $GW$ run. Also a cbas must be set before starting escf.

14.1 Single particle spectra based on the \(GW\) approximation

14.1.1 Theoretical background

14.1.2 \(GW\) features

14.1.3 General recipe for \(G_0W_0\) calculations