We show that a rigid scissors-like *GW* self-consistency approach, labeled here $\Delta \xafGW0$, can be trivially implemented at zero additional cost for large scale one-shot *G*_{0}*W*_{0} calculations. The method significantly improves one-shot *G*_{0}*W*_{0} and for large systems is very accurate. $\Delta \xafGW0$ is similar in spirit to ev*GW*_{0} where the self-consistency is only applied on the eigenvalues entering Green’s function, while both *W* and the eigenvectors of Green’s function are held fixed. $\Delta \xafGW0$ further assumes that the shift of the eigenvalues is rigid scissors-like so that all occupied states are shifted by the same amount and analogously for all the unoccupied states. We show that this results in a trivial modification of the time-dependent *G*_{0}*W*_{0} self-energy, enabling an *a posteriori* self-consistency cycle. The method is applicable for our recent stochastic-*GW* approach, thereby enabling self-consistent calculations for giant systems with thousands of electrons. The accuracy of $\Delta \xafGW0$ increases with the system size. For molecules, it is up to 0.4-0.5 eV away from coupled-cluster single double triple (CCSD(T)), but for tetracene and hexacene, it matches the ionization energies from both CCSD(T) and ev*GW*_{0} to better than 0.05 eV. For solids, as exemplified here by periodic supercells of semiconductors and insulators with 6192 valence electrons, the method matches ev*GW*_{0} quite well and both methods are in good agreement with the experiment.

## I. INTRODUCTION

The *GW* approximation^{1} to many-body perturbation theory is often used to calculate electron removal or addition energies and the related (inverse) photoemission spectra of molecules, nanostructures, and bulk materials.^{2–15} *GW* is part of a family of methods that describe the probability amplitude of a quasiparticle (QP) to propagate between two space-time points $r,t$ and $r\u2032,t\u2032$ with Green’s function $Gr,r\u2032,t,t\u2032$, the poles of which are the QP energies. Green’s function is obtained perturbatively from reference (non-interacting) Green’s function, $G0r,r\u2032,t,t\u2032$, via a Dyson equation (all equations use atomic units),

where $\Sigma r1,r2,t$ is the time-dependent self-energy. Reference Green’s function is typically^{16,17} given by the Kohn-Sham^{18} (KS) density function theory (DFT).^{19} In *GW*, the self-energy Σ is approximated as

where $Wr,r\u2032,t$ is the screened Coulomb interaction, usually evaluated within the random phase approximation (RPA).

A solution of the equations above requires in principle a self-consistent procedure since both $\Sigma r,r\u2032,t$ and $Wr,r\u2032,t+$ depend on $Gr,r\u2032,t$. In practice, the self-consistency is often abandoned and the most common *GW* treatment is based on “one-shot” scheme,^{16,17} *G*_{0}*W*_{0}, where the right hand side of Eq. (2) becomes $G0r,r\u2032,tW0r,r\u2032,t$ and *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*).^{8,20–25} A fully self-consistent solution is computationally extremely demanding,^{10,21,26–33} and in many situations, it yields results that are worse than *G*_{0}*W*_{0};^{8,34,35} heuristically, since *GW* is an approximation, an iteration of that approximation may enhance the deviation from the exact result.

To simplify the problem, a static and Hermitian approximation to the self-energy^{20,36} is sometimes used in the so-called QP self-consistent *GW* (qp*GW*), which iteratively updates the Σ and QP wave-functions (i.e., the Dyson orbitals). The qp*GW* approach is still computationally expensive and cannot be applied for large systems; furthermore, it tends to overestimate *E*_{g}^{35–37} and the ionization potentials^{15,38,39} due to an overly strong screened Coulomb interaction.^{15,35,36} Alternatively, the *W* term is kept “frozen” and self-consistency is sought only in Green’s function.^{8,27,28,35,38} This method is termed eigenvalue self-consistent *GW* (ev*GW*_{0}), and it was applied successfully to bulk systems^{8,35,40} and to organic molecules.^{25,27,28} Despite the lack of self-consistency in *W*, which leads to overscreening (since it is based on an underestimated KS DFT gap)^{35} and non-negligible starting point dependence,^{41,42} ev*GW*_{0} provides a reasonable route for improving the “one-shot” predictions or gaps in extended systems, notwithstanding the fact that for fairly small molecular systems the starting-point dependence is very large.

While overall self-consistency methods provide an improvement over one-shot *G*_{0}*W*_{0}, there is an issue of scaling and overall cost. Even one-shot *G*_{0}*W*_{0} has usually been quite expensive, scaling as $O(Ne3)\u2212O(Ne4)$ with the number of electrons *N*_{e}, and self-consistency methods at least multiply the cost by the number of iterations. Thus, in spite of their potential, *G*_{0}*W*_{0} and its self-consistent improvements were only applicable to fairly small systems unless drastically loose convergence criteria are made and/or huge computational efforts are used, and even then large systems with many thousands of electrons were not feasible.

Recently, however, we developed an effectively linear-scaling stochastic approach to *G*_{0}*W*_{0} with moderate computational cost.^{43–45} The approach has been verified to give the same results as deterministic converged-basis (or converged plane-wave) *G*_{0}*W*_{0} simulations and can handle systems with more than 10 000 electrons.^{45} But as stochastic *G*_{0}*W*_{0} obtains directly only a column of the self-energy matrix rather than the full-matrix, most self-consistent improvements cannot be applied to it.

Here, we enhance the power of the stochastic approach, by showing that a variant of ev*GW*_{0} whereby the energies of the HOMO and LUMO are shifted rigidly (a scissors approach, labeled below $\Delta \xafGW0$) can be easily applied on top of *G*_{0}*W*_{0}, without any additional cost. Thus stochastic *GW* can be enhanced beyond *G*_{0}*W*_{0}; the combined approach is labeled stochastic $\Delta \xafGW0$.

Specifically, we first develop the method showing how the application of a scissors operator for Green’s function is an *a posteriori* step. The $\Delta \xafGW0$ formulation requires only the knowledge and use of two quasiparticle states, rather than the whole quasiparticle spectrum. As long as one has access to the relevant two diagonal matrix elements (usually the HOMO and LUMO elements) of the self-energy at all frequencies or all times, then, irrespective of the system size, the computational cost of the self-consistency step is negligible (i.e., seconds on a single-core machine), so $\Delta \xafGW0$ costs no more than *G*_{0}*W*_{0}.

Equally important, we show that stochastic $\Delta \xafGW0$ improves with increasing system size. First, we show for large molecules that the ionization energies from $\Delta \xafGW0$ converge to the ev*GW*_{0} predictions and both methods approach for large systems the more exact coupled-cluster single double triple (CCSD(T)) results.

Next, we perform stochastic *G*_{0}*W*_{0} calculations for periodic semiconductors and insulators using large supercells with 6192 valence electrons. For solids, $\Delta \xafGW0$ consistently increases *E*_{g} toward better agreement with ev*GW*_{0} and the experiment (with resulting mean absolute error relative to the experiment of better than 0.1 eV). In all cases, the self-consistency is reached in very few iterations without any additional cost beyond the *G*_{0}*W*_{0} calculation.

## II. THEORY

### A. Green’s function self-consistency in the time domain

The QP energy of the *i*th state $\epsilon QPi$ is calculated using the usual form of the perturbative *GW* approximation in which the Kohn-Sham eigenvalues (*ε*^{0}) are corrected by the QP shift (Δ) using a fixed point equation,

where

Here, $vxc$ is the Kohn-Sham exchange-correlation potential for the DFT density, $\Sigma \u0303i\omega $ is the Fourier transform of the matrix-element of $\Sigma ^(t)$,

and $\Sigma ^(t)$ is given by Eq. (2).

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

where Tr denotes a trace over all KS states, *μ* is the chemical potential, *θ* is the Heaviside step function that guarantees forward and backward time propagation for particles and holes, respectively, and *ĥ*^{0} is the KS Hamiltonian,

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;^{43,44} this is equivalent to using the RPA approximation for *W*_{0}.

In the time-domain, the self-energy matrix element for the *i*th state is

Finally, after Fourier transformation combined with time-ordering,^{43,44} the “one-shot” QP energy is calculated through Eq. (4).

In the ev*GW*_{0} procedure, 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

where $\Delta ^$ contains all the many-body contributions to the spectrum.

As common in ev*GW*_{0} self-consistency,^{8,35,38} for the purpose of Eq. (9), one ignores the non-Hermiticity of the self-energy and treats it as diagonal so that in the KS basis the correction operator is approximated as

Hence, all the KS energies in the exponent in Eq. (9) are shifted to the QP energies obtained from Eq. (3). Equation (10) is the fundamental equation of ev*GW*_{0}.

Since $\Delta ^$ in Eq. (9) is by construction diagonal in the KS basis set, we view it as a function of the KS Hamiltonian $\Delta \xaf\u01250$ that interpolates all QP shifts. Green’s function is therefore

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

### B. Efficient implementation

In many cases, $\Delta \xaf$ is well described by a low degree polynomial with discontinuity only at the bandgap energies, $\epsilon H$ and $\epsilon L,$ corresponding to the highest occupied (*H*) and lowest unoccupied (*L*) states, respectively. The zeroth order term in this polynomial is a scissors operator^{46,47} which shifts the occupied and unoccupied states down and up in energy, respectively,

We call this approximation $\Delta \xafGW0$ and use it in Sec. III for molecules and periodic systems. Furthermore, rather than choosing highest occupied and lowest unoccupied states, $Re\Delta \xaf$ may be taken based on an *average* shift of all occupied and unoccupied states, respectively.

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

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

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 uses^{46,47} of the scissors-operator in *GW* which require repeated evaluations of the self-energy. Furthermore, the $\Delta \xafGW0$ approach is applicable to any implementation which yields $\Sigma \u0303i\omega $ for the HOMO and LUMO states. Thus, this method is naturally suited for the stochastic *G*_{0}*W*_{0} method^{43–45} which provides the self-energy over the full-time domain and therefore over a wide range of frequencies (spanning several hundred eV). Note that while in principle one could use higher order polynomials of Δ, and such polynomials would generally require going beyond a post-processing method, i.e., would require additional *GW* calculations.^{48} Therefore, the simple equation (14), for a scissors-type expression, is our main result.

## III. RESULTS AND DISCUSSION

### A. Molecules

We first test our approach on ionization potentials *I* [taken as the negative of the HOMO energy, −*ε*(*H*)] for a set of small molecules listed in Table I. A ground state DFT calculation is performed using a Fourier real-space grid, ensuring (through the Martyna-Tuckerman approach)^{49} that the potentials are not periodic. The exchange-correlation interaction is described by the local density approximation (LDA)^{50} with Troullier-Martins pseudopotentials;^{51} the DFT eigenvalues are converged up to <10 meV with respect to the spacings of the real space grids (given in Table I).

. | . | I (eV)
. | ||||||||
---|---|---|---|---|---|---|---|---|---|---|

System . | h (a_{0})
. | LDA . | G_{0}W_{0}
. | $\Delta \xafGW0$ . | evGW_{0}
. | qpGW
. | CCSD(T) . | Expt. . | $\delta (\Delta \xafGW0)$ . | δ(evGW_{0})
. |

Nitrogen | 0.35 | 10.44 | 15.08 (0.05) | 15.93 (0.05) | 15.32^{a} | 16.01^{a} | 15.57^{b} | 15.58 | 0.85 | 0.43 |

Ethylene | 0.35 | 6.92 | 10.50 (0.04) | 10.87 (0.04) | 10.24^{a} | 10.63^{a} | 10.67^{b} | 10.68 | 0.37 | −0.09 |

Urea | 0.30 | 6.10 | 9.53 (0.08) | 10.48 (0.08) | 9.81^{a} | 10.45^{a} | 10.05^{b} | 10.28 | 0.95 | 0.49 |

Naphthalene | 0.35 | 5.71 | 8.10 (0.09) | 8.39 (0.09) | 8.15^{c} | … | 8.25^{c} | 8.14 | 0.29 | 0.19 |

Tetracene | 0.35 | 4.89 | 6.79 (0.08) | 6.94 (0.08) | 6.84^{c} | … | 7.02^{c} | 6.97 | 0.15 | 0.16 |

Hexacene | 0.35 | 4.52 | 6.15 (0.06) | 6.33 (0.06) | 6.19^{c} | … | 6.32^{c} | 6.33 | 0.18 | 0.15 |

. | . | I (eV)
. | ||||||||
---|---|---|---|---|---|---|---|---|---|---|

System . | h (a_{0})
. | LDA . | G_{0}W_{0}
. | $\Delta \xafGW0$ . | evGW_{0}
. | qpGW
. | CCSD(T) . | Expt. . | $\delta (\Delta \xafGW0)$ . | δ(evGW_{0})
. |

Nitrogen | 0.35 | 10.44 | 15.08 (0.05) | 15.93 (0.05) | 15.32^{a} | 16.01^{a} | 15.57^{b} | 15.58 | 0.85 | 0.43 |

Ethylene | 0.35 | 6.92 | 10.50 (0.04) | 10.87 (0.04) | 10.24^{a} | 10.63^{a} | 10.67^{b} | 10.68 | 0.37 | −0.09 |

Urea | 0.30 | 6.10 | 9.53 (0.08) | 10.48 (0.08) | 9.81^{a} | 10.45^{a} | 10.05^{b} | 10.28 | 0.95 | 0.49 |

Naphthalene | 0.35 | 5.71 | 8.10 (0.09) | 8.39 (0.09) | 8.15^{c} | … | 8.25^{c} | 8.14 | 0.29 | 0.19 |

Tetracene | 0.35 | 4.89 | 6.79 (0.08) | 6.94 (0.08) | 6.84^{c} | … | 7.02^{c} | 6.97 | 0.15 | 0.16 |

Hexacene | 0.35 | 4.52 | 6.15 (0.06) | 6.33 (0.06) | 6.19^{c} | … | 6.32^{c} | 6.33 | 0.18 | 0.15 |

Six molecules were tested, ranging in size from N_{2} to hexacene. In all cases, stochastic *G*_{0}*W*_{0}^{43–45} was used to calculate the self-energy. $\Delta \xafGW0$ converges in 3-4 iterations after the initial *G*_{0}*W*_{0} calculation; the self-energy curves are illustrated for hexacene in Fig. 2. Since our self-consistency procedure is performed *a posteriori*, its computational cost is negligible (less than a second on a single core machine). In all cases, $\Delta \xafGW0$ increases the ionization potentials above the one-shot values. This is a desired outcome since *G*_{0}*W*_{0} underestimates *I*.^{24,42,44}

Table I shows the ionization energies at several levels, including *G*_{0}*W*_{0}, $\Delta \xafGW0$, deterministic self-consistent *GW* on two levels, CCSD(T), and also the experiment.

The *G*_{0}*W*_{0} column shows our stochastic *G*_{0}*W*_{0} values that are fully converged with grid size (i.e., in a plane-wave language, corresponding to a fully converged plane-wave cutoff) and are therefore close to the converged deterministic plane-wave basis results benchmarked in Ref. 52 as well as basis set extrapolated values obtained with deterministic codes, with differences of typically less than 0.1 eV, which are mostly due to the use of an LDA functional rather than PBE exchange-correlation functional.^{44}

Our partially self-consistent $\Delta \xafGW0$ is then compared in the table to other partially self-consistent techniques in the literature, suitable for small molecules. The most direct comparison should be to ev*GW*_{0}, where, as mentioned, Green’s function is also converged (leaving *W*_{0} frozen at its initial value), for an initial Hamiltonian with a local potential. There are large deviations for the small molecules (0.6–0.7 eV), and there are several contributing factors to these differences. First, for the three small molecules about 0.1-0.2 eV difference is due to the non-completeness of the basis set (see Ref. 24 for the effect of basis-set extrapolation). For nitrogen and urea, another factor of ≈0.1 eV is due to our use of LDA rather the PBE in the quoted ev*GW*_{0} results.^{44} Taken together, this still leaves about an error of ≈0.4 eV between our results and ev*GW*_{0} for the small molecules.

The difference between the present $\Delta \xafGW0$ and ev*GW*_{0} is made more explicit by the last two columns in Table I, which shows the degree of correction due to self-consistency (i.e., the difference between the final result and the one-shot *G*_{0}*W*_{0}) for our calculation and ev*GW*_{0}. For the small molecules, $\Delta \xafGW0$ gives a correction which is higher than ev*GW*_{0} by more than 0.4 eV. But for larger systems, the deviation is getting much smaller. Specifically, the deviation from $\Delta \xafGW0$ decreases to 0.1–0.25 eV for the three larger molecules (three acenes) where the deterministic basis-set results for the acenes were already extrapolated to the infinite basis-set limit. A small portion of the difference is due to our use of LDA rather than PBE. Overall the much better agreement of the rigid-shift approach with ev*GW*_{0} for large systems is due to the improved assumption of a rigid shift.

A separate question from the agreement of the rigid-shift approximation with ev*GW*_{0} is how accurate the overall method is, i.e., does the method agree with more exact non-GW approaches. Here, the comparison with the CCSD(T) ionization energy (extrapolated to the limit of infinite basis-set) is quite encouraging. For the small molecules, there is a large deviation from CCSD(T), but it decreases dramatically to 0.15 eV and below once we go to the larger acenes.

Interestingly, $\Delta \xafGW0$ is generally closer to the CCSD(T) results than ev*GW*_{0}. This is probably due to error cancellation; i.e., ev*GW*_{0} based on a local potential somewhat underestimates the CCSD(T) results, and $\Delta \xafGW0$ generally somewhat overestimates ev*GW*_{0} due to the rigid shift assumption. Nevertheless, Table I does show that for the ionization energy, as the system size grows the essentially exact CCSD(T) is better and better predicted by both the ev*GW*_{0} approximation and our $\Delta \xafGW0$ further approximation to ev*GW*_{0}. Likewise, $\Delta \xafGW0$ agrees well with the experiment for the larger molecules (the three acenes), which makes sense since except for the effect of vibrational energy change CCSD(T) should reproduce the experiment.

To understand the improvement of ionization energies from $\Delta \xafGW0$ with the system size, we examine the dependence of the state dependence shift, Δ(*i*), on the state energy *ε*^{0}(*i*). The large errors for small molecules stem from the strong dependence of Δ(*i*) on *ε*^{0}. Namely, the QP shift is far from being constant for all occupied (or unoccupied) states, so Eq. (12) is not fulfilled. For the smallest molecule we tested, N_{2}, *G*_{0}*W*_{0} shifts the lowest valence state by −7.34 ± 0.06 eV, while the HOMO energy is decreased on this level by only −4.63 ± 0.05 eV; i.e., the difference in QP shifts across the valence band is very high, 2.60 eV. But for larger systems, the scissors operator approximation in Eq. (12) is getting gradually better. Thus, for hexacene, the difference in QP shifts between the bottom valence and the HOMO state is only 0.87 eV, three times smaller than that for *N*_{2}.^{53}

The improvement of the rigid-shift approximation with the system size is important since it implies that for large system sizes, for which other self-consistent techniques are too expensive, the $\Delta \xafGW0$ approach would be sufficiently accurate.

We conclude this section by a cautionary note. The success of both ev*GW*_{0} based on a local potential and $\Delta \xafGW0$ is for the ionization potential. For the LUMO, i.e., for the electron affinity, the results of both ev*GW*_{0} and $\Delta \xafGW0$ are not as good (deviations of more than 0.5 eV from the exact results even for the larger systems). For the LUMO, a time dependent DFT (TDDFT)-based kernel may give better one-shot results,^{43,54} and the application of rigid-shift self-consistency for the LUMO would be described in a future publication.

### B. Periodic systems

Next we study self-consistency for periodic solids, using an extension of the stochastic GW approach to treat periodic boundary conditions.^{45} The self-consistent $\Delta \xafGW0$ method is demonstrated on large supercells with 1712 atoms (corresponding, for the system we study below, to 6 × 6 × 6 conventional cells with 6912 valence electrons). The stochastic approach handles well such large system sizes,^{43–45} yielding again self-energies for selected electronic states. Here we are interested in the fundamental bandgap, *E*_{g} = *ε*(*L*) − *ε*(*H*), so we need both the highest occupied and lowest unoccupied states. For large extended systems, the complications associated with a TDDFT vs. RPA based kernel are tamed in comparison with molecules, so it makes sense to study the gap, based solely on RPA.^{57}

Table II shows for several solids the one-shot and $\Delta \xafGW0$ and compares them again to previous self-consistent *GW* calculations and to experiment. The results demonstrate first that the one-shot correction yields bandgaps that are mostly (though not always) lower than experimental values, in agreement with previous calculations.^{8,20,35,36}

. | . | E_{g} (eV)
. | ||||||||
---|---|---|---|---|---|---|---|---|---|---|

System . | h (a_{0})
. | LDA . | G_{0}W_{0} (present)
. | G_{0}W_{0}^{a}
. | $\Delta \xafGW0$ . | evGW^{a}
. | qpGW^{a}
. | Expt. . | $\delta (\Delta \xafGW0)$ . | δ(evGW_{0})
. |

Si | 0.446 | 0.50 | 1.24 (0.04) | 1.12 | 1.32 (0.04) | 1.20 | 1.28 | 1.17^{b} | 0.08 | 0.08 |

SiC | 0.293 | 1.36 | 2.40 (0.04) | 2.27 | 2.50 (0.04) | 2.43 | 2.64 | 2.42^{c} | 0.10 | 0.16 |

AlP | 0.368 | 1.46 | 2.45 (0.04) | 2.44 | 2.58 (0.04) | 2.59 | 2.77 | 2.51^{d} | 0.13 | 0.15 |

C | 0.336 | 4.15 | 5.60 (0.04) | 5.50 | 5.71 (0.04) | 5.68 | 5.99 | 5.48^{e} | 0.11 | 0.18 |

BN | 0.380 | 4.35 | 6.11 (0.04) | 6.10 | 6.28 (0.04) | 6.35 | 6.73 | 6.1–6.4^{f} | 0.17 | 0.25 |

. | . | E_{g} (eV)
. | ||||||||
---|---|---|---|---|---|---|---|---|---|---|

System . | h (a_{0})
. | LDA . | G_{0}W_{0} (present)
. | G_{0}W_{0}^{a}
. | $\Delta \xafGW0$ . | evGW^{a}
. | qpGW^{a}
. | Expt. . | $\delta (\Delta \xafGW0)$ . | δ(evGW_{0})
. |

Si | 0.446 | 0.50 | 1.24 (0.04) | 1.12 | 1.32 (0.04) | 1.20 | 1.28 | 1.17^{b} | 0.08 | 0.08 |

SiC | 0.293 | 1.36 | 2.40 (0.04) | 2.27 | 2.50 (0.04) | 2.43 | 2.64 | 2.42^{c} | 0.10 | 0.16 |

AlP | 0.368 | 1.46 | 2.45 (0.04) | 2.44 | 2.58 (0.04) | 2.59 | 2.77 | 2.51^{d} | 0.13 | 0.15 |

C | 0.336 | 4.15 | 5.60 (0.04) | 5.50 | 5.71 (0.04) | 5.68 | 5.99 | 5.48^{e} | 0.11 | 0.18 |

BN | 0.380 | 4.35 | 6.11 (0.04) | 6.10 | 6.28 (0.04) | 6.35 | 6.73 | 6.1–6.4^{f} | 0.17 | 0.25 |

Moving to $\Delta \xafGW0$, in all cases studied, self-consistency is quickly achieved within 2 or 3 iterations and the self-energy curves are illustrated in Fig. 3. The fundamental bandgaps are enlarged by as much as 0.17 eV compared to *G*_{0}*W*_{0}. For all the solids in Table II, the difference between the QP shifts for the bottom and top valence states is small (correlating slightly with *E*_{g}) and similar in magnitude for all states.^{58} The underlying assumption that all states in either the occupied or unoccupied subspace are shifted (approximately) by a constant is fulfilled well.

Gaps predicted by $\Delta \xafGW0$ are similar to previously published ev*GW* results.^{59} Much of the difference is due to the use of distinct atomic potentials in previous calculations, which are known to cause energy differences up to ∼0.2 eV.^{60} More importantly for our purpose, here are the energy differences due to the self-consistency (i.e., the differences of either ev*GW*_{0} or $\Delta \xafGW0$ energies from *G*_{0}*W*_{0}). These differences (which are less sensitive to the particular type of potentials used) are very close, to within about 0.04 eV in average.

In Table II, we also list experimental bandgaps. Such comparison has to be performed with caution since our calculations do not account for renormalization of *E*_{g} due to zero-point motion,^{61} which decreases the gap by as anywhere from 0.05 eV (silicon) to 0.37 eV (diamond).^{62} Yet it is encouraging to see that $\Delta \xafGW0$ enlarges *E*_{g} and that our predictions are quite close to experimental values.

## IV. SUMMARY AND CONCLUSIONS

We derived here a general form of Green’s function self-consistency in the time domain and introduced its simplified form, which we label $\Delta \xafGW0$. The starting point is the scissors correction, where the differences between Kohn-Sham eigenvalues and quasiparticle energies are approximately just two constants, one for occupied and the other for unoccupied states. We approximate this scissors-like correction by the corrections to the HOMO and LUMO energies. But unlike previous self-consistency schemes (even those relying on scissors-like operators), our approach does not require repeated many-body calculation and uses only two self-energy matrix elements. Practically, it is merely an *a posteriori* treatment of the time-dependent self-energy matrix. Yet, it captures the major changes in the structure of Green’s function and hence of the self energy for large systems.

$\Delta \xafGW0$ 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 can treat extremely large systems with thousands of electrons. The combined method is best labeled as stochastic $\Delta \xafGW0$.

We tested stochastic $\Delta \xafGW0$ on molecules and on periodic semiconductors and insulators with large periodic supercells with 6192 electrons. The predicted ionization potentials and fundamental bandgaps are overall 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, for which no other computationally affordable approach existed up to now.

$\Delta \xafGW0$ is an approximation on top of another, successful approximation, ev*GW*_{0}. For small systems, both approximations have faults (dependence on initial DFT Hamiltonian in the final results of ev*GW*_{0}, compounded by strong state dependence of the self-energy diagonal elements violating the underlying assumption on $\Delta \xafGW0$). But the purpose of the new approach is to address very large systems, for which these two issues are much less severe, and for which $\Delta \xafGW0$ is automatically suited, since, unlike even ev*GW*_{0}, it does not require the computation of new operators when the quasi-particle energies are updated.

The stochastic partially self-consistent $\Delta \xafGW0$ approach presented here is both accurate for large systems and efficient, opening the door to many future applications in chemistry, physics, and nano- and material sciences.

## ACKNOWLEDGMENTS

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 from the Binational Science Foundation, under Grant No. 2015687, and for support from the Israel Science Foundation—FIRST Program, under Grant No. 1700/14. The calculations were performed as part of XSEDE computational Project No. TG-CHE170058.^{68}

## REFERENCES

Actually, a first order polynomial could also be incorporated, with some additional storage steps, as will be detailed in a future publication.

The *G*_{0}*W*_{0} shifts for hexacene are −2.73 ± 0.08 and −1.63 ± 0.06 eV for the bottom valence and HOMO states.

Of course, just like for molecules, 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.

For BN, we observe that in the zeroth iteration (*G*_{0}*W*_{0}) the bottom valence state and *ε*_{H} are shifted by −0.73 ± 0.04 and −0.21 ± 0.03 eV, respectively. However, for Si, the shifts are only 0.13 ± 0.03 and −0.16 ± 0.02 eV.