Ab initio Sternheimer-GW method for quasiparticle calculations using plane waves

We report on the extension and implementation of the Sternheimer-GW method introduced by Giustino et al. [Phys. Rev. B 81, 115105 (2010)] to the case of first-principles pseudopotential calculations based on a plane-waves basis. The Sternheimer-GW method consists of calculating the GW self-energy operator without resorting to the standard expansion over unoccupied Kohn-Sham electronic states. The Green's function is calculated by solving linear systems for frequencies along the real axis. The screened Coulomb interaction is calculated for frequencies along the imaginary axis by using the Sternheimer equation. Analytic continuation to the real axis is performed using Padé approximants. The generalized plasmon-pole approximation is avoided by performing explicit calculations at multiple frequencies using Frommer's multishift solver. We demonstrate our methodology by reporting tests on common insulators and semiconductors, including Si, diamond, LiCl, and SiC. Our calculated quasiparticle energies are in agreement with the results of fully converged calculations based on the sum-over-states approach. As the Sternheimer-GW method yields the complete self-energy $\Sigma (\mathbf{r},{\mathbf{r}}^{\prime },\omega )$ and not only its expectation values on Kohn-Sham states, this work opens the way to nonperturbative GW calculations and to direct calculations of spectral functions for angle-resolved photoemission spectroscopy. As an example of the capabilities of the method we calculate the G${}_0$W${}_0$ spectral functions of silicon and diamond.

Figure 1

Parallel execution time for a Sternheimer-GW calculation of silicon. The time refers to a calculation of the complete self-energy ${\Sigma }_{\mathbf{k}}(\mathbf{G},{\mathbf{G}}^{\prime },\omega )$ with the parameters given in Sec. 3c. We report the total execution time (full black bars), the time required for calculating the screened Coulomb interaction (straight cross-hatched blue bars), the Green's function (oblique cross-hatched red bars), and the frequency convolution in Eq. (18) (oblique hatched black bars).

Figure 2

Quasiparticle corrections to the band edge states at high-symmetry points in silicon: Sternheimer-GW (solid blue line) and sum-over-states approaches as implemented in sax (red disks and dashed line). The corrections are shown as a function of the energy of the highest unoccupied state included in the sum-over-states calculation. The zero of the energy is set to the top of the valence band. (a) Band edges at $\Gamma$, (b) band edges at $X$, and (c) band edges at $L$.

Figure 3

Quasiparticle corrections to the band edge states at high-symmetry points in diamond: Sternheimer-GW (solid blue line) and sum-over-states approaches as implemented in sax (red disks and dashed line). The corrections are shown as a function of the energy of the highest unoccupied state included in the sum-over-states calculation. The zero of the energy is set to the top of the valence band. (a) Band edges at $\Gamma$, (b) band edges at $X$, and (c) band edges at $L$.

Figure 4

(a) Quasiparticle spectral function $A_{\mathbf{k}}\left(\omega \right)$ of silicon for $\mathbf{k}$ along $\Gamma$-$X$, calculated using the Sternheimer-GW method within the diagonal ${\text{G}}_0{\text{W}}_0$ approximation [Eq. (19)]. The “discrete” structure visible in the more dispersive band is simply a visualization artifact, resulting from our choice of computing the self-energy at 20 equally spaced $\mathbf{k}$ points. (b) DFT/LDA band structure of silicon along $\Gamma$-$X$ (black dashed lines), and the corresponding quasiparticle band structure obtained from (a) (blue solid lines). The units of the color bar are eV${}^{-1}$.

Figure 5

(a) Quasiparticle spectral function $A_{\mathbf{k}}\left(\omega \right)$ of diamond for $\mathbf{k}$ along $\Gamma$-$L$, calculated using the Sternheimer-GW method within the diagonal ${\text{G}}_0{\text{W}}_0$ approximation [Eq. (19)]. The “discrete” structure visible in the more dispersive band is simply a visualization artifact, resulting from our choice of computing the self-energy at 20 equally spaced $\mathbf{k}$ points. (b) DFT/LDA band structure of diamond along $\Gamma$-$L$ (black dashed lines), and the corresponding quasiparticle band structure obtained from (a) (blue solid lines). The units of the color bar are eV${}^{-1}$.

Figure 6

(a) Quasiparticle spectral function $A_{\mathbf{k}}\left(\omega \right)$ of silicon for $\mathbf{k}$ along $\Gamma$-$X$, calculated using the Sternheimer-GW method within the diagonal ${\text{G}}_0{\text{W}}_0$ approximation [Eq. (19)]. The energy range extends to $-40$ eV in order to show the plasmaron band structure. The spectral function is given in a logarithmic scale in order to enhance the plasmaron satellites. (b) Zoom on the plot in (a) around the binding energy 15–40 eV, showing the plasmaron band structure. In this case we use a linear scale with values given by the color bar. We note that the satellite energy is incorrect in the ${\text{G}}_0{\text{W}}_0$ approximation (Refs. 13 and 84, 85, 86). The units of the color bar are eV${}^{-1}$.