SternheimerGW: a program for calculating GW quasiparticle band structures and spectral functions without unoccupied states

TL;DR

SternheimerGW is a computational tool that calculates GW quasiparticle band structures and spectral functions without requiring unoccupied states, using a linear response approach for improved accuracy and efficiency.

Contribution

It introduces a method based on Sternheimer equations to compute GW properties without summing over unoccupied states, enhancing computational efficiency and applicability.

Findings

Accurate spectral functions for solids achieved without unoccupied states.

Efficient parallel implementation integrated with Quantum Espresso.

Applicable to indirect band gap semiconductors and angle-resolved spectra.

Abstract

The SternheimerGW software uses time-dependent density-functional perturbation theory to evaluate GW quasiparticle band structures and spectral functions for solids. Both the Green's function G and the screened Coulomb interaction W are obtained by solving linear Sternheimer equations, thus overcoming the need for a summation over unoccupied states. The code targets the calculation of accurate spectral properties by convoluting G and W using a full frequency integration. The linear response approach allows users to evaluate the spectral function at arbitrary electron wavevectors, which is particularly useful for indirect band gap semiconductors and for simulations of angle-resolved photoelectron spectra. The software is parallelized efficiently, integrates with version 6.3 of Quantum Espresso, and is continuously monitored for stability using a test farm.

Tables

Table 1. Table 1: Combined code coverage of integration and unit tests, i.e., the fraction of F90 lines and subroutines of the code that are executed at least once by the test set. The higher tiers are covered exclusively by the integration tests; for the lower tiers unit tests play a large role.

	main	phys	algo	data
F90 lines
integration tests	100.0%	92.1%	70.3%	31.4%
unit tests	0.0%	0.0%	41.3%	74.3%
combined	100.0%	92.1%	86.7%	92.3%
subroutines
integration tests	100.0%	97.3%	64.6%	30.0%
unit tests	0.0%	0.0%	53.3%	80.0%
combined	100.0%	97.3%	87.7%	93.3%

Table 2. Table 2: Comparison of the band gap of a selected set of materials with reference results from the literature. All values are given in eV. Numerical parameters used for the calculations are listed in Table  3 .

	Present	Ref. 114	Ref. 115	Ref. 116	Ref. 117
BN	6.27	6.10	6.20		6.19
C	5.53	5.50	5.62	5.62	5.59
CdS	2.20	2.06	2.18		2.31
GaAs	1.23	1.30	1.31	1.21	1.10
Ge	0.78		0.75	0.63	0.50
MgO	7.13	7.25	7.17	6.71	7.08
Si	1.17	1.12	1.11	1.26	1.17

Table 3. Table 3: Convergence parameters used in the calculation of the G ​ W 𝐺 𝑊 GW corrections reported in Table 2 : The pseudopotentials (PP) are from the SG15 library 118 or the stringent ones from pseudo-dojo.org ; 119 the energy cutoffs are employed for the DFT self-consistent field cycle (scf), the exchange (x) or the correlation (c) self energy; the integration grids for the self-energy and the density response are homogeneous grids with N q subscript 𝑁 𝑞 N_{q} and N k subscript 𝑁 𝑘 N_{k} points; the coarse frequency mesh for W 𝑊 W are obtained from a ( N ω − 1 ) subscript 𝑁 𝜔 1 (N_{\omega}-1) -node Gauss-Laguerre quadrature along the imaginary axis including the origin, and extending up to ℏ ​ ω ^ Planck-constant-over-2-pi ^ 𝜔 \hbar\hat{\omega} ; and the dense mesh for W 𝑊 W as well as the coarse one for Σ Σ \Sigma are automatically constructed by the code. The row with the dagger contains the loose values used for the convergence test and the row with the asterisk is used for the convolution along the real frequency axis.

		energy cutoffs (Ry)					coarse, $W$		dense, $W$		coarse, $\Sigma$		threshold
material	PP	$E_{\text{cut}}^{\text{scf}}$	$E_{\text{cut}}^{\text{x}}$	$E_{\text{cut}}^{\text{c}}$	$N_q^{1/3}$	$N_k^{1/3}$	$N_{\omega }$	$\hslash \hat{\omega }$ (eV)	$N_{\omega }$	$\hslash \hat{\omega }$ (eV)	$N_{\omega }$	$\hslash \hat{\omega }$ (eV)	$t_W$	$t_G$
BN	SG15	60	50	40	8	8	31	104	80	320	20	150	${10}^{-4}$	${10}^{-5}$
C	SG15	70	50	24	8	8	31	104	50	280	30	180	${10}^{-4}$	${10}^{-5}$
CdS	Dojo	60	55	25	8	8	31	104	150	150	30	180	${10}^{-4}$	${10}^{-5}$
GaAs	Dojo	90	70	40	8	8	31	104	130	130	20	150	${10}^{-4}$	${10}^{-5}$
Ge	Dojo	90	70	25	8	8	31	104	130	130	20	150	${10}^{-4}$	${10}^{-5}$
MgO	SG15	70	65	30	6	8	31	104	50	280	30	180	${10}^{-4}$	${10}^{-5}$
Si	SG15	20	19	15	8	8	31	104	50	280	40	200	${10}^{-4}$	${10}^{-6}$
BN†	SG15	60	30	10	1	4	21	67	50	50	2	100	${10}^{-4}$	${10}^{-5}$
Si∗	SG15	16	15	12	8	8	201	768	800	100	351	35	${10}^{-10}$	${10}^{-5}$

Equations

$$\phi_{n\boldsymbol{k}}(\boldsymbol{r})=\frac{1}{\sqrt{\Omega}}\sum_{\boldsymbol{G}}u_{n\boldsymbol{k},\boldsymbol{G}}\exp\bigl{[}\text{i}(\boldsymbol{k}+\boldsymbol{G})\cdot\boldsymbol{r}\bigr{]}~{},$$

$${\forall }_{\mathbf{k},{\mathbf{G}}^{\prime }}{\mathcal{L}}_{\mathbf{k}}^{\text{G}}(\omega )G_{\mathbf{G}{\mathbf{G}}^{\prime }}(\omega )=-{\delta }_{\mathbf{G}{\mathbf{G}}^{\prime }}.$$

$${\mathcal{L}}_{\mathbf{k}}^{\text{G}}(\omega )={\hat{H}}_{\mathbf{k}}-\hslash \omega ,$$

$${ϵ}_{\mathbf{q},\mathbf{G}{\mathbf{G}}^{\prime }}(\omega )={\delta }_{\mathbf{G}{\mathbf{G}}^{\prime }}-\Delta n_{\mathbf{q},\mathbf{G}{\mathbf{G}}^{\prime }}(\omega ),$$

$$\Delta n_{\mathbf{q},\mathbf{r}{\mathbf{G}}^{\prime }}(\omega )=2\frac{1}{N_{\mathbf{k}}}\sum _{n\mathbf{k},\pm }^{\text{occ.}}{u}_{n\mathbf{k},\mathbf{r}}^{\ast }\Delta u_{n\mathbf{k}+\mathbf{q},\mathbf{r}{\mathbf{G}}^{\prime }}(\pm \omega ),$$

$$\forall_{n\boldsymbol{k},\boldsymbol{G}}\quad\mathcal{L}^{\text{C}}_{\boldsymbol{k}+\boldsymbol{q}}\Bigl{(}\frac{\varepsilon_{n\boldsymbol{k}+\boldsymbol{q}}}{\hbar}+\omega\Bigr{)}\Delta u_{n\boldsymbol{k}+\boldsymbol{q},\boldsymbol{G}\boldsymbol{G}^{\prime}}(\omega)=\delta u_{n\boldsymbol{k},\boldsymbol{G}\boldsymbol{G}^{\prime}}~{}.$$

$${\mathcal{L}}_{\mathbf{k}}^{\text{C}}(\omega )={\hat{H}}_{\mathbf{k}}-\hslash \omega +{\hat{P}}^{\text{occ}},$$

$$\delta u_{n\boldsymbol{k},\boldsymbol{G}\boldsymbol{r}^{\prime}}=\bigl{(}1-\hat{P}^{\text{occ}}\bigr{)}e^{i\boldsymbol{G}\cdot\boldsymbol{r}^{\prime}}u_{n\boldsymbol{k},\boldsymbol{r}^{\prime}},$$

$$\bar{W}_{\boldsymbol{q},\boldsymbol{G}\boldsymbol{G}^{\prime}}(\omega)=v_{\boldsymbol{q}+\boldsymbol{G}}~{}\bigl{(}\epsilon_{\boldsymbol{q},\boldsymbol{G}\boldsymbol{G}^{\prime}}^{-1}(\omega)-\delta_{\boldsymbol{G}\boldsymbol{G}^{\prime}}\bigr{)}~{},$$

$${\Sigma }_{\mathbf{k},\mathbf{r}{\mathbf{r}}^{\prime }}^{\text{c}}(\omega )=\sum _{{\omega }^{\prime }\mathbf{q}}\frac{\alpha w_{{\omega }^{\prime }}{w}_{\mathbf{q}}}{2\pi }{G}_{\mathbf{k}-\mathbf{q},\mathbf{r}{\mathbf{r}}^{\prime }}({\omega }^{\prime }){\bar{W}}_{\mathbf{q},\mathbf{r}{\mathbf{r}}^{\prime }}({\omega }^{\prime }-\omega ).$$

$${\Sigma }_{\mathbf{k},\mathbf{G}{\mathbf{G}}^{\prime }}^{\text{x}}=-\sum _{n\mathbf{q},{\mathbf{G}}^{\prime \prime }}^{\text{occ.}}\frac{w_{n\mathbf{k}-\mathbf{q}}}{\Omega }{u}_{n\mathbf{k}-\mathbf{q},{\mathbf{G}}^{\prime }-{\mathbf{G}}^{\prime \prime }}^{\ast }{u}_{n\mathbf{k}-\mathbf{q},\mathbf{G}-{\mathbf{G}}^{\prime \prime }}{v}_{\mathbf{q}+{\mathbf{G}}^{\prime \prime }},$$

$$v_{\boldsymbol{q}}=\frac{e^{2}}{4\pi\epsilon_{0}}\frac{4\pi}{q^{2}}\bigl{[}1-\cos(qR_{\text{cut}})\bigr{]}~{}.$$

$$v_{\boldsymbol{q}}=\frac{e^{2}}{4\pi\epsilon_{0}}\frac{4\pi}{q^{2}}\Bigl{[}1-\exp(\sqrt{q_{x}^{2}+q_{y}^{2}}z_{\text{cut}})\Bigr{]}\cos(q_{z}z_{\text{cut}})~{}.$$

$$A_{\boldsymbol{k},\boldsymbol{G}\boldsymbol{G}^{\prime}}(\omega)=-\frac{1}{\pi}\text{Im}\bigl{(}G_{\boldsymbol{k},\boldsymbol{G}\boldsymbol{G}^{\prime}}(\omega)\bigr{)},$$

$$A_{\boldsymbol{k},\boldsymbol{G}\boldsymbol{G}^{\prime}}(\omega)=-\frac{1}{\pi}\text{Im}\Bigl{(}\hbar\omega-H_{\boldsymbol{k}}-\Sigma^{\text{xc}}_{\boldsymbol{k}}(\omega)+V_{\boldsymbol{k}}^{\text{xc}}\Bigr{)}^{-1}_{\boldsymbol{G}\boldsymbol{G}^{\prime}}~{}.$$

$$A_{n\boldsymbol{k}}(\omega)=\frac{-\text{Im}{\Sigma^{\text{c}}_{n\boldsymbol{k}}(\omega)}}{\pi\Bigl{\lvert}\hbar\omega-\varepsilon_{n\boldsymbol{k}}-\Sigma^{\text{xc}}_{n\boldsymbol{k}}(\omega)+V_{n\boldsymbol{k}}^{\text{xc}}\Bigr{\rvert}^{2}}~{}.$$

$$\varepsilon_{n\boldsymbol{k}}^{\text{QP}}=\varepsilon_{nk}+Z_{n\boldsymbol{k}}\Bigl{(}\text{Re}\Sigma^{\text{xc}}_{n\boldsymbol{k}}(\varepsilon_{nk}/\hbar)-V_{n\boldsymbol{k}}^{\text{xc}}\Bigr{)}~{},$$

Full text

††journal: Computer Physics Communications

Program Summary

Program Title: SternheimerGW

Licensing provisions: GNU General Public License v3.0

Programming language: Fortran 2003

Nature of problem: The lack of the exchange-correlation discontinuity in density-functional theory (DFT) leads to a systematic underestimation of the band gap between conduction and valence states. Many-body perturbation theory in the $GW$ approximation provides an effective solution to this problem, as well as other limitations faced by DFT in the description of electronic excitations. However, the $GW$ method comes with its own set of limitations: (i) The perturbation of the system is typically expressed in terms of the unoccupied states, and achieving numerical convergence with respect to these states is often cumbersome. Since the underlying DFT codes rely only on the occupied states, their default behavior is often ill-suited to provide a sufficient amount of empty states. Combined with the large number of parameters that need to be converged, this limits the use of $GW$ codes by non-expert users and automatic scripts. (ii) Currently, $GW$ codes require that a homogeneous $\boldsymbol{k}$-point mesh is used in the calculation. Hence, features close to the band edges are often only accessible via prohibitively expensive dense $\boldsymbol{k}$-point meshes or interpolation techniques. (iii) To evaluate the frequency convolution of the Green’s function $G$ and the screened Coulomb interaction $W$, many current $GW$ calculations rely on approximations such as the plasmon-pole approximation, the analytic continuation, or the contour deformation. The relative merits and accuracy of the various approximations are not fully understood.

Solution method: In SternheimerGW, we address (i) by replacing the summation over unoccupied states with a linear response equation. The solution depends on the occupied states only, and it employs linear response algorithms already provided by the Quantum ESPRESSO suite to compute phonons and related properties. As an additional benefit, transforming the problem in a linear response equation removes the restriction to homogeneous $\boldsymbol{k}$-point meshes (ii), so that the $GW$ self-energy for any arbitrary point can be evaluated. Finally (iii), we provide the possibility to perform a full frequency convolution along the real frequency axis. This feature can serve as a benchmark for approximate integrations using models or analytic continuation.

1 Introduction

Kohn-Sham density function theory (DFT)3, 4, 5 is a powerful and extremely popular formalism for studying the total energy of an interacting electron system in its ground state. When used in the study of electronic excitations, such as the calculation of electron band structures and wavefunctions, DFT exhibits some well known deficiencies, for example the lack of the exchange-correlation discontinuity in the exchange-correlation potential,6, 7 which leads to an underestimation of the band gap in insulators, and the excessive delocalization of $d$ and $f$ electrons. Applying many-body perturbation theory corrections in the $GW$ approximation8, 9 allows one to correct some of these deficiencies using the non-local and frequency-dependent electron self-energy $\Sigma(\omega)$. Early numerical implementations of this method demonstrated an improved band gap for diamond10, 11 and similar improvement were subsequently shown for other typical semiconductors.12, 13 Since then, $GW$ has been used to accurately describe the band gaps for a variety of systems ranging from solids to interfaces and molecules.14, 15, 16, 17 It yields accurate band offsets18, defect energies,19 and improved effective masses.20 Recent developments focused on the self-consistency of the $GW$ method for molecules21, 22, 23 and solids,24 as well as total energy calculations.25, 26

More recent developments of the $GW$ method relate to photoelectron spectroscopy. Progress in high-energy resolution angle-resolved photoelectron spectroscopy (ARPES) had led to a renewed interest in the spectral function of materials, in particular the signatures of electron correlations and bosonic excitations. These features were originally examined for the homogeneous electron gas,27, 28, 9, 29, 30 but have recently been investigated both experimentally and theoretically for charge carriers in semiconductors coupling to plasmons,31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 1, 43 Fröhlich polarons44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55 or hybridizations of these excitations.56, 57 While the $GW$ approximation does produce satellite features in the spectral function, their strength is overestimated and their energy is blue-shifted.27, 28, 29, 58, 32 This mismatch is intrinsic to the $GW$ approximation and cannot be overcome by self-consistency.59, 60, 23, 41 Nevertheless, the $GW$ spectral function can provide a starting point for the cumulant expansion, in which the satellite and quasiparticle spectral function are convoluted and the position and weight of the satellite features are corrected.61, 62, 14, 63, 64, 65, 66, 54, 67

At present, there are two major limitations to obtain accurate spectral functions with the $GW$ approximation. First, the frequency dependence of the screened Coulomb interaction is commonly approximated by a Dirac delta function in the plasmon-pole model.12, 68, 69 To increase the accuracy, more recent approaches employ an analytic continuation70, 71, 72 or contour deformation.73, 74, 75 A detailed assessment of the accuracy of these techniques is still lacking and may prove crucial for reliable calculations of quasiparticle lifetimes.76, 77 Secondly, most $GW$ implementations rely on a summation over the empty Kohn-Sham states.12 This requires a very accurate description of the unoccupied subspace and therefore imposes additional constraints on the construction of pseudopotentials. Furthermore, converging the eigenvalues occasionally demands the inclusion of a considerable number of empty states78, 79, 80, and the convergence with respect to these states is interdependent with other convergence parameters.79, 81 To address these issues, several techniques have been proposed82, 83, 84, 85, 86, 87 which improve the rate of convergence of the summations. However, some of these methods introduce additional convergence parameters. Alternatively, one can obtain the dielectric matrix in density-functional perturbation theory88, 89 by solving linear Sternheimer equations.90 This approach has been subsequently optimized by iterative diagonalization to determine the most significant eigenvalues91, 92, 93 or employing an optimal representation of the polarization.94, 95 One can extend this approach, so that both the Green’s function $G$ and the screened Coulomb interaction $W$ can be obtained by solving linear response equations.96, 94 This approach is referred to as the Sternheimer $GW$ method and has been successfully applied to large systems with a single $\boldsymbol{k}$ point.94, 95, 97

In this work, we develop an implementation for solids building on the methodology described in Ref. 96 and 98. Our implementation, called SternheimerGW, features the plasmon pole model, the analytic continuation, as well as the real frequency integration on the same footing, and is therefore suitable to assess the accuracy of these techniques. The frequency dependent self-energy allows us to access spectral properties from first principles.77 Two illustrative examples of the application of SternheimerGW are shown in Fig. 1. In Fig. 1a we show the complete wavevector-dependent spectral function of silicon. For these results the frequency-dependent self-energy obtained with SternheimerGW is used as a starting point for a cumulant expansion, to obtain both quasiparticle band structure and plasmon replicas with the correct energy and intensity.40, 1 In Fig. 1b, we show the use of SternheimerGW to evaluate the Fermi surface of the transition metal dichalcogenide NbS2.2 In this system, the change of density of states at the Fermi energy due to many-body perturbation theory corrections is important to obtain more accurate predictions for the superconducting transition temperature.99

2 Description of the code

2.1 Kohn-Sham reference system

To evaluate the many-body perturbation correction, we start from an unperturbed reference system within Kohn-Sham density functional theory. The important ingredients of this reference system for our calculation are the exchange-correlation potential $V^{\text{xc}}$ and the Kohn-Sham Hamiltonian $\hat{H}_{\boldsymbol{k}}$ along with its eigenvalues $\varepsilon_{n\boldsymbol{k}}$ and wavefunctions $\phi_{n\boldsymbol{k}}(\boldsymbol{r})$. Here $n$ is a band index, $\boldsymbol{r}$ is the position vector, and $\boldsymbol{k}$ a Bloch vector in the Brillouin zone of the crystal. We employ a planewaves basis set, and express the wavefunction $\phi_{n\boldsymbol{k}}(\boldsymbol{r})$ in terms of the expansion coefficients $u_{n\boldsymbol{k},\boldsymbol{G}}$:

[TABLE]

where $\boldsymbol{G}$ denotes reciprocal lattice vectors and $\Omega$ is the volume of the unit cell. In a planewaves basis only the kinetic energy part of the Hamiltonian and the non-local part of the pseudopotentials are represented in reciprocal space; the local Kohn-Sham potential is applied in real space via a fast Fourier transform.

2.2 Green’s function

From the Kohn-Sham Hamiltonian, we obtain the electronic Green’s function, which is one of the ingredients to evaluate the $GW$ self-energy. The Green’s function is given by the linear equation:

[TABLE]

Here, $\omega$ is a complex frequency, chosen to satisfy the time ordering, $\mathcal{L}^{\text{G}}$ is the linear operator for the Green’s function:

[TABLE]

and $\delta_{\boldsymbol{G}\boldsymbol{G}^{\prime}}$ is the Kronecker $\delta$. Most $GW$ codes expand Eq. (2) in the basis of the eigenstates of the Hamiltonian and solve it by summing over many unoccupied states. In contrast, in SternheimerGW we utilize an iterative Krylov subspace solver100 to explicitly find the solution to the linear equation (2). We employ the fact that shifted linear problems described by the linear operator $\mathcal{L}^{\text{G}}$ differing only in the frequency $\omega$ span the same Krylov subspace. This property allows us to solve the linear problem for all frequencies at the computational complexity of a single frequency.

2.3 Dielectric matrix

The dielectric matrix of a system can be evaluated from the electron density response to a planewave perturbation:

[TABLE]

where $\omega$ is the frequency. The calculation of the density response is performed by starting from the linear variation of the wavefunctions in real space, and then carrying out a Fourier transform into reciprocal space:

[TABLE]

where the factor 2 accounts for the spin degeneracy and we are assuming a uniform grid of $N_{\boldsymbol{k}}$ $\boldsymbol{k}$-vectors in the Brillouin zone. The change of the wavefunction coefficients is obtained from the Sternheimer equation, which can be written as:

[TABLE]

Since the linear operator of the Coulomb response,

[TABLE]

is very similar to the one for the Green’s function (3), the same Krylov subspace solvers can be employed. The only difference with respect to (3) is that now we have an additional projector on the occupied states $\hat{P}^{\text{occ}}$; this term improves the condition number of the linear operator. The right-hand side of (6) is given by a Fourier transform of:

[TABLE]

and describes a driving field corresponding to a periodic perturbation.

2.4 Correlation self-energy

For the correlation self-energy, we combine the inverse of the dielectric matrix with the truncated Coulomb potential (see Sec. 2.6):

[TABLE]

where $\bar{W}$ is the correlation part of the screened Coulomb interaction. To overcome numerical issues associated with the inversion of the dielectric matrix for $q\approx 0$,101 we evaluate that element for a small but finite $\boldsymbol{q}$102 and set the off-diagonal elements to zero if $\boldsymbol{q}+\boldsymbol{G}=\boldsymbol{0}$ or $\boldsymbol{q}+\boldsymbol{G}^{\prime}=\boldsymbol{0}$.103 Subsequently, we perform a Fourier transform of $G$ and $\bar{W}$ to real space, because the self-energy can be expressed as a product there:

[TABLE]

$w_{\boldsymbol{q}}$ and $w_{\omega^{\prime}}$ are weights for $\boldsymbol{q}$ and $\omega$ integration respectively. The factor $\alpha$ depends on whether the frequency integration is performed along the real axis ($\alpha=\text{i}$) or along the imaginary axis ($\alpha=-1$). We transform the resulting self-energy back to reciprocal space to evaluate the expectation values with the Kohn-Sham wavefunctions.

2.5 Exchange self-energy

To evaluate the exchange self-energy, we evaluate the Fock potential by summing over all occupied wavefunctions:

[TABLE]

where the weight $w_{n\boldsymbol{k}-\boldsymbol{q}}$ contains the weight of the $\boldsymbol{q}$ integration. The Coulomb potential is truncated as described in the following section.

2.6 Truncation

Owing to the long-range nature of the Coulomb potential, the Fourier transform of this quantity diverges for small arguments $q\rightarrow 0$, leading to instabilities in the convergence with respect to the grid in the Brillouin zone. Furthermore, when describing systems of reduced dimensionality, such as slabs or nanowires, it is advantageous to truncate the Coulomb potential along the non-periodic direction in order to avoid spurious interactions between periodic images. In order to address these issues, it is common practice to truncate104, 105, 106 the Coulomb potential at a certain distance. This approach yields a Fourier transform with a finite value in the limit $q\rightarrow 0$. In SternheimerGW, we implement three different truncation schemes: First, one can truncate the Coulomb potential at a distance $R_{\text{cut}}$ (spherical truncation) resulting in a potential:

[TABLE]

Secondly, one can truncate the Coulomb interaction at a height $z_{\text{cut}}$ (slab truncation), which yields the following potential:104

[TABLE]

Lastly, one can truncate the Coulomb interaction in the Wigner-Seitz supercell.106 In this latter case, since an analytic expression for the potential is not known, we tabulate the results of the truncation and look up the values for the relevant vectors as necessary.

2.7 Analytic continuation

Quantities such as the dielectric matrix and the self-energy exhibit significant structure along the real axis, for example due to electron-hole excitations. In order to improve numerical stability and accuracy, it can be advantageous to evaluate quantities along the imaginary axis and perform an analytic continuation to the real axis. In SternheimerGW, there are two levels where such an analytic continuation may be employed. On the one hand, one can evaluate the dielectric function on the imaginary axis, perform an analytic continuation to the real axis, and subsequently convolute it with the Green’s function according to Eq. (10). On the other hand, one can resolve the convolution along the imaginary axis, thereby obtaining a self-energy at imaginary frequency values. The self-energy at real frequencies is then obtained via analytic continuation. To perform the analytic continuation, one can use a Padé expansion107 or the adaptive Antoulas-Anderson (AAA) algorithm.108 In general, determining the Padé approximant using a continued fraction expansion suffers from numerical instabilities, since a very high precision (number of significant digits) is required. In contrast, the AAA algorithm relies on a singular-value decomposition, and it refines iteratively the data points used to construct the approximant, until the largest deviation falls below a specified threshold.

2.8 Spectral function

The spectral function is related to the imaginary part of the retarded Green’s function:

[TABLE]

and represents a $\boldsymbol{k}$-resolved many-body density of states. Since the Green’s function is the resolvent of the many-body Hamiltonian, this leads to:

[TABLE]

Approximating the $GW$ wavefunctions by their Kohn-Sham counterpart $\phi_{n\boldsymbol{k}}(\boldsymbol{r})$, the diagonal part of the spectral function in the Kohn-Sham basis can be expressed as:

[TABLE]

The peaks of the spectral function correspond to the quasiparticle eigenvalues. Linearizing the equation in the vicinity of the eigenvalue $\varepsilon_{n\boldsymbol{k}}$ results in the following quasiparticle eigenvalue:

[TABLE]

where Z_{n\boldsymbol{k}}=1-\hbar^{-1}\text{Re}\bigl{(}\partial\Sigma_{n\boldsymbol{k}}^{\text{c}}/\partial\omega\bigr{)}.

2.9 Software design

Architecture

In Fig. 2, we illustrate the multitier architecture of SternheimerGW. The code builds on functionality provided by the Quantum Espresso software package109, 110 and provides additional functionality handling low-level operations in the data tier. On top of this layer, we have an algo tier handling the interface with Quantum Espresso. The algo tier also provides integration grids, linear solvers, Coulomb truncation, and analytic continuation functionalities. In the phys tier, we evaluate the $GW$ self-energy and matrix elements, employing the functionality of the lower lying tiers. The highest tier provides the user interface and the infrastructure for testing the software.

User input

The typical approach for specifying the user interface of a code requires that changes to input variables are manually translated into an update of the documentation. As this scheme poses the risk that the documentation becomes outdated, we employ an inverse scheme in SternheimerGW. The core file controlling the user input is a human- and machine-readable YAML file documenting all input variables (see example in Fig. 3a). Starting from the documentation, we run a script to generate the Fortran file that processes the input file (Fig. 3b). As a consequence, introducing a new input variable requires an update of the documentation. This strategy allows us to make sure that the code and the documentation are always in sync. Furthermore, we process the YAML file to generate the online documentation of SternheimerGW (Fig. 3c) on the website http://sternheimergw.org.

Continuous integration and testing

To ensure the long-term stability of the code and its alignment with the Quantum Espresso suite, we employ two techniques. First, we examine whether the code compiles and reproduces a set of reference results using a buildbot test farm.111, 112 These integration tests involve 12 different calculations to check most of the functionality of SternheimerGW. Secondly, we perform unit testing of some individual packages of the code, by benchmarking the parallelization, the algebra, the analytic continuation, the linear solver, and the input/output (io) package independently from the rest of the code. For the unit testing, at variance with the tests of the complete code, we compare our implementation with analytic results. The unit tests are implemented in the pFUnit framework.113 Overall, all these tests cover at least 80% of the code in each tier (see Table 1).

2.10 Parallelization

$GW$ calculations are computationally much more demanding than Kohn-Sham DFT calculations. To allow for a reasonable time to solution, we employ two different parallelization levels in SternheimerGW. As we inherit the distribution logic from Quantum Espresso, we refer to these levels as pool and image parallelization. In Fig. 4, we compare the performance of these parallelization levels for a bulk system of MgO for different parts of the code. With the pool parallelization, we distribute different $\boldsymbol{k}$ or $\boldsymbol{q}$ points to different processors. This method is very efficient for both $G$ and $W$, if the number of these points is a multiple of the number of CPUs. The image parallelization distributes $\boldsymbol{G}$ vectors used in the linear response and Fourier transform to different CPU. For the calculation of $W$, this parallelization is very efficient because the solution of the linear-response equation for each $\boldsymbol{G}$ vector is independent of the other $\boldsymbol{G}$ vectors. When evaluating $\Sigma^{\text{c}}$ the speedup of this parallelization saturates when the number of grid points along the $z$ axis reduces to one per CPU. The benefit of the image parallelization is that it reduces the memory footprint of the code.

3 Examples

3.1 Benchmark results

To assess the accuracy of SternheimerGW, we verify the implementation by comparing to reference calculations from the literature.114, 115, 116, 117 We use seven reference systems in the diamond (C, Si, and Ge), the zincblende (BN, CdS, and GaAs), or rocksalt structure (MgO). All lattice constants are taken from Ref. 114, except for that of Ge which is from Ref. 120.

The DFT ground state calculation was run with Quantum ESPRESSO v6.3 using a $12\times 12\times 12$ $\boldsymbol{k}$-point mesh and the PBE exchange-correlation potential.121 The pseudopotentials are verified according to the $\Delta$-test;122 as a compromise between speed and accuracy, we employ the ONCV123 pseudopotentials from the SG15 library118 for compounds not including $d$ electrons, and the stringent pseudopotentials v0.4 from pseudo-dojo.org119 otherwise.

In Table 2, we present the comparison between the reference calculations and the results obtained with SternheimerGW v0.15. The SternheimerGW calculations were converged using the algorithm described in Sec. 3.2 resulting in the parameters listed in Table 3. Overall, we find that our results agree well with the literature, within an error bar of 0.1 eV usually quoted in the literature for $GW$ calculations.

3.2 Converging a GW calculation

Predictive $GW$ calculations require the convergence of several numerical parameters. An efficient strategy is required to obtain these values at a reasonable computational cost. In order to tackle this multi-dimensional optimization, a possible strategy is to keep all convergence parameters at a loose setting, and only focus on the convergence behavior of a single parameter. To determine reasonable values for the loose setting it is a good strategy to focus on parts of the calculation that can be evaluated separately from the rest. First, one considers the exchange contribution to the $GW$ self-energy, which is the most challenging quantity to converge in SternheimerGW, but at the same time its calculation is significantly less demanding than the correlation part. The analysis of the exchange self-energy yields estimates for the wavefunction cutoff $E_{\text{cut}}^{\text{scf}}$ used in the DFT calculation (Fig. 5a), the exchange cutoff $E_{\text{cut}}^{\text{x}}$ (Fig. 5b), and the number $N_{\boldsymbol{q}}$ of points in the $\boldsymbol{q}$-point grid used to evaluate the Brillouin-zone integrals (Fig. 5c). The wavefunction and the exchange cutoff are somewhat intertwined, and occasionally the former needs to be increased to allow for convergence of the latter. Secondly, we investigate the convergence of the head of the dielectric function ($\boldsymbol{q}=0$ point) with respect to the number of $\boldsymbol{k}$ points $N_{\boldsymbol{k}}$ (Fig. 5d) and the coarse frequency mesh $N_{\omega,W}^{\text{coarse}}$ used for the analytic continuation (Fig. 5e).

Guided by the initial studies of the exchange self-energy and the dielectric matrix, we construct the loose setting for the convergence study. In general, we will choose the smallest possible values in the smooth region of convergence. For the correlation, we need to set three further numerical parameters: For the dense frequency mesh $N_{\omega,W}^{\text{dense}}$ (Fig. 5f), we choose a mesh twice as dense as the coarse mesh; for the correlation cutoff $E_{\text{cut}}^{\text{c}}$ (Fig. 5g), we select typically a third of the value used for the exchange self-energy; and we use just two frequencies $N_{\omega,\Sigma}^{\text{coarse}}=2$ (Fig. 5h) to perform the analytic continuation of the correlation self-energy from the imaginary to the real frequency axis.

With the constructed loose mesh, we then perform the convergence study individually increasing every single one of the parameters until the quasiparticle energies of interest are converged. The advantage of this strategy is that the computational cost of the complete convergence study is significantly smaller than that of the subsequent $GW$ calculation. In Fig. 5, we demonstrate that the estimate based on the loose setting provides a conservative estimate of the stringent parameters needed to obtain converged results.

3.3 Spectral function

One noteworthy feature of SternheimerGW is the possibility of calculating the $\boldsymbol{k}$-resolved spectral function at any arbitrary point in the Brillouin zone. This feature allows us to compute $GW$ band structures without the need for interpolation techniques. The frequency resolution provides insight into the electronic lifetimes, and can be directly compared to ARPES experiments. In Fig. 6, we show a representative calculation for silicon. We can see the quasiparticle bands and the weaker replica of the bands due to the plasmon satellite. We note that the energy and broadening of the plasmon satellites are not accurately described in $GW$. An accurate description of these features requires the inclusion of higher order effects via the cumulant expansion.61, 62, 14, 63, 64, 65, 66, 54, 67

3.4 Real frequency integration

As outlined in Sec. 2.4, SternheimerGW provides the option to perform the frequency convolution along the real or the imaginary frequency axis. In Fig. 7, we demonstrate that these two implementations give nearly identical results for the spectral function at high-symmetry points of silicon. The main identifiable differences is a small shift of the valence band top to lower energies and an increased noise when using integration along the real axis. Both of these differences originate from the small parameter $\eta=0.15$ eV that is used to shift the frequency slightly off the real axis. This shift is necessary in order to avoid the singularities on the real axis, so as to obtain numerically stable results. On the one hand, if $\eta$ is decreased to 0.08 eV the resulting spectral function becomes very noisy (see Fig. 8a). On the other hand, if $\eta$ is increased to 0.3 eV, the energy of the quasiparticle becomes inaccurate (see inset of Fig. 8a).

In general, the integration along the imaginary axis with subsequent analytic continuation is one to two orders of magnitude cheaper because the Green’s function and the screened Coulomb interaction have less structure there. As the results of both methods are very similar, we chose the integration along the imaginary axis as the default in SternheimerGW, and we provide the integration along the real axis as an optional feature for validation.

In Fig. 8b, we show that obtaining accurate quasi-particle lifetimes requires going beyond the plasmon-pole approximation, because in this approximation there is only one sharp pole at the plasmon frequency (cf. Fig.9a). In contrast, analytic continuation techniques can produce multiple poles and thereby capture the broadening of the quasiparticle spectra. In Fig. 9b, we show that the Padé expansion107 can reproduce the shape of the dielectric function more accurately than the plasmon-pole approximation, since it naturally incorporates multiple poles. However, this Padé algorithm is numerically not stable when using a large number of frequencies, therefore it cannot be used for the convolution along the real frequency axis. To address this issue, we employ the AAA analytic108 continuation technique, whereby only significant frequencies are selected via a singular value decomposition. Subsequently, poles with weak residues are removed so that only the most important contributions remain. In Fig. 9c, we show that the AAA method reduces the number of poles compared to the Padé expansion while still reproducing the overall shape and magnitude of the dielectric matrix.

4 Conclusion

We presented the SternheimerGW code, which calculates many-body $GW$ self-energies, quasiparticle band structures, and spectral functions by solving linear-response Sternheimer equations. The linear response scheme allows us to calculate the $GW$ self-energy at arbitrary $\boldsymbol{k}$ points. This is particularly useful for computing band structures without resorting to interpolation techniques, and for comparison to ARPES experiments. SternheimerGW goes beyond the plasmon-pole approximation by performing full-frequency integration along the real or the imaginary axis.

SternheimerGW is aligned with the version 6.3 of Quantum Espresso. The dedicated website http://www.sternheimergw.org provides users with installation guidelines, tutorials, and a comprehensive documentation of all input variables. The code is efficiently parallelized and has been tested extensively on high-performance computing architectures. Updated versions of the code are distributed in the GitHub repository https://github.com/QEF/SternheimerGW.

Acknowledgments

The development of SternheimerGW has received funding from the Leverhulme Trust (Grant RL-2012-001), the Graphene Flagship (Horizon 2020 Grant No. 785219 - GrapheneCore2), and the UK Engineering and Physical Sciences Research Council (Grant No. EP/M020517/1). Furthermore, the authors acknowledge the use of the University of Oxford Advanced Research Computing (ARC) facility (http://dx.doi.org/10.5281/zenodo.22558), the ARCHER UK National Supercomputing Service under the “T-Dops” project, the Cambridge Service for Data Driven Discovery (CSD3) funded by EPSRC Tier-2 (Grant EP/P020259/1), the DECI resource “Cartesius” based in the Netherlands at SURFsara and “Abel” based in Oslo with support from the PRACE AISBL, and PRACE for awarding us access to MareNostrum at BSC-CNS, Spain.