On the basis set selection for molecular core-level $GW$ calculations

The $GW$ approximation has been recently gaining popularity among the method for simulating molecular core-level X-ray photoemission spectra. Traditionally, $GW$ core-level binding energies have been computed using either the cc-pV$n$Z or def2-$n$ZVP basis set families, extrapolating the obtained results to the complete basis set limit, followed by a an element-specific relativistic correction. Despite of achieving good accuracy, these binding energies are chronically underestimated. By using first-row elements and standard techniques known to offer good cost-accuracy ratio in other theories, we show that the cc-pV$n$Z and def2-$n$ZVP families show large contraction errors and lead to unreliable complete basis set extrapolations. On the other hand, we demonstrate that uncontracted versions of these basis sets offer vastly improved convergence. Even faster convergence can be obtained using core-rich, property-optimized, basis sets families like pcSseg-$n$, pcJ-$n$ and ccX-$n$Z. Finally, we also show that the improvement over the core properties does not degrade the calculation of the valence excitations, and thus offer a balanced description of both core and valence regions.

The GW approximation (GWA) 5 to the self-energy is a Green's-function-based method that can be used to obtain accurate molecular CLBEs at a reasonable cost. [6][7][8][9] The CLBEs obtained with this approach include some level of orbital relaxation effects, however, the quality of the results might depend strongly on the starting point. For example, it has been known that a large fraction of exact exchange is needed at the one-shot G 0 W 0 level, mainly due to a large self-interaction error (SIE) present when standard "pure" exchange-correlation density functional approximations (DFAs) are adopted. 6,7,9 A possible workaround when starting from a pure Kohn-Sham result is to use the more demanding evGW partial self-consistent approach, 7-10 although some starting-point dependency will remain.
Besides these unique requirements, CLBEs obtained with the GWA will also depend on the quality of the basis set. The accuracy of GWA-based CLBEs has been studied recently by means of two benchmark datasets. 6,8,10 Both studies extrapolated the CLBEs to the complete basis-set (CBS) limit in order to eliminate the dependency on the basis set. In Ref. 6, the authors used the def2-TZVP and def2-QZVP 11 basis sets and omitted relativistic corrections; while in Refs. 7 and 8, the authors used Dunning's correlation-consistent basis sets cc-pVnZ, 12 from triple-to sextuplezeta quality, and included relativistic corrections. Both studies showed that GW can be used to predict CLBEs with a mean absolute error of about 0.3 eV for first row elements, but that this error can only be "inconsistently" reduced by adding relativistic corrections. 10 Both benchmark studies also relied on the extrapolation of the CLBEs, obtained with standard basis sets, to the CBS limit. Standard CBS extrapolation techniques were developed to extrapolate the correlation energy of small molecules. [13][14][15] This is usually a small-energy regime with wellknown basis set convergence behavior, 16,17 and asymptotic behavior-E c /Z ∼ −0.02073 ln Z + 0.0372 18,19 -for atoms. In contrast, the energy involved for ionizing the core is rather largeε 1s ∼ Z 2 /2-and does not follow the same convergence trend. This may explain why Golze and co-workers used a a linear regression fit with respect to the inverse number of basis functions 8 compared with Bruneval and co-workers, who used the more common exponential convergence formula for GW valence binding energies. 20,21 Even when the correlation energy of molecules and the core ionization energy have similar asymptotic behaviors, the basis set requirements for a response method, like the GWA, might be very different compared with a single-point energy calculation. It has been shown that basis sets optimized for energies often fail to describe other molecular properties due to a poor description of the orbitals in regions important for the given property but otherwise energetically unimportant. 22,23 This numerical issue prevents the application of systematic improvements and of reliable extrapolation approaches. A natural consequence is that extrapolation using non-optimal basis sets might result in values far away from the real CBS limit.
The importance on the selection of the basis set have already been addressed for other response methods, like the equation-of-motion coupled-cluster 23,24 or the linear-response time-dependent density functional theory formalism. 25 Basis set convergence has also been studied within the GWA applied to molecules, mostly focused to valence ionization potentials, 20,21,26 usually showing smooth, albeit slow, convergence. However, a systematic study comparing CLBEs obtained with basis sets other than the cc-pVnZ and def2-nZVP families is not available in the literature.
In this study, we use the CORE65 dataset, comprised of 65 CLBEs of first-row elements, in order to test the effect of the basis set core-level GW calculations. We asses two types of basis sets that are known to offer very good results in other theories: uncontracted versions of traditional energy-optimized basis sets, and core-rich property-optimized basis sets.
In turn, we show that both cc-pVnZ and def2 basis set families are not well suited to describe first-row 1s CLBEs within the GWA, often resulting in very slow convergence and unreliable CBS extrapolations. Findings of this work are consistent with earlier studies dealing with other linearresponse methods, by showing that while uncontracted versions of these basis sets offer vastly improved results, the CLBEs seem to converge to a different value as compared to the CLBEs obtained with other core-rich basis sets like pcJ-n, 27 pcSseg-n, 22 and the very recent ccX-nZ 23 families. The difference is small (around 100 meV), however, this value is close to the intrinsic accuracy of the GWA for CLBEs and, as such, it is worth noticing.

Overview of the GW Approximation
The central object of the GWA is the one-particle Green's function G describing particle and hole scattering in the interacting many-body system. In order to obtain the electron addition and removal energies from such a Green's function, a non-local and dynamic effective potential, the self-energy Σ, is often introduced. The self-energy Σ substitutes the mean-field exchange-correlation operator, and the GW quasiparticle (QP) energies ε GW n can be obtained as corrections to the mean-field energies ε n : Here, σ denotes the spin index, and V xc nσ and Σ nσ denote the nth diagonal element of the corresponding matrix representation in the molecular orbital basis. Equation 1 is non-linear and must be solved iteratively.
The self-energy operator Σ is given in terms of the Green's function G σ and the screened Coulomb interaction W : with η being a positive infinitesimal. In practice, the GWA is often performed as a one-shot perturbative approach known as G 0 W 0 . In this case, G σ 0 is the non-interacting mean-field Green's function, and W 0 is obtained using the random phase approximation (RPA) as In the preceding equations, µ is the Fermi-level of the system, v(r, r ) is the bare Coulomb inter-action, and (r, r , ω) is the RPA dynamical dielectric function: In the RPA, the irreducible polarizability χ 0 has a simple sum-over-states representation 28,29 where the index i runs over the occupied orbitals while the index a runs over the virtual ones.
The contour deformation (CD) technique is used in this work for the accurate integration of Σ, as it has been shown to have a good cost-effective profile for the evaluation of CLBEs. Further details about CD-GW and its implementation using local orbitals can be obtained from References 7, 30, and 9 for example.

Computational Details
The adequacy of the cc-pVnZ and def2-nZVP basis set families for core-level GW calculations was evaluated by computing the 65 CLBEs from the CORE65 benchmark dataset. 8  In order to find other appropriate basis sets families, the 65 CLBEs were also obtained with the pcJ-n, pcSseg-n, and ccX-nZ families, as well as their uncontracted versions.
Since our CD-GW implementation fits the four-center ERIs, an adequate fitting basis is needed.
Here, we do not directly assess the impact of various fitting basis sets on the quality of CLBEs but we do use automatically generated fitting bases using the "AutoAux" procedure 35 for each orbital basis set.

The cc-pVnZ and def2-nZVP families
We first begin studying the two families used in the CLBEs benchmark studies in the literature, the cc-pVnZ and def2-nZVP families, respectively. Both basis set families were developed with total energies in mind. In energy-optimized basis sets, contraction of the core basis functions is usually performed in order to achieve high efficiency at the price of a small loss in accuracy. 36 The removal of a core electron, however, induces a significant relaxation of the core orbital incompatible with a Given that the flexibility around the core seems to play a very important role in determining the quality of the GW CLBEs, the use of core-rich basis sets might be a better choice than the use of either un-cc-pVnZ or un-def2-nZVP. 23 As mentioned earlier, we could not find a systematic study of GW CLBEs using these basis sets. The only data available was found in the Supplemental Material of reference, 10 where the CLBEs obtained with the core-augmented cc-pCVnZ basis set was compared with CLBEs obtained with the standard cc-pVnZ basis sets. In that study, the authors noted that the errors of the cc-pCVnZ were not systematic. A similar conclusion was reached in Reference 23, where the authors noted that uncontracting the cc-pCVnZ family leads to larger deviations with respect to the CBS limit, a result which suggests an uncontrolled error cancellation in the contraction scheme. Here, we have extended the results presented in Reference 10 and have included the pcJ-n, pcSseg-n, and ccX-nZ core-rich basis set families. Figure 4 shows CLBEs deviations with respect to experiment for all basis set families tested. The same color code as in the previous sections is being used: blue hues for contracted sets, and red hues for uncontracted sets. Note that relativistic corrections have not been included.
The three core-rich basis set families predict converged CLBEs in very good agreement to the ones obtained using the uncontracted def2-nZVP or cc-pVnZ basis sets. Moreover, there is no longer a mismatch between the CLBEs obtained with contracted and uncontracted versions of these basis sets. The fact that the same limit is achieved by all uncontracted basis sets, as well as with the contracted pcJ-n, pcSseg-n, and ccX-nZ, strongly supports the idea that this limit is in fact closer to the CBS than the extrapolations given in Reference 8.
Another aspect worth noting is the manner in which the CLBEs converge for each basis set. As noted before, the def2-nZVP basis sets show the wrong shift going from triple-to quadruple-ζ, the cc-pVnZ family converges very slowly from below, the pcSseg-n family converges from above, the pcJ-n family converges from below, and the ccX-nZ family is very well converged even with double-ζ quality.

Relative core-level binding energies
Chemical shifts show faster convergence with any of the basis set used in this study, as shown in Figure 5. The contracted pcJ-n and ccX-nZ tend to show smaller dispersions overall, however, this effect is rather small and barely noticeable. A special case is obtained with cc-pV5Z for Nitrogen 1s chemical shifts, as Figure 5 shows that these chemical shifts are too small as compared to the consensus of the remaining bases. We have not been able to identify the cause for this discrepancy, but note that the chemical shifts obtained with the uncontracted version are in very good agreement with the rest of the results.
Since the relative binding energies are less sensitive to the basis set chosen, the statistics presented by Golze et al. 8 does not change by more than 0.05 eV when compared to our ccX-QZ results.

Valence binding energies
Improvement in convergence of core-level binding energies should not, in principle, come at the expense of a poor description of the valence region. We evaluated the accuracy achieved with each one of the contracted basis set families by means of a small subset of the GW100 dataset. 37,38 The subset is made exclusively of molecules containing B, C, N, O, and/or F (the ccX family has only been generated for these elements). The result is a test set with 10 vertical ionization potentials (N 2 , F 2 , C 4 , CF 4 , BF, BN, CO 2 , CO, O 3 , C 6 F 6 ) and 6 vertical electron affinities corresponding to the experimentally-bound anions (F 2 , C 4 , BN, CO, O 3 , C 6 F 6 ). In order to directly compare to existing benchmark data, 38,39 the results presented in this section used the "pure" PBE functional at the G 0 W 0 level without basis set extrapolation. Table 1 shows the vertical ionization potentials (VIP), in eV, obtained with the different basis sets tested. It is immediately apparent that the core-rich basis sets converge as fast as the traditional cc-pVnZ and def2-nZVP families. In contrast with the CLBEs case, the VIPs computed with all basis set families tend to the same value.
Slightly larger variations in the rate of convergence for each basis set family can be seen in the vertical electron affinities (VEA) series shown in Table 2. Here, it is evident that the corerich basis sets offer better VEAs at the double-and triple-ζ level than the corresponding cc-pVnZ and def2-nzvp sets, especially for molecules containing fluorine atoms. The good performance of the core-rich basis sets for VEAs is partly explained by the presence of some more diffuse basis functions. Table 3 shows the smallest exponent of each angular momentum for all double-ζ basis sets. To facilitate, all exponents are given as the fraction of the corresponding def2-SVP ones. It is evident that the ccX-nZ family is consistently more diffuse than their counterparts, hence it good performance for GW VEAs calculations, but the story is not as evident for the other core-rich families. The remaining difference is explained by the additional variation flexibility core-rich basis set have in comparison to standard ones. For example, the uncontracted def2-SVP VEAs for F 2 and C 6 F 6 (-1.40 eV and -1.45 eV, respectively) show a slight improvement toward the basis set limit with respect to the standard def2-SVP. It is reassuring to see that the core-rich basis sets are very well balanced to describe both coreand valence-excitations at the GW level.

Summary
Convergence of the core-level binding energies obtained in GW calculations have often been deemed hard to achieve with respect to the basis set quality. Here, we have shown that the def2nzvp and cc-pVnZ basis sets families are shown to have substantial contraction errors when used to compute GW core-level binding energies of molecules containing first row atoms. Furthermore, we have shown that basis set extrapolations using such bases might be unreliable, as in the case of Exposing all the core basis functions by uncontracting either basis set adds enough flexibility to the calculation, shifting the binding energies to higher values and showing much faster convergence behavior. The same behavior can be observed with partially uncontracted basis sets, like pcJ-n, pcSseg-n, and ccX-n families. In particular, both pcJ-n and ccX-n seem very well suited to study GW core-level binding energies.
The quality of GW valence charged excitations does not degrade by using one of the core-rich basis sets. In fact, vertical electron affinities achieve faster convergence rates.
We therefore recommend the use of either pcJ-n or ccX-nZ families for core-level and valence GW calculations. However, it is important to note that both basis sets do impact the speed of the calculation. A more cost-effective choice might be the use of an uncontracted standard basis set, like cc-pVnZ and def2-nzvp. Note that standard J, JK, an RI fitting bases are not recommended when using uncontracted basis sets.
Another fact worth noting is that the current study focused only on the elements C to F (present in the CORE65 database); therefore the conclusions of this paper might not necessarily hold when applying GWA to molecules containing heavier elements.