The increasingly popular GW method is becoming a convenient tool to determine vertical ionization energies in molecular systems. However, depending on the formalism used and the range of orbitals investigated, it may be hampered by a steep computational scaling. To alleviate this issue, correlated natural virtual orbitals (NVOs) based on second-order Møller–Plesset (MP2) and direct MP2 correlation energies are implemented, and the resulting correlated NVOs are tested on GW quasiparticle energies. Test cases include the popular GW variants G0W0 and evGW0 as well as more elaborate vertex corrections. We find that for increasingly larger molecular systems and basis sets, NVOs considerably improve efficiency. Furthermore, we test the performance of the truncated (frozen) NVO ansatz on the GW100 test set. For the latter, it is demonstrated that, using a carefully chosen truncation threshold, NVOs lead to a negligible loss in accuracy while providing speedups of one order of magnitude. Furthermore, we compare the resulting quasiparticle energies to very accurate vertical ionization energies obtained from coupled-cluster theory with singles, doubles, and noniterative triples [CCSD(T)], confirming that the loss in accuracy introduced by truncating the NVOs is negligible compared to the methodical errors in the GW approximation. It is also demonstrated that the choice of basis set impacts results far more than using a suitably truncated NVO space. Therefore, at the same computational expense, more accurate results can be obtained using NVOs. Finally, we provide improved reference CCSD(T) values for the GW100 test set, which have been obtained using the def2-QZVPP basis set.

## I. INTRODUCTION

Correlated natural virtual orbitals (NVO)^{1} are an emerging class of orbitals that diagonalize a correlated density matrix spanning the full virtual orbital space. Most prominently, second-order Møller–Plesset perturbation theory (MP2) has been used for the construction of the virtual part of the density matrix.^{2,3} Another viable alternative that can provide correlated density matrices at a rather low cost is the random-phase approximation (RPA).^{4} It has further been shown that using only a truncated subset of these MP2-based NVOs is a well-suited approach for obtaining correlation energies or properties from higher-order correlated methods, allowing for truncation of the space spanned by the virtual orbitals. MP2-based NVOs have since been applied to coupled-cluster (CC) theory with singles–and–doubles (CCSD),^{5} with perturbative triples [CCSD(T)],^{6–11} and in equation-of-motion-based CC (EOM-CC) and linear-response CC methods.^{12–14} Recently, MP2-based NVOs have also been applied to two- and four-component relativistic electron-correlation methods.^{15}

Another interesting application for correlated NVOs is the *GW* approximation for molecules.^{16–18} In the *GW* approximation, the self-energy is calculated to correct a set of Kohn–Sham (KS) orbital energies to yield quasiparticle energies. The latter are closely related to vertical ionization energies (VIEs). Obtaining the correlation part of the self-energy is however not a trivial task. Within the random-phase approximation, a naive ansatz scales as *N*^{6} with molecular size, which is roughly the same scaling behavior as CCSD. While certain techniques have been developed to lower this high computational effort,^{19–24} they are usually specific to a single variant of *GW*, namely, with the response calculated in the direct RPA approximation (dRPA), neglecting all exchange terms. Furthermore, these approximations are often limited to only a few quasiparticle states, either drastically losing either efficiency or accuracy if full quasiparticle spectra are to be extracted. Also, specific vertex corrections are not accessible using these low-scaling *GW* variants. In many molecular applications, a *GW* calculation is only the beginning of a larger series of calculations, including follow-up calculations. Typical examples of follow-up calculations are simulations of optical excitation spectra using the Bethe–Salpeter equation,^{22,25–29} transport simulations,^{30–33} or computations of further molecular properties.^{34–36} When these types of calculations are performed *a posteriori*, it is questionable to only use partial quasiparticle spectra, and instead full ones are desirable. Therefore, the *GW* step is often the time-limiting step, especially if more than a few quasiparticle states are to be computed. Reducing the virtual expansion space can then significantly speed up the *GW* step, which has been recognized early as a viable path to obtaining fast yet reliable results. Previous techniques employed to reduce the virtual expansion space in *GW*, for example, focused on directly optimizing the virtual orbital space.^{37} In solid-state applications, attempts to reduce the expansion space in a *GW* calculations have also been applied successfully.^{38,39} A straightforward reduction of the virtual space is however connected to many caveats, as the *GW* self-energy generally converges only slowly with respect to the number of empty states. Therefore, truncation-based approaches are usually accompanied by extrapolation techniques or other error compensation approaches.^{37,38} Still, the convergence of truncated virtual spaces is still problematic.^{40,41} Like coupled-cluster theory, *GW* methods would see significant benefits from being able to truncate the virtual space drastically. In this study, we will demonstrate that correlated NVOs provide a convenient route toward truncating the virtual space. Viable paths to use truncated NVO spaces in *GW* calculations will be outlined, and the impact of certain truncation parameters will be discussed based on carefully assessed test sets.

## II. THEORY

Unlike in correlated wave-function theory, as, for example, in coupled-cluster or configuration-interaction theory, the starting point for a *GW* calculation is usually a Kohn–Sham reference Slater determinant, not a Hartree–Fock (HF) reference determinant. Therefore, in Sec. II, all orbitals are assumed to be obtained from a converged KS calculation. While certainly there is a strong dependence of the *GW* results on the functional used to optimize the orbitals,^{42,43} in this study, we will show that the performance of NVOs is similar for all functional approximations.

### A. Natural virtual orbitals from second-order perturbation theory

^{44}First defined only for occupied orbitals, it was shown that post-HF calculations using natural orbitals could be brought to convergence significantly faster than with canonical orbitals. As shown by Davidson,

^{1}natural orbitals minimize the least square error between the canonical density matrix and its truncated approximate

^{45–47}natural transition orbitals (NTOs),

^{48}and MP2 natural orbitals.

^{49,50}While NTOs are mainly used for the calculation of excited state properties,

^{51}PNOs and MP2-NOs are primarily used for ground state calculations. Whereas PNOs are constructed for pairs of occupied HF orbitals, MP2-NOs are constructed from the total first-order wave function.

^{49,50,52,53}is the most central step. This virtual–virtual block is given by

*ɛ*

_{pq}=

*ɛ*

_{p}+

*ɛ*

_{q}denotes the sum of two orbital energies from two canonical spin-orbitals, and (

*pq*|

*rs*) is a two-electron integral in Mulliken notation. We denote occupied spin-orbitals with the indices

*i*,

*j*,

*k*, …, virtual spin-orbitals with

*a*,

*b*,

*c*, …, and general spin-orbitals with

*p*,

*q*,

*r*, …. Note that in our implementation, the computed occupation numbers (i.e., eigenvalues) refer to the spin-orbital-based density matrix in Eq. (2) not only in the case of an unrestricted Hartree–Fock (UHF) or Kohn–Sham (UKS) reference determinant but also in the case of a restricted Hartree–Fock (RHF) or Kohn–Sham (RKS) reference determinant.

The NVOs are then obtained from the eigenvectors **N** of the virtual–virtual density matrix **D**^{MP2} of Eq. (2). We keep all vectors $N\u0303$ with occupation numbers larger than some prescribed threshold. In the next step, we build semi-canonical NVOs by diagonalizing the Fock matrix in the truncated space. The essential steps for the determination of the semi-canonical NVOs, collected in the matrix $C\u0303$, are summarized in Table I. The orbital energies of the semi-canonical NVOs are the diagonal elements of the matrix $\epsilon \u0303$. The matrix **C** contains the original HF or KS virtual spin-orbitals.

(1) Computation of the virtual–virtual MP2 density matrix | D^{MP2} |

(2) Diagonalization of the density matrix | N^{†}D^{MP2}N = n |

(3) Selection of a reduced number of eigenvectors | $N\u2192N\u0303$ |

(4) Construction of Fock matrix in truncated space | $F\u0303=N\u0303\u2020\epsilon N\u0303$ |

(5) Diagonalization of the Fock operator in the truncated space | $U\u2020F\u0303U=\epsilon \u0303$ |

(6) Construction of semi-canonical orbitals | $C\u0303=CN\u0303U$ |

(1) Computation of the virtual–virtual MP2 density matrix | D^{MP2} |

(2) Diagonalization of the density matrix | N^{†}D^{MP2}N = n |

(3) Selection of a reduced number of eigenvectors | $N\u2192N\u0303$ |

(4) Construction of Fock matrix in truncated space | $F\u0303=N\u0303\u2020\epsilon N\u0303$ |

(5) Diagonalization of the Fock operator in the truncated space | $U\u2020F\u0303U=\epsilon \u0303$ |

(6) Construction of semi-canonical orbitals | $C\u0303=CN\u0303U$ |

Note that the density matrix becomes block diagonal if the wave function is chosen to be an eigenfunction of the spin operator $S\u0302z$, as it is done in restricted or unrestricted HF/KS calculations. Then, the off-diagonal *α*–*β* blocks are zero, i.e., (**D**_{αβ} = **D**_{βα} = **0**). Accordingly, the density matrix can be diagonalized separately for the *α*–*α* and *β*–*β* blocks. The computational effort can be reduced even further in the closed-shell RHF/RKS case by exploiting the symmetry **D**_{αα} = **D**_{ββ}. Otherwise, steps (1)–(5) must be carried out separately for the two spin cases. If the wave function is neither an eigenfunction of $S\u0302z$ nor $S\u03022$, as it is the case when spin–orbit effects are included, also the off-diagonal blocks must be included in the calculation of **D**^{MP2}. In the latter case, the resulting complex NVOs are suitable to be applied in two- or four-component relativistic calculations, or calculations in external magnetic fields.^{22,54–57} Direct MP2-based NVOs (denoted as dNVOs) can be obtained similarly by simply neglecting the exchange integrals when generating the MP2 density matrix using Eq. (2).^{58} Note that a prefactor of two has to be applied additionally in this case. Equation (2) also determines the overall scaling of the NVO generation step, scaling as *N*^{5}. For dNVOs, the scaling can be reduced to *N*^{4} using Laplace transformations.^{59} While MP2 density matrices are a convenient choice in localized orbital simulations, they are less common in periodic frameworks. For the latter, density matrices based on the RPA may be better suited. Virtual–virtual RPA density matrices, for example, occur in self-consistent RPA,^{60} and their calculation scales similar to Laplace-transformed direct MP2.^{61}

### B. Correlation self-energies from the *GW* method within the random-phase approximation

**Ω**is a matrix containing all positive dRPA excitation energies (see, for example, Ref. 62). The latter can be obtained from the eigenvalue problem

**X**and

**Y**satisfy the condition

**X**

^{†}

**X**−

**Y**

^{†}

**Y**=

**1**. Within the dRPA, the orbital rotation matrices

**A**and

**B**are given as

*ɛ*

_{p}refers to a KS orbital energy. Within the direct RPA, the interaction kernel

*v*

_{pq,rs}= (

*pq*|

*sr*) is just the bare Coulomb interaction.

*GW*approximation. As in our previous work,

^{36,63}the

*GW*self-energy

^{64}

*ɛ*

_{q}of the spin-orbital (or spinor)

*ϕ*

_{q}and the Fermi level

*ɛ*

^{F}. 0

^{+}describes an infinitesimally small, positive number. Note that the self-energy as given in Eq. (7) contains an exchange part as well as a correlation part. Using Eq. (4), and subtracting the Coulomb potential

*v*from the screened interaction

*W*, the correlation part of the self-energy can be rewritten as

*η*= 2

*δ*, and where the fluctuation potential

*V*

_{m}(

**x**) is determined for all excited states

*m*as

*pq*|

*ρ*

_{m}) is computed as

*GW*correlation self-energies are tedious to compute, needing to fully solve Eq. (4). This severely limits the accessible molecular systems for

*GW*.

### C. Introducing natural virtual orbitals in the *GW* method

*GW*quasiparticle energy can be calculated as

*ω*. At the zeroth iteration

*ω*corresponds to the KS orbital energy

*ɛ*

_{p}. Σ

^{X}is the Hartree–Fock exchange energy of the p-th orbital, and $VpKS$ the KS potential of the p-th orbital.

*Z*is a linearization factor used to circumvent self-consistent solutions of Eq. (16), which would otherwise be necessary already for

*G*

_{0}

*W*

_{0}.

^{17}To reintroduce

*G*and

*W*explicitly, we rewrite Eq. (16) as

Depending on where NVOs are used for the calculation of the correlation self-energy, three different variants of the NVO-*GW* method can be derived as outlined in Eq. (17). In the first variant, we restrict ourselves by only representing the screened interaction *W* using NVOs. This approach is equivalent to using NVOs exclusively for the charge fluctuation *ρ*_{m}(**x**) by substituting the virtual orbitals (or spinors) *ϕ*_{c}(**x**) by their NVO counterparts $\varphi \u0303c(x)$ in Eq. (12), as well as in the eigenvalue problem in Eq. (4). This transformation of *W* is crucial to all approaches if NVOs are to be used since the solution of the eigenvalue problem in Eq. (4) is associated with the greatest computational effort. Note that in Eq. (13), in addition to the charge fluctuation *ρ*_{m}(**x**), also the energy eigenvalues *ω*_{m} in $Dp,q,m\xb1$ must be replaced by the eigenvalues obtained from Eq. (4) using NVOs. If besides the screened interaction *W* also Green’s function *G* is represented by NVOs, marked as variant 2 in Eq. (17), the whole correlation self-energy Σ_{c} is computed in the NVO space. In this case, also the sum over all virtual orbitals *ϕ*_{c}(**x**) in the second term on the right-hand side of Eq. (13) is replaced by NVOs as well as the corresponding orbital energies *ɛ*_{c} in $Dp,c,m\u2212$. In the third variant, in addition to the self-energy, also the orbitals or spinors $\varphi px$, for which the self-energy shall be computed, are transformed into a basis of natural virtual orbitals. This variant is relevant only if the considered orbital $\varphi px$ in Eq. (13) denotes an unoccupied orbital *ϕ*_{c}(**x**) because, otherwise, it is equivalent with the second variant. While the first and second variants require the parallel use of two basis sets during the *GW* calculation, the third approach considerably simplifies the technicalities of implementing the approach. In the latter, the virtual canonical orbitals are simply replaced by the NVOs throughout all necessary steps. However, the eigenvalues of the NVOs are considerably worse approximations of the quasiparticle energies, and therefore we mainly opt to use NVOs to evaluate Eqs. (4) and (12), that is, variant 1. The remaining terms in the correlation self-energy given in Eq. (13) are obtained using the canonical KS orbitals. Variant 1 is generally also compatible with self-consistent variants as for example, ev*GW*_{0}, during which *G* is updated but *W* is kept at its initial state. Introducing self-consistency in *W* through full eigenvalue self-consistent *GW* (ev*GW*) is however troublesome using NVOs. Starting from the second iteration, the eigenvalues used in the construction of *W* are no longer the NVO eigenvalues, but the ev*GW* ones, which are, in general, not compatible with the NVO eigenfunctions. This leads to significant errors in the unoccupied quasiparticle states. Due to the iterative nature, the errors in the unoccupied quasiparticle states subsequently also disturb the occupied quasiparticle states. Quasiparticle self-consistent *GW* (qs*GW*)^{65} is currently beyond the scope of this work, as it would require to also update the NVOs during the iterative process. We note in passing that, contrary to ev*GW*, the use of NVO for qs*GW* is nevertheless well-defined, as long as they are updated in each iteration. From a methodical scaling point of view, classical *GW* variants scale as *N*^{6} with system size. Therefore, using NVOs pays off given that they can be determined with at most *N*^{5} cost. For low-scaling *GW* algorithms, scaling as *N*^{4} or better, NVO is less rewarding. Nevertheless, the possibility of applying low-scaling *GW* algorithms is mostly limited to direct RPA response, while NVOs in combination with Eqs. (4) and (13) are generally applicable. We further note that the construction of the (d)MP2 density matrix can be made local, reducing its effective scaling drastically.^{66}

### D. Vertex corrections

Vertex corrections are the next step on the *GW* hierarchy of methods, surpassing the simple approach of using solely Coulomb interactions from the direct random-phase approximation.^{67–70} Unfortunately, calculating vertex corrections is computationally demanding, yielding rather complicated expressions. While it is yet to be determined if vertex corrections improve VIEs obtained from the *GW* method,^{71} current investigations of transition-metal monoxides have revealed some importance for certain molecular systems.^{72} The high costs of calculating vertex corrections, however, severely hamper further investigations, being unfeasible for molecular systems of chemically relevant size. As NVOs provide a convenient route to reducing the cost of obtaining vertex corrections, we will expand our discussion by including NVOs in their computation.

*et al.*use ordered second-order diagrams to obtain the vertex correction Δ

*ɛ*

_{p}(

*ω*) as

*W*in Eq. (18) refers to the full static screened interaction as described in Eq. (9), where

*ω*= 0 is set. The latter choice is termed test-charge–test-charge (TC–TC) interaction, and labeled

*W*

^{TC–TC}in Ref. 68. Evaluating Eq. (18) scales as

*N*

^{6}, providing a steep increase in computational demands for

*GW*. Using the resolution-of-the-identity approximation for the construction of

*W*

^{TC–TC}can reduce this to

*N*

^{5},

^{22}and accordingly using NVOs can further accelerate the construction of vertex corrections. In our current implementation, we use NVOs for all orbital functionals as well as in the construction of

*W*

^{TC–TC}, yielding an efficient way to calculate vertex corrections.

## III. COMPUTATIONAL DETAILS

The NVO generation scheme outlined in Table I has been implemented into a locally modified version of TURBOMOLE 7.7, making use of OpenMP parallelization.^{73,74} The corresponding necessary modifications of the *GW* code of TURBOMOLE have been carried out in the course of this work. To test the performance of NVOs within the *GW* method, VIEs are obtained with NVOs using the linearized *G*_{0}*W*_{0} and the ev*GW*_{0} approaches. As reference VIE values, the canonical results (obtained using canonical orbitals), will be used. This allows us to assess the deviations caused directly by truncating the NVO space. Furthermore, we compare the obtained VIE values, with and without using truncating the NVO space, to accurate CCSD(T) reference values. *GW*-based VIEs will be evaluated for three density-functional approximations (DFAs), namely, the global hybrid PBE0,^{75} the range-separated hybrid *w*B97X-D,^{76} and the local hybrid TMHF.^{77} Hybrid functionals, especially range-separated and local hybrid functionals, have earlier shown to yield excellent results for VIEs in combination with the *G*_{0}*W*_{0} method.^{42,43} This will allow one to conclude the transferability of NVOs between different reference functionals. All DFT calculations were tightly converged, with energy and density-matrix thresholds being set to 10^{−9} *E*_{h} and 10^{−8}, respectively. Fine integration grids of size 4 or better were used to obtain fully converged results for all tested molecules and functionals.^{78} To determine the dependence of the NVO approximation on a given basis-set expansion, all *G*_{0}*W*_{0} and ev*GW*_{0} computations were carried out using the triple-*ζ* def2-TZVPP basis set.^{79} Additionally, for an initial investigation, the *G*_{0}*W*_{0} computations were also performed using the quadruple-*ζ* def2-QZVPP basis sets.^{79} All DFT, NVO, and *GW* calculations were performed using a locally modified development version of TURBOMOLE 7.7.^{73} The resolution-of-the-identity approximation has been used in all NVO and *GW* calculations, and the corresponding def2-TZVPP and def2-QZVPP auxiliary basis sets optimized for post-HF calculations (denoted as cbas in TURBOMOLE jargon), were used.^{80} See Ref. 80 and references therein for a detailed description of the auxiliary basis sets.

For the GW100 test set,^{81} we compare the computed VIEs obtained as *GW* quasiparticle energies with reference data obtained at the coupled-cluster level, including single and double substitutions as well as a perturbative correction for connected triple substitutions, that is, at the CCSD(T) level. The CCSD(T) reference values have been carefully re-assessed, and changes and new values are outlined in Sec. IV. All EOMIP-CCSD and EOMIP-CCSDT calculations have been performed using the CFOUR program.^{82} All of our primary reference data [geometries, UHF/def2-TZVPP, UHF/def2-QZVPP, CCSD(T)/def2-TZVPP, and CCSD(T)/def2-QZVPP energies, number of orbitals excluded from the correlation treatment] are published as supplementary data in the supplementary material.

## IV. COUPLED-CLUSTER REFERENCE DATA

The CCSD(T)/def2-TZVPP reference data have already been published in an earlier work,^{83} and in the present paper, we specifically report revised data for four molecules of the test set. Furthermore, CCSD(T)/def2-QZVPP reference data for the GW100 test set have been computed, and a full overview of the calculated VIEs is provided in the supplementary material. In our previous work,^{83} the CCSD(T) computations were based on unrestricted Hartree–Fock (UHF) reference determinants, both for the cationic and neutral systems, and for all systems, we had attempted to find the UHF solution with the lowest energy. Recently, Bruneval, Dattani, and van Setten (BDvS) have repeated our CCSD(T)/def2-TZVPP calculations using a restricted Hartree–Fock (RHF) reference determinant for the neutral systems instead.^{84} In the course of the latter work, BDvS found UHF solutions with lower energy for the cations of uracil (CAS Reg. No. 66-22.8), cytosine (71-30-7), and hydrogen cyanide (74-90-8). We have adopted these solutions in our present work. Furthermore, BDvS pointed out in their work that our previously reported value for carbonyl selenide (1603-84-5, OCSe) had erroneously been computed for carbonyl sulfide (in the OCSe geometry). In the present work, we report the correct values for the carbonyl selenide molecule.

After the above revisions, only two clear outliers remain when comparing our CCSD(T)/def2-TZVPP data with the data of BDvS. These outliers are sulfur dioxide (7446-09-5, SO_{2}) and magnesium oxide (1309-48-4, MgO), which will be discussed below. Excluding SO_{2} and MgO from the comparison, the mean deviation of our values from those of BDvS amounts to 0.003 eV. The mean absolute deviation is 0.02 eV, and the sample standard deviation is 0.05 eV (note that data for the cesium dimer are missing from DBvS’s supplementary data, such that the comparison is made for 97 molecules).

The problem with SO_{2} is that we use the lowest-energy UHF solution for the $SO2+$ cation, which amounts to −546.881 967 *E*_{h}, whereas BDvS use the solution with *E*_{UHF} = −546.861 914 *E*_{h}. Admittedly, our lowest-energy solution is heavily spin-contaminated $(\u27e8S\u03022\u27e9=1.13\u210f2)$ and the corresponding CCSD(T) energy does not seem accurate. The CCSD(T) energy obtained with BDvS’s UHF solution is much closer to the CCSD(T) energy obtained with a restricted open-shell Hartree–Fock (ROHF) reference determinant for $SO2+$ (Table II). Also the EOMIP-CCSD (equation-of-motion-ionization-potential coupled-cluster-singles-and-doubles) value^{85} favors BDvS’s CCSD(T) result.

Type of computation . | Hartree–Fock . | CCSD . | CCSD(T) . | VIE (eV) . |
---|---|---|---|---|

UHF (this work) | −546.881 967 | −547.461 26 | −547.488 62 | 13.486 |

UHF (BDvS)^{84} | −546.861 914 | −547.494 98 | −547.532 27 | 12.298 |

ROHF | −546.844 345 | −547.493 36 | −547.532 89 | 12.281 |

EOMIP-CCSD | −547.496 39 | 12.373 | ||

EOMIP-CCSDT | 12.225 |

Type of computation . | Hartree–Fock . | CCSD . | CCSD(T) . | VIE (eV) . |
---|---|---|---|---|

UHF (this work) | −546.881 967 | −547.461 26 | −547.488 62 | 13.486 |

UHF (BDvS)^{84} | −546.861 914 | −547.494 98 | −547.532 27 | 12.298 |

ROHF | −546.844 345 | −547.493 36 | −547.532 89 | 12.281 |

EOMIP-CCSD | −547.496 39 | 12.373 | ||

EOMIP-CCSDT | 12.225 |

The problem with MgO is that we use the lowest-energy UHF solution for the neutral molecule, whereas BDvS use the closed-shell restricted Hartree–Fock (RHF) solution (Table III). Also in this case, the EOMIP-CCSD value of 8.172 eV for the VIE favors BDvS’s CCSD(T) result (7.909 eV). We acknowledge that our value (7.487 eV) is very likely to be too small by about 0.4 eV and suggest using the EOMIP-CCSDT (EOMIP-CCSD-and-triples) values for the molecules SO_{2} and MgO instead. Those VIE values amount to 12.225 and 7.877 eV for SO_{2} and MgO, respectively, as computed by the CFOUR program.^{82} In the def2-QZVPP basis set, the corresponding EOMIP-CCSDT values are 12.403 and 7.939 eV, respectively (cf supplementary material).

## V. RESULTS

### A. Ionization energies using NVOs

To test the capabilities of NVOs within the *GW* approximation, as a first test, we calculate quasiparticle energies for the GW100 test set. The latter consists of VIEs of 100 small molecules.

Type of computation . | Hartree–Fock . | CCSD . | CCSD(T) . | VIE (eV) . |
---|---|---|---|---|

UHF (this work) | −274.464 527 | −274.813 34 | −274.820 75 | 7.487 |

RHF (BDvS)^{84} | −274.374 647 | −274.811 94 | −274.836 26 | 7.909 |

EOMIP-CCSD | 8.172 | |||

EOMIP-CCSDT | 7.877 |

Type of computation . | Hartree–Fock . | CCSD . | CCSD(T) . | VIE (eV) . |
---|---|---|---|---|

UHF (this work) | −274.464 527 | −274.813 34 | −274.820 75 | 7.487 |

RHF (BDvS)^{84} | −274.374 647 | −274.811 94 | −274.836 26 | 7.909 |

EOMIP-CCSD | 8.172 | |||

EOMIP-CCSDT | 7.877 |

Figure 1 outlines the general convergence of the correlation self-energy as NVOs are subsequently added, starting from 20% to 100%, with the latter value exactly reproducing the untruncated canonical value. First, it is visible that simply truncating the space spanned by the canonical orbitals in *GW* calculations is not a viable choice. The convergence toward the canonical result is slow, with any practically useable truncation parameter leading to unacceptable errors. Even including 80%–90% of the canonical orbitals does not yield converged results. Contrary, very good results are obtained if NVOs are used to calculate the self-energy Σ. In Fig. 1, NVOs are employed to evaluate Eqs. (4) and (12), consistent with variant 1 outlined in Sec. II A. Converged results are obtained if 30%–40% of the total number of virtual orbitals are included. Before this point, the error curve is less steep and the overall deviation therefore already largely reduced. This effect gets more prominent as more virtual orbitals are discarded. Note that the basis set size does not seem to be too important. The def2-QZVPP basis set shows deviations similar to those of the def2-TZVPP basis set, with the latter showing slightly reduced deviations in the case of a very low number of virtual orbitals taken into account.

Next, we will analyze the deviation introduced by NVOs dependent on the system size and the cardinality of the correlation-consistent basis set used.^{86,87} Therefore, the deviations of the VIEs in the acene series from benzene to hexacenes are shown in Fig. 2 for a 50% truncated virtual space. Furthermore, Fig. 2 also shows the deviations obtained for benzene with increasing basis set size.

Figure 2 exhibits that the deviation of the VIE normalized to Σ_{C} obtained in the fully intact virtual space decreases with the system size. The latter effect is to be expected as the virtual expansion size linearly increases in the row, providing more flexible combinations from which the NVOs can be generated. Moreover, Fig. 2 outlines how the deviation of the VIEs converges with an increasing number of basis functions for the benzene molecule. Again, raising the cardinal number of the basis set (adding a considerable amount of basis functions) provides a more flexible basis from which the NVOs can be generated. Discarding half of the virtual orbitals in each case therefore still leads to an increasingly efficient description of the virtual space by the NVOs, significantly increasing the accuracy of NVO calculations. Overall, Fig. 2 points out at large molecules combined with large basis sets being a primary target for NVOs, with the latter working exceptionally well in these cases.

Table IV outlines the timings needed to build NVOs from scratch for the acenes shown in Fig. 2. Even on moderate computational resources, building NVOs for systems with up to 1000–1500 Cartesian basis functions (BFs) only takes seconds to a few minutes. Compared to classical *N*^{6} scaling *GW* calculations, this is a vanishing effort. Still, for extensive systems, the overall *N*^{4}–*N*^{5} effort needed to construct canonical NVOs will be significant once 2000–3000 BF are surpassed. In this case, local methods to construct the (d)MP2 density matrix need to be employed to reduce the overall scaling in the NVO construction.

Acene . | N_{BF}
. | T_{MP2} (s)
. |
---|---|---|

C_{6}H_{6} | 306 | 1.3 |

C_{10}H_{8} | 480 | 11.3 |

C_{14}H_{10} | 654 | 50.9 |

C_{18}H_{12} | 828 | 161 |

C_{22}H_{14} | 1002 | 403 |

C_{26}H_{16} | 1176 | 910 |

Acene . | N_{BF}
. | T_{MP2} (s)
. |
---|---|---|

C_{6}H_{6} | 306 | 1.3 |

C_{10}H_{8} | 480 | 11.3 |

C_{14}H_{10} | 654 | 50.9 |

C_{18}H_{12} | 828 | 161 |

C_{22}H_{14} | 1002 | 403 |

C_{26}H_{16} | 1176 | 910 |

### B. Ionization energies using NVOs for the GW100 test set

To analyze the overall deviation stemming from truncating the virtual orbitals space, we determine the error in the quasiparticle energies for the GW100 test set using various settings that determine how many virtual orbitals are neglected (frozen). Figures 3 and 4 outline a more detailed statistical analysis of the usage of NVOs in the GW100 test set for the def2-TZVPP and def2-QZVPP basis sets, respectively. Both plots show the obtained deviations for a 50% reduced virtual space and for three different occupation number thresholds that have been used to choose NVOs of lesser importance. Box plots were used to give a compact representation of the data distribution with respect to the full virtual orbital space. The horizontal orange line indicates the median and the upper and lower limits of the box (upper and lower quartiles) indicate that 25% of the data values are located above and below the quartile. The range of the box is the interquartile range (IQR). The whiskers extend to the data with values less than 1.5 times the interquartile range plus the lower or upper quartile.

First, comparing Figs. 3 and 4 suggests that a required occupation number threshold is dependent on the basis set size. For example, for the def2-TZVPP basis set, a threshold of less than 10^{−4} below which NVOs are discarded leads to basically converged results. Contrary, the def2-QZVPP basis set requires a threshold of less than 5 × 10^{−5} to reduce the deviations found for the VIEs to a similar accuracy. For comparison, Figs. 3 and 4 also outline results with truncated canonical orbital spaces, where the number of active virtual orbitals has been chosen to match those of the NVOs. Figures 3 and 4 confirm the initial assumption that truncating canonical orbitals is a rather error-prone approach. For any truncation threshold, truncating canonical virtual orbitals yields a broad error distribution, being unsuitable to obtain quasiparticle energies. While truncating the NVO space is superior to truncating canonical virtual orbitals, both plots exhibit that the exact means of the generation of the NVOs are not too important. Quasiparticle energies obtained from NVOs generated from a direct-MP2 density matrix (dNVOs) are practically indistinguishable from their standard counterparts. Especially for larger molecules, where high-performance algorithms for Laplace-transformed direct MP2 have been developed,^{58,59} this can be translated in considerable gains while practically no accuracy is lost. Concerning the choice of the density functional approximation, no dependence of the NVOs on the underlying functional can be observed. Figure 5 depicts the fraction of the discarded virtual orbitals in dependence on the occupation number thresholds used in Figs. 3 and 4.

As outlined by Fig. 5, the tighter truncation thresholds for def2-QZVPP still lead to a higher relative number of NVOs that are discarded, translating into more pronounced speedups. To further check for the validity of NVOs in *GW* calculation, we also investigate the behavior in the ev*GW*_{0} method. In this method, specifically, the screened exchange *W* is kept at the KS level. Accordingly, *W* can be calculated using the NVOs obtained from the KS reference, and we update only the Greens function in every iteration.^{88} Deviations found for the ev*GW*_{0} method, shown in Fig. 6, are in principle identical to the previously observed deviations for the *G*_{0}*W*_{0} approximation. The same conclusions drawn before for the truncated virtual spaces and difference density functional approximations therefore also apply to ev*GW*_{0}. Using (d)NVOs can therefore help to considerably accelerate ev*GW*_{0}. Notably, Fig. 6 shows that there is a single outlier for the smaller thresholds. This outlier is in every case helium, which reacts very sensitively to the omission of any orbitals in the *GW* treatment. For atomic cases, especially for helium, NVOs should therefore not be used.

Figure 7 depicts the statistical errors obtained for the VIE of the energetically lowest occupied orbital. This orbital corresponds to the 1s orbital when no effective core potential (ECP) is used, and to the first orbital not included in the ECP otherwise. For these vertical core ionization energies, a generally amplified magnitude of errors is observed. The latter is caused by the significantly raised ionization energies needed to expel core or non-valence electrons. Non-valence electron VIEs exhibit an increased sensitivity to the chosen occupation number threshold, with only the tighter thresholds yielding acceptable results. Notably, using canonical orbitals, the VIEs of these orbitals are strongly overestimated. In contrast, as outlined in Figs. 4–6, valence electron VIEs are strongly underestimated with truncated canonical orbitals. The exceedingly large error of VIEs when frozen canonical orbitals are used, reaching up to 10–20 eV, basically renders this approach useless for core ionization energies. Using NVOs, these errors can be largely reduced. The overall accuracy with frozen NVOs is then again comparable to the overall accuracy that can be reached with *GW* methods while speeding up the calculation by an order of magnitude or more. Notably, the outliers in Fig. 7 can be grouped into two classes. The first class is constituted by elements where ECPs have been employed, as e.g., I_{2}. The second class is built by elements with 1s VIEs exceeding 10 000 eV, for example, Br_{2} or GeH_{4}. For orbitals close to an effective core potential, using truncated NVOs is therefore not recommendable. For the second class, while larger than usual absolute errors are obtained, the relative error, which is also connected to the overall accuracy of *GW*, is very good. The presented approach can therefore also be recommended for core ionization energies.

Finally, the question remains of how much the VIEs are impacted at all by the use of a truncated set of virtual orbitals. Figure 8 compares the VIEs obtained from CCSD(T)/def2-TZVPP, *G*_{0}*W*_{0}@DFT/def2-TZVPP, and *G*_{0}*W*_{0}@DFT/def2-QZVPP to the reference VIEs obtained from CCSD(T)/def2-QZVPP.

In Fig. 8, it is visible that NVOs are not lowering the overall accuracy of the *G*_{0}*W*_{0} method when compared to the results using the full virtual space. The overall error when compared to CCSD(T)/def2-QZVPP reference values is hardly affected at all, with the NVO/def2-TZVPP combination even exhibiting a slightly improved performance compared to a full virtual space. The latter effect can however safely be attributed to error cancellation. It is therefore important to point out two things. First, using NVOs one could accelerate a calculation using the same basis set, without a significant loss in accuracy. Second, using a truncated NVO space, quadruple-*ζ* results can be obtained basically at triple-*ζ* cost or less. As shown in Fig. 8, the accuracy improvement of going from triple-*ζ* to quadruple-*ζ* far outperforms any error introduced by the NVOs. Contrary, truncating the canonical orbital space significantly deteriorates the obtained VIEs. Figure 8 estimates that the average errors are 2–3 times larger in the latter case. Note that even CCSD(T), when used with a smaller basis set, is considerably in error. Deviations of up to 0.25 eV are found for the latter basis set when compared to the reference values, with a mean error of ∼0.1 eV.

### C. Vertex corrections

Subsequently, to standard *G*_{0}*W*_{0} and ev*GW*_{0} calculations, NVOs are especially interesting to speed up the calculation of otherwise rather expensive higher-order vertex corrections to the *G*_{0}*W*_{0} or ev*GW*_{0} quasiparticle energies. Using NVOs, the virtual space can be truncated, which leads to large gains when evaluating Eq. (18). Especially the evaluation of the first term in Eq. (18), featuring a sum over two virtual indices, will benefit from truncating the virtual space. Of course, the question of how much NVOs will affect the results remains. It turns out that for vertex corrections the deviation introduced by using NVOs is again reduced quickly and smoothly with the occupation number threshold as outlined in Fig. 9.

Given that the magnitude of the vertex corrections Δ*ɛ*_{p} is usually much smaller than the initial correlation self-energy obtained from *G*_{0}*W*_{0} or ev*GW*_{0}, smaller deviations are obtained overall. Still, truncating the virtual space while using canonical orbitals leads to significant deviations. In the latter case, the incorporated error is of the same order of magnitude as the vertex corrections itself, yielding unacceptable deviations. Contrary, truncating the virtual space is possible without sacrificing much accuracy when NVOs are used. With NVOs the error is well under control, especially for lower occupation numbers. An occupation number threshold of 10^{−4} again seems to provide a reasonable balance between performance and accuracy, with the mean error on the vertex correction approaching 0.01–0.02 eV. The error distribution for vertex corrections is therefore in line with that of the quasiparticle energies themselves. Using the same thresholds for both, vertex corrections and quasiparticle energies, is therefore possible and advisable, yielding a consistent yet efficient and accurate method. As a final remark, it is interesting to note that for the lowest occupation number thresholds of 5 × 10^{−5} and 10^{−4}, the single largest outlier exhibiting double the error of the second largest one is again helium. The latter already showed to be troublesome for ev*GW*_{0} and is, in general, known to be a challenging case for *GW* methods.^{89} We, however, note that this is only a minor setback, as truncating orbital spaces for a single atom is hardly ever necessary.

## VI. CONCLUSION

We have outlined the principal capabilities of using the truncated natural virtual orbital spaces to accelerate the evaluation of quasiparticle energies from the *GW* approximation. For VIEs, we find that the NVO space can be considerably truncated without a significant loss in accuracy in the quasiparticle energies obtained from either *G*_{0}*W*_{0} or ev*GW*_{0}. The truncation thresholds used to construct a set of NVOs are furthermore shown to be transferable between different classes of density functional approximations. However, it is not automatically transferable between basis sets, and most likely an optimal value needs to be determined for at least every cardinal number of orbital basis sets. We recommend truncating the space spanned by the NVOs by neglecting orbitals with occupation numbers of less than 10^{−4} for triple-*ζ* basis sets, and a slightly tighter threshold of 5 × 10^{−5} for quadruple-*ζ* basis sets. Using these thresholds leads to tightly converged quasiparticle energies, and also to tightly converged vertex corrections. The remaining error of using NVOs has then been shown to be insignificant compared to the overall error of the *GW* method. This also holds for non-valence orbitals, for example, 1s orbitals of lighter elements. For the latter, NVOs yield excellent quasiparticle energies, while largely reducing the cost of the associated *GW* calculation. It is furthermore interesting to note that the exact way the correlated NVOs are obtained is not too important. NVOs obtained from direct MP2 lead to practically indistinguishable result, while speeding up the process by an order of magnitude. For periodic *GW* calculations, where MP2 density matrices are not commonly available, one could resort to RPA density matrices, which have similar properties and are expected to work as well as dMP2 density matrices for NVO generation due to their principal similarity. Given the robustness of the NVO scheme, exhibiting no significant difference between MP2 and dMP2 density matrices, we also expect RPA density matrices to provide results of similar quality. These results are especially encouraging for vertex corrections, which are generally considered to be computationally demanding. The latter can be accelerated using NVOs to yield methods that are routinely applicable to molecular systems, allowing future investigations to take them into account when needed.

For electron attachment energies (not reported here), we find NVOs to be not suitable within the *GW* ansatz when the full virtual space is replaced by NVOs. The failure can be attributed to the fact that the virtual eigenvalues of the Fock matrix in the truncated space of the NVOs should not be used as starting points for the computation of the respective quasiparticle energies. The problem can be remedied by not replacing those virtual orbitals that are relevant for the electron attachment energies by NVOs. This must also be taken into account when computing excitation energies through the Bethe–Salpeter equation, in which quasiparticle energies of occupied and virtual levels are used as input.

## SUPPLEMENTARY MATERIAL

See the supplementary material for the geometries used in the GW100 test set, the Hartree–Fock and CCSD(T) energies in the def2-TZVPP and def2-QZVPP basis sets for the cations as well as the neutral systems, and energies for the molecules SO2, MgO, and their cations. Hartree–Fock and CCSD(T) energies in the def2-TZVPP and def2-QZVPP are additionally available as machine-readable ASCII files. Furthermore, the G0W0 and evGW0 results are also reported as machine-readable ASCII files.

## ACKNOWLEDGMENTS

L.M. and W.K. gratefully acknowledge financial support by the Deutsche Forschungsgemeinschaft (DFG) through the Research Training Group 2450 on “Tailored Scale-Bridging Approaches to Computational Nanoscience.” C.H. gratefully acknowledges the Volkswagen Stiftung for financial support.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Laurenz Monzel**: Data curation (lead); Formal analysis (lead); Software (lead); Visualization (lead); Writing – original draft (lead). **Christof Holzer**: Data curation (lead); Formal analysis (lead); Software (lead); Visualization (lead); Writing – original draft (lead). **Wim Klopper**: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Project administration (equal); Supervision (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available within the article and its supplementary material.