Simple eigenvalue-self-consistent $\bar{\Delta}GW_{0}$

TL;DR

The paper introduces a simplified eigenvalue self-consistent $GW$ method, $ar{ riangle}GW_{0}$, which improves accuracy with system size and is computationally efficient for large-scale systems like semiconductors and insulators.

Contribution

It derives a new simplified eigenvalue self-consistency scheme for $GW$, called $ar{ riangle}GW_{0}$, applicable to large systems with improved accuracy.

Findings

Accuracy increases with system size, approaching CCSD(T) results.

Method is computationally comparable to one-shot $G_{0}W_{0}$.

Effective for large periodic systems with minimal errors.

Abstract

We derive a general form of eigenvalue self-consistency for $GW_{0}$ in the time domain and use it to obtain a simplified postprocessing eigenvalue self-consistency, which we label $\bar{\Delta}GW_{0}$. The method costs the same as a one-shot $G_{0}W_{0}$ when the latter gives the full frequency-domain (or time-domain) matrix element of the self-energy. The accuracy of $\bar{\Delta}GW_{0}$ increases with system size, as demonstrated here by comparison to other $GW$ self-consistency results and to CCSD(T) predictions. When combined with the large-scale stochastic $G_{0}W_{0}$ formulation $\bar{\Delta}GW_{0}$ is applicable to very large systems, as exemplified by periodic supercells of semiconductors and insulators with 2048 valence electrons. For molecules the error of our eventual partially self-consistent approach starts at about 0.2eV for small molecules and decreases to 0.05eV for large ones, while for the periodic solids studied here the mean-absolute-error is only 0.03eV.

Tables

Table 1. Table 1: Ionization potentials I 𝐼 I for small molecules as calculated by different levels of theory. The G 0 ​ W 0 subscript 𝐺 0 subscript 𝑊 0 G_{0}W_{0} and Δ ¯ ​ G ​ W 0 ¯ Δ 𝐺 subscript 𝑊 0 \bar{\Delta}GW_{0} estimates were obtained using an LDA starting point and calculated using a stochastic implementation Neuhauser  et al. ( 2014 ); Vlcek  et al. ( 2017 ) with the statistical errors reported in parentheses. Experimental values are from Ref.   NIS, . For several of the acenes the CCSD(T) results are estimates extrapolated to the infinite basis-set limit. Rangel  et al. ( 2016 ) Molecular geometries were taken from Refs.  NIS, and Rangel  et al. , 2016 .

system	$h$ ($a_0$)	$I\left(\mathrm{eV}\right)$
		LDA	$G_{0}W$		$\overline{\Delta }G{W}_0$		ev$G{W}_0$	qp$GW$	CCSD(T)	Exp.
nitrogen	0.35	10.44	15.08	(0.05)	15.93	(0.05)	15.32111Ref. Caruso et al.,2016	16.01a	15.57222Ref. Krause et al.,2015	15.58
ethylene	0.35	6.92	10.50	(0.04)	10.87	(0.04)	10.24a	10.63a	10.67b	10.68
urea	0.30	6.10	9.53	(0.08)	10.48	(0.08)	9.81a	10.45a	10.05b	10.28
naphtalene	0.35	5.71	8.10	(0.09)	8.39	(0.09)	8.15333Ref. Rangel et al.,2016	-	8.25c	8.14
tetracene	0.35	4.89	6.79	(0.08)	6.94	(0.08)	6.84c	-	7.02c	6.97
hexacene	0.35	4.52	6.15	(0.06)	6.33	(0.06)	6.19c	-	6.32c	6.33

Table 2. Table 2: Fundamental band gaps ( E g subscript 𝐸 𝑔 E_{g} ) for a sample of solids. For each system a 4 × 4 × 4 4 4 4 4\times 4\times 4 conventional super-cell (with 2048 valence electrons) was used. Values labeled by * are for the minimum fundamental gap, which is not accessed in our calculations (see text).

system	$h$ ($a_0$)	$E_g$ (eV)
		LDA	$G_0{W}_0$		$\overline{\Delta }G{W}_0$		ev$GW$ (Ref. Shishkin and Kresse, 2007b)	qp$GW$ (Ref. Shishkin and Kresse, 2007b)	Exp.
Si	0.446	0.56	1.29	(0.04)	1.35	(0.04)	1.20*	1.28*	1.3444Ref. Tiago et al., 2004 (1.17*)a
SiC	0.293	1.37	2.29	(0.04)	2.35	(0.04)	2.43	2.64	2.42 555Ref. Bimberg et al.,1981
AlP	0.368	1.46	2.41	(0.03)	2.50	(0.03)	2.59	2.77	2.52 c
C	0.336	4.16	5.40	(0.06)	5.47	(0.06)	5.50	5.99	5.48 d
BN	0.380	4.48	6.21	(0.06)	6.41	(0.07)	6.10	6.73	6.1 - 6.4 e

Equations

$$G\left(\mathbf{r},{\mathbf{r}}^{\prime },t,t^{\prime }\right)=G_0\left(\mathbf{r},{\mathbf{r}}^{\prime },t,t^{\prime }\right)+\iint G_0\left(\mathbf{r},{\mathbf{r}}_1,t,t_1\right)\times$$

$$\Sigma \left({\mathbf{r}}_1,{\mathbf{r}}_2,t_1-t_2\right)G\left({\mathbf{r}}_2,{\mathbf{r}}^{\prime },t_2,t^{\prime }\right)d{\mathbf{r}}_{1}d{t}_{1}d{\mathbf{r}}_{2}d{t}_2,$$

$$\Sigma \left(\mathbf{r},{\mathbf{r}}^{\prime },t\right)=iG\left(\mathbf{r},{\mathbf{r}}^{\prime },t\right)W\left(\mathbf{r},{\mathbf{r}}^{\prime },t^{+}\right),$$

$${\epsilon }^{QP}\left(i\right)={\epsilon }^0\left(i\right)+\Delta \left(i\right),$$

$$\Delta \left(i\right)={\tilde{\Sigma }}_i\left[\omega ={\epsilon }^{QP}\left(i\right)\right]-⟨{\varphi }_i\left|v_{xc}\right|{\varphi }_i⟩.$$

$${\tilde{\Sigma }}_i\left[\omega \right]=\int ⟨{\varphi }_i\left|\hat{\Sigma }\left(t\right)\right|{\varphi }_i⟩e^{i\omega t}dt,$$

$$i{G}_0\left(\mathbf{r},{\mathbf{r}}^{\prime },t\right)$$

$$\left[\theta \left(t\right){\theta }_{\beta }\left({\hat{h}}^0-\mu \right)-\theta \left(-t\right){\theta }_{\beta }\left(\mu -{\hat{h}}^0\right)\right]\},$$

$${\hat{h}}^0=-\frac{1}{2}{\nabla }^2+v_{\mathrm{ext}}+v_{\mathrm{H}}+v_{\text{xc}},$$

$$⟨{\varphi }_i\left|{\hat{\Sigma }}_0\left(t\right)\right|{\varphi }_i⟩=$$

$$i\iint {\varphi }_i\left(\mathbf{r}\right)G_0\left(\mathbf{r},{\mathbf{r}}^{\prime },t\right)W_0\left(\mathbf{r},{\mathbf{r}}^{\prime },t^{+}\right){\varphi }_i\left({\mathbf{r}}^{\prime }\right)d\mathbf{r}d{\mathbf{r}}^{\prime }.$$

$$iG\left(\mathbf{r},{\mathbf{r}}^{\prime },t\right)$$

$$\left[\theta \left(t\right){\theta }_{\beta }\left({\hat{h}}^0-\mu \right)-\theta \left(-t\right){\theta }_{\beta }\left(\mu -{\hat{h}}^0\right)\right]\},$$

$$\hat{\Delta }\simeq \sum _i\left|{\varphi }_i⟩\mathrm{Re}\Delta \left(i\right)⟨{\varphi }_i\right|.$$

$$iG\left(\mathbf{r},{\mathbf{r}}^{\prime },t\right)$$

$$\left[\theta \left(t\right){\theta }_{\beta }\left({\hat{h}}^0-\mu \right)-\theta \left(-t\right){\theta }_{\beta }\left(\mu -{\hat{h}}^0\right)\right]\}.$$

$$\mathrm{Re}\bar{\Delta }\left[{\epsilon }^0\left(i\right)\right]\approx \{\begin{array}{ll}\Delta \left(H\right)& {\epsilon }^0\left(i\right)\le {\epsilon }^0\left(H\right)\\ \Delta \left(L\right)& {\epsilon }^0\left(i\right)\ge {\epsilon }^0\left(L\right)\end{array}.$$

$$\bar{\Delta }\left(t\right)\equiv \{\begin{array}{ll}\Delta \left(H\right)& t<0\\ \Delta \left(L\right)& t>0\end{array}.$$

$$⟨{\varphi }_i\left|\hat{\Sigma }\left(t\right)\right|{\varphi }_i⟩=e^{-i\bar{\Delta }\left(t\right)t}⟨{\varphi }_i\left|{\hat{\Sigma }}_0\left(t\right)\right|{\varphi }_i⟩.$$

$$E_g={\epsilon }_L-{\epsilon }_H.$$

$$\bar{\Delta }\left(\epsilon \right)=a_2{\left(\epsilon \right)}^2+a_1\epsilon +a_0.$$

Full text

I Introduction

The $GW$ approximation Hedin (1965) to many-body perturbation theory is often used to calculate electron removal or addition energies and related (inverse) photoemission spectra of molecules, nanostructures, and bulk materials.Aryasetiawan and Gunnarsson (1998); Rieger et al. (1999); Steinbeck et al. (1999); Onida et al. (2002); Rubio and Louie (2005); Friedrich and Schindlmayr (2006); Shishkin and Kresse (2007a); Trevisanutto et al. (2008); Rostgaard et al. (2010a); Tamblyn et al. (2011); van Setten et al. (2012); Stefanucci and van Leeuwen (2013); Govoni and Galli (2015); Kaplan et al. (2016a) $GW$ is part of a family of methods that describe the probability amplitude of a quasiparticle (QP) to propagate between two space-time points $\left({\bf r},t\right)$ and $\left({\bf r}^{\prime},t^{\prime}\right)$ with a Green’s function $G\left({\bf r},{\bf r}^{\prime},t,t^{\prime}\right)$, the poles of which are the QP energies. The Green’s function is obtained perturbatively from a reference (non-interacting) Green’s function, $G_{0}\left({\bf r},{\bf r}^{\prime},t,t^{\prime}\right)$, via a Dyson equation (all equations use atomic units):

[TABLE]

where $\Sigma\left({\bf r}_{1},{\bf r}_{2},t,t^{\prime}\right)$ represents the self-energy. The reference Green’s function is typicallyHybertsen and Louie (1985, 1986a) given by the Kohn-Sham Kohn and Sham (1965) (KS) density function theory (DFT).Hohenberg and Kohn (1964) In $GW$ , the self-energy $\Sigma$ is approximated as:

[TABLE]

where $W\left({\bf r},{\bf r}^{\prime},t\right)$ is the screened Coulomb interaction, usually evaluated within the random phase approximation (RPA).

A solution of the equations above in requires, in principle, a self-consistent procedure since both $\Sigma\left({\bf r},{\bf r}^{\prime},t\right)$ and $W\left({\bf r},{\bf r}^{\prime},t^{+}\right)$ depend on $G\left({\bf r},{\bf r}^{\prime},t\right)$. In practice, the self-consistency is often abandoned and the most common $GW$ treatment is based on “one-shot” schemeHybertsen and Louie (1985, 1986a), which we label as $G_{0}W_{0}$, since the right hand side of Eq. (2) becomes $G_{0}\left({\bf r},{\bf r}^{\prime},t\right)W_{0}\left({\bf r},{\bf r}^{\prime},t\right)$, where $W_{0}$ is obtained from a random phase approximation that uses the KS eigenstates. The one-shot approach improves significantly the KS-DFT results, yet it depends on the choice of the reference system and often underestimates the QP gaps ($E_{g}$) and the ionization potentials ($I$).Faleev et al. (2004); Shishkin and Kresse (2007a); Caruso et al. (2012a); Bruneval and Marques (2012); Bruneval (2012); van Setten et al. (2015a); Knight et al. (2016) A fully self-consistent solution is computationally extremely demandingStan et al. (2006, 2009); Rostgaard et al. (2010b); Blase and Attaccalite (2011); Deslippe et al. (2012); Nguyen et al. (2012); Caruso et al. (2012b, 2013); Koval et al. (2014); Wang (2015) and in many situations it yields results that are worse than $G_{0}W_{0}$.Holm and von Barth (1998); Shishkin and Kresse (2007a); Bruneval and Gatti (2014)

To simplify the problem, a static and Hermitian approximation to the self-energy Faleev et al. (2004); van Schilfgaarde et al. (2006) is sometimes used in the so-called QP self-consistent $GW$ (qp$GW$), which iteratively updates $\Sigma$ and QP wave-functions (i.e., Dyson orbitals). The qp$GW$ approach is still computationally expensive and cannot be applied for large systems; further it tends to overestimate $E_{g}$ van Schilfgaarde et al. (2006); Kotani et al. (2007); Bruneval and Gatti (2014) and the ionization potentials Kaplan et al. (2015, 2016b); Caruso et al. (2016), due to an overly strong screened Coulomb interaction. van Schilfgaarde et al. (2006); Bruneval and Gatti (2014); Kaplan et al. (2016b) Alternatively, the $W$ term is kept “frozen” and self-consistency is sought only in the Green’s function.Shishkin and Kresse (2007a, b); Stan et al. (2009); Blase and Attaccalite (2011); Bruneval and Gatti (2014); Kaplan et al. (2015) This method is termed eigenvalue self-consistent $GW$ (ev$GW_{0}$) and it was applied successfully to bulk systemsNorthrup et al. (1987); Shishkin and Kresse (2007a); Bruneval and Gatti (2014) and to organic molecules,Stan et al. (2009); Blase and Attaccalite (2011); Knight et al. (2016) with remarkable success. Even though it is cheaper than other self-consistency methods, $\Sigma$ has to be recalculated in each iteration, making ev$GW_{0}$ out of reach for nanoscale systems with thousands of occupied electronic states.

Here a time domain formulation (Sec. II) is used to derive a simplified ev$GW_{0}$ formalism, labeled $\bar{\Delta}GW_{0}$, where the self-consistency is only a postprocessing step. Hence, as long as one has access to the matrix element of the self-energy at all frequencies or all times, then, irrespective of system size, the computational cost of the self-consistency is negligible (i.e., seconds on a single-core machine) so $\bar{\Delta}GW_{0}$ costs not more than $G_{0}W_{0}.$

We specifically combine $\bar{\Delta}GW$ with our stochastic $G_{0}W_{0}$ approach, which has a nearly linear scalingNeuhauser et al. (2014); Vlcek et al. (2017) and enables $G_{0}W_{0}$ for extremely large systemsVlček et al. (2016, 2018). The stochastic $G_{0}W_{0}$ method has automatically the necessary ingredient for $\bar{\Delta}GW_{0}$, as it produces the matrix element of the self-energy at all times.

The combined method (stochastic $G_{0}W_{0}$ with $\bar{\Delta}GW_{0}$) is first tested in Sec. III on molecules, and we find that $\bar{\Delta}GW$ becomes more accurate as the system size increases. Next, we perform stochastic $G_{0}W_{0}$ calculations for periodic semiconductors and insulators using large supercells with 2048 valence electrons. For solids $\bar{\Delta}GW_{0}$ gives $E_{g}$ in excellent agreement with experiment and a mean absolute error of 0.03 eV.

In all cases the self-consistency is reached in very few iterations without any additional cost on top of the $G_{0}W_{0}$ step.

II Theory

II.1 Green’s function self-consistency in the time domain

The QP energy of the $i^{{\rm th}}$ state $\varepsilon^{QP}\left(i\right)$ is calculated using the usual form of the perturbative $GW$ approximation in which the Kohn-Sham eigenvalues ($\varepsilon^{0}$) are corrected by the QP shift ($\Delta$) using a fixed point equation:

[TABLE]

where

[TABLE]

Here, $v_{xc}$ is the Kohn-Sham exchange-correlation potential for the DFT density, $\tilde{\Sigma}_{i}\left[\omega\right]$ is the Fourier transform of the matrix-element of $\hat{\Sigma}(t)$

[TABLE]

and $\hat{\Sigma}(t)$ is given by Eq. (2).

Starting from a KS DFT reference point, the initial self-energy is constructed from the KS propagator

[TABLE]

where ${\rm Tr}$ denotes a trace over all KS states, $\mu$ is the chemical potential, $\theta$ is the Heaviside step function that guarantees forward and backward time propagation for particles and holes, respectively, and $\hat{h}^{0}$ is the KS Hamiltonian

[TABLE]

where we introduced the kinetic energy and the external and Hartree potentials. In the rest of the paper we employ real time-dependent Hartree propagation to calculate the screened Coulomb interaction Neuhauser et al. (2014); Vlcek et al. (2017); this is equivalent to using the RPA approximation for $W_{0}$.

In the time-domain, the self-energy matrix element for the $i^{{\rm th}}$ state is

[TABLE]

Finally, after Fourier transformation combined with time-ordering Neuhauser et al. (2014); Vlcek et al. (2017) the “one-shot” QP energy is calculated through Eq. (4).

In the ev$GW_{0}$ procedure, the Green’s function is reconstructed in each iteration, employing the QP energies from the previous iteration. In the time domain this corresponds to writing the propagator as

[TABLE]

where $\hat{\Delta}$ contains all the many-body contributions.

As common in ev$GW_{0}$ self-consistency Shishkin and Kresse (2007b); Bruneval and Gatti (2014); Kaplan et al. (2015), the fact that the true self-energy operator is non-Hermitian and non-diagonal is disregarded and, for the purposes of Eq. (9), it is expressed in the KS basis as

[TABLE]

Hence, all the KS energies in the exponent in Eq. (9) are shifted to the QP energies obtained from Eq. (3). Eq. (10) is basically the fundamental eqaution of ev$GW_{0}$.

We now note that the operator $\hat{\Delta}$ in Eq. (9) is by construction diagonal in the KS basis set and we thus express it as a function of the KS Hamiltonian $\bar{\Delta}\left(\hat{h}^{0}\right)$ that interpolates all QP shifts. The Green’s function is therefore:

[TABLE]

This simple expression allows for a further approximation described below that significantly reduces the computational cost associated with self-consistent treatment.

II.2 Efficient and inexpensive implementation

In many cases $\bar{\Delta}$ is well described by a low degree polynomial with discontinuity at the band gap energies, $\varepsilon\left(H\right)$ and $\varepsilon\left(L\right),$ corresponding to the highest occupied ($H$) and lowest unoccupied ($L$) states, respectively. The zeroth order term in this polynomial corresponds to a scissors operatorFilip and Giustino (2014); Qian et al. (2015) which shifts the occupied and unoccupied states down and up in energy, respectively

[TABLE]

We call this approximation $\bar{\Delta}GW_{0}$ and use it in Sec. III for molecules and periodic systems.

Combining Eqs. (12) and (11) leads to a modified Green’s function which acquires an additional phase shift that is different for positive and negative times. Therefore, in the time domain we can define:

[TABLE]

In each iteration, the updated self-energy matrix element is then calculated as

[TABLE]

Next, the self-energy matrix element is transformed to the frequency domain and used in Eq. (3) to calculate a new estimate of the QP energy. The new QP energy is used iteratively to update Eqs. (12) and (13). The full cycle is illustrated in Fig. 1.

Note that this form of self-consistency is trivial and is a postprocessing step with no additional cost, unlike previous usesFilip and Giustino (2014); Qian et al. (2015) of the scissors-operator in $GW$ which require repeated evaluations of the self-energy. Further, the $\bar{\Delta}GW_{0}$ approach is applicable to any implementation which yields $\Sigma\left(\omega\right)$. It is naturally suited for the stochastic $G_{0}W_{0}$ method Neuhauser et al. (2014); Vlcek et al. (2017) which provides the self-energy on the full-time domain and therefore on a wide range of frequencies (spanning several hundred eV).

III Results and Discussion

III.1 Molecules

We first test our approach on ionization potentials $I$ (taken as $-\varepsilon_{H}$) for a set of small molecules listed in Table. 1. A ground state DFT calculation is performed using a Fourier real-space grid, ensuring (using the Martyna-Tuckerman approach)Martyna and Tuckerman (1999) that the potentials are not periodic. The exchange-correlation interaction is described by local density approximation (LDA) Perdew and Wang (1992) with Troullier-Martins pseudopotentialsTroullier and Martins (1991); the DFT eigenvalues are converged up to$<$10 meV with respect to the spacings of the real space grids (given in Table 1).

The systems listed in the table are ordered according to the number of valence electrons; N2 and hexacene are the smallest and the largest molecules studied here. In all cases, the stochastic $G_{0}W_{0}$ approachNeuhauser et al. (2014); Vlcek et al. (2017) was used to calculate the self-energy. We compare our calculations with reference values taken from experiment and from CCSD(T). The geometries of the acene molecules are taken from the $GW$ and CCSD(T) benchmark in Ref. Rangel et al., 2016.

Compared to the LDA eigenvalues, one-shot $G_{0}W$ predictions for the ionization potentials are much closer to the CCSD(T) values with a mean absolute error of 0.29 eV. In all cases, the value of $I$ is underestimated, in agreement with previous benchmark studiesvan Setten et al. (2015b); Vlcek et al. (2017); Rangel et al. (2016). As the system size increases the difference between the $G_{0}W$ and CCSD(T) values decreases so that the one-shot correction is an increasingly better approximation.

The simplified eigenvalue self-consistency converges in 3-4 iterations after the initial $G_{0}W$ calculation; the initial and final self-energy curves are illustrated for hexacene in Fig. 2. Since our self-consistency procedure is merely a postprocessing step, its computational cost is negligible (less than a second on single core machine).

We first compare our $\bar{\Delta}GW$ results with previous eigenvalue and quasiparticle self-consistent $GW$ treatments (ev$GW$ and qp$GW$, respectively). All methods consistently increase ionization potentials above the one-shot values. The $\bar{\Delta}GW_{0}$ estimates are higher than ev$GW$, but appear to be closer to qp$GW$. As there are very few published qp$GW$ results for molecules, it is not possible to assess whether this is a general trend for molecules.

For all the studied molecules our method yields results in good agreement with CCSD(T). The improvement is only modest for the smallest molecules, since the QP shift strongly depends on $\varepsilon^{0}$, i.e., it is not constant for all occupied (or unoccupied) states and the assumption of Eq. (12) is not fulfilled. For instance, for N2 $G_{0}W_{0}$ shifts the lowest valence state by -7.34 $\pm$0.06 eV, but the HOMO energy is decreased by -4.63$\pm$0.05 eV. In contrast, the shifts are closer for hexacene: $-2.73\pm 0.08$ and $-1.63\pm 0.06\,{\rm eV}$  for the bottom valence and HOMO states, respectively.

The mean absolute difference between the $\bar{\Delta}GW_{0}$ and CCSD(T) values is 0.20 eV for molecules, but for the largest systems (tetracene and hexacene) it is only $0.05$ eV. This indicates that (i) the scissors operator approximation in Eq. (12) is more appropriate for larger molecules and (ii) $\bar{\Delta}GW_{0}$ self-consistency is more accurate when $G_{0}W_{0}$ is already a good approximation, i.e., in our case it gives results that are sufficiently close to the CCSD(T) values.

III.2 Periodic systems

Next we study self-consistency for several periodic solids listed in Table 2 where we calculate the fundamental band gaps

[TABLE]

The stochastic approach is extended here to treat periodic boundary conditionsVlček et al. (2018). We again employ LDA with Troullier-Martins pseudopotentials and Fourier real-space grids with a spacing $h$ which is sufficiently small that the eigenvalues are converged to $<10$ meV (see Table 2). The method is demonstrated on large supercells with 512 atoms (corresponding to $4\times 4\times 4$ conventional cells with 2048 valence electrons).

As mentioned earlier, the $G_{0}W_{0}$ treatment of such large systems is enabled by the stochastic approach Neuhauser et al. (2014); Vlcek et al. (2017); Vlček et al. (2018), but our self-consistency scheme is applicable to any $G_{0}W_{0}$ implementation that yields the full-frequency or full-time matrix element of the self-energy.

The results in Table 2 show that the one-shot correction yields band gaps that are lower than experimental values, in agreement with previous calculations. Faleev et al. (2004); van Schilfgaarde et al. (2006); Shishkin and Kresse (2007a, b); Bruneval and Gatti (2014) In all cases studied, self-consistency is quickly achieved within 3 or 4 iterations. The resulting fundamental band gaps are enlarged by as much as 0.20 eV,** **and are quite close to experiment. The effect of using $\bar{\Delta}GW_{0}$ on the self-energy curves is illustrated in Fig. 3.

Comparison to previous results (Table 2) shows that the $\bar{\Delta}GW_{0}$ fundamental gaps are overall at least as good as the full eigenvalue self-consistency predictions. In contrast, qp$GW$ band gaps are too high and overestimate experiment by $\sim$10% (also see Ref. van Houcke et al. (2017)).

The $\bar{\Delta}GW_{0}$ results reproduce well the experimental values, with the exception of bulk silicon. Note, however, that we employ $4\times 4\times 4$ conventional cells with $\Gamma$ point sampling. For silicon, this cell, while very large, is still not large enough to reach the bulk limit. We can still compare our result with the experimental $\Gamma-X$ gap which is higher (1.3 eV Tiago et al. (2004)). Therefore, overall, for the set of solids investigated the simplified self-consistency of $\bar{\Delta}GW_{0}$ yields gaps with excellent mean absolute error of 0.03 eV with respect to experiment.

In all the investigated cases the difference between the QP shifts for the bottom and top valence states are small and correlate slightly with $E_{g}$. For BN we observe that in the zeroth iteration ($G_{0}W_{0}$) the bottom valence state and $\varepsilon_{H}$ are shifted by $-0.73\pm 0.04$ and $-0.21\pm 0.03\,{\rm eV}$. However, for Si the shifts are $0.13\pm 0.03$ and $-0.16\pm 0.02$ eV.

The excellent performance of $\bar{\Delta}GW$ is surprising but not fortuitous. The structure of the self-energy curve is dominated by plasmon polesHybertsen and Louie (1986b); Godby and Needs (1989); Aryasetiawan and Gunnarsson (1998); Larson et al. (2013), and the energy of these poles is proportional to the band gap. The main goal of the iterative treatment is to capture the necessary changes in the plasmon energy. We accomplish this goal by employing a relative shift of occupied vs. unoccupied states that acts like a scissors operator (Eq. (12)) that opens up the band gap and leads to the desired increase in the plasmon frequency.

To test the dependence of the self-energy on energy and its implications for the self-consistency, we further performed a set of complementary calculations for two nanorystals, Si35H36 and Si705H300, studied by stochastic $GW$ in the pastNeuhauser et al. (2014); Vlček et al. (2018). For these crystals we did several calculations at different Kohn-Sham energies at fitted the QP correction by a quadratic polynomial:

[TABLE]

For the smaller system the QP correction terms are $a_{2}=0.01$ eV*-1*, $a_{1}=-0.22$ and $a_{0}=-1.13$ eV so $\bar{\Delta}$ is far form being a constant. For the large Si nanocrystal, which is already bulk–like, we find however that $a_{2}=0.00$ eV*-1*, $a_{1}=0.04$ and $a_{0}=-1.52$ eV, indicating only a weak linear dependence. We conclude that the approximation of rigid shifts of all occupied and all unoccupied states is well justified for solids, for which the variation of the QP shift across the occupied states is much smaller than for molecular systems, making the scissors-like assumption appropriate.

IV Summary and Conclusions

In this paper we derived a general form of Green’s function self-consistency in the time domain and introduced its simplified form, which we label $\bar{\Delta}GW_{0}$. The underlying assumption of our method is that the differences between Kohn-Sham eigenvalues and quasiparticle energies are approximately just two constants, one for occupied and one for unoccupied states. We approximate this scissors-like correction by the corrections to the HOMO and LUMO energies.

Our approach is merely an a-posteriori treatment of the time-dependent self-energy matrix. Hence, $\bar{\Delta}GW_{0}$ has essentially no additional computational cost beyond that of a one-shot $G_{0}W_{0}$ calculation. In conjunction with the nearly linear scaling stochastic $G_{0}W_{0}$, it is easily applicable to extremely large systems with thousands of electrons. The combined method is best labeled as stochastic $\bar{\Delta}GW_{0}$ or just abbreviated as stochastic $GW_{0}$.

We tested stochastic $\bar{\Delta}GW_{0}$ on molecules and on periodic semiconductors and insulators with large periodic supercells containing 2048 electrons. The predicted ionization potentials and fundamental band gaps are overall much better than one-shot $G_{0}W_{0}$ values when compared to high-level methods and/or experiments. Our simplified self-consistency treatment is especially appropriate for large molecules and periodic systems, and for the latter it yields a mean absolute error of only 0.03 eV.

The stochastic partially self-consistent $\bar{\Delta}GW_{0}$ approach presented here is both accurate and efficient, opening the door to many future applications in chemistry, physics nano- and material sciences.

Acknowledgements.

V.V., E.R. and D.N. were supported by the Center for Computational Study of Excited State Phenomena in Energy Materials (C2SEPEM) at the Lawrence Berkeley National Laboratory, which is funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division under contract No. DEAC02-05CH11231 as part of the Computational Materials Sciences Program. R.B. is grateful for support by the Binational Science Foundation, Grant 2015687 and for support from the Israel Science Foundation – FIRST Program, Grant No. 1700/14. The calculations were performed as part of the XSEDE computational Project No. TG-CHE170058 Towns et al. (2014).