Robust relativistic many-body Green’s function based approaches for assessing core ionized and excited states

The GW method has become a staple to assess the valence and core excited states of complex molecular systems, being a helpful theoretical tool for the well established fields of UV/Vis and x-ray spectroscopy. The latter is an especially valuable tool in the characterization of complex materials. The usage of high-energy radiation allows for the direct excitation of electrons occupying states close to the nuclei. Since the core energy states are highly element specific, x-ray spectroscopy can be employed to obtain information about the local (electronic) structure of a system. For example, it allows one to gain new insights into complex bonding situations in lanthanoid containing molecular systems.¹ Likewise, the description of these experiments using quantum mechanical methods is complicated and remains an active field of development, as the high-energy regime represents a difficult case, where small errors may be largely amplified.

While time-dependent density functional theory (TD-DFT) often provides a reasonable description, its applications to the x-ray regime regularly require the introduction of large energetic shifts,^2,3 which partly arise from the amplification of the self-interaction error in density functional approximations (DFAs).^4–8 Similarly, the ambiguity in the choice of the DFA needs to be studied carefully.^2,9,10 Contrary, wave function-based methods require a significant computational effort in order to systematically improve upon the more commonly employed DFT based approaches. Therefore, wavefunction based methods that are able to assess core excited states are limited to systems with only a few light atoms. Alternatives to TD-DFT and wavefunction based methods to assess core excited states with high accuracy are, therefore, needed to gain insights into the light–matter interactions of core states of more complex molecular systems. In a recent article, we combined the damped response formalism with the Bethe–Salpeter equation (BSE) or, more specifically, the composite GW-BSE approach. In the GW-BSE approach, first the quasiparticle energies are obtained from the GW method. In the second step, charge-neutral excited states are obtained from the BSE, which is solved in the static screened approximation.² While the initial study was limited to the description of small molecules and atoms at the one- and two-component levels, that is, including scalar relativistic effects and spin–orbit coupling (SOC), it could be shown that the GW-BSE method is a powerful tool to assess core excited states. A scalar-relativistic study carried out by Yao and co-workers later came to similar conclusions.¹¹ Both studies indicated that the GW-BSE method provides a reasonable alternative to other established techniques, especially with respect to experimental reference results, thereby largely eliminating the need for an artificial shift of the spectrum as needed in TD-DFT.^2,11,12 The computational scaling of the BSE is likewise the same as for TD-DFT, making it even more interesting in the framework of core excitations. Still, a straightforward application of the GW-BSE method to more realistic systems, that is, systems with more than a few atoms, provides additional challenges to the method. This is especially true if relativistic effects as, for example, spin–orbit coupling are taken into account.^13–15 Here, the calculation of the quasiparticle energies from the GW method may quickly become a bottleneck. On one hand, common techniques, such as the contour deformation (CD) approach introduced in the 1980s, scale unfavorably if low-lying states, such as core orbitals, are of concern.^16–18 On the other hand, the analytic continuation of the GW self-energy yields significantly deteriorated quasiparticle energies for core states.¹² 

The aim of the manuscript is twofold. First, we introduce a more pragmatic yet accurate approximation to the GW method for the calculation of the quasiparticle energies for core states in one- and two-component frameworks, based on the work of Duchemin and Blase.¹⁹ This work is then combined with reduced scaling GW algorithms,^12,19–28 with a special focus on current relativistic implementations.^14,15 Given the high interest in open-shell relativistic systems, such as many lanthanoid containing systems, we will furthermore expand our method toward them. The GKS determinants of systems with unpaired electrons usually do not exhibit Kramers symmetry and are, therefore, especially challenging for TD-DFT and the GW-BSE method.²⁹ Contour deformation based GW could, however, recently be extended to such systems, allowing us to step forward in this direction.¹⁵ The resulting strictly $O(N4)$ scaling relativistic frequency-sampled contour deformation GW (fsCD-GW) algorithm is able to efficiently assess quasiparticle states far away from the Fermi level, yielding core ionization energies with high accuracy. Furthermore, we hope that the fsCD-GW method will provide an entry point for even lower scaling GW algorithms.^30–35 fsCD-GW is finally combined with the BSE in the core–valence separation (CVS) approximation. The latter projects out, that is, freezes, valence orbitals, or spinors, to directly access core excited states, while including their screening effects in the dielectric matrix. Subsequently, we apply the resulting fsCD-GW method as well as the CVS-BSE method to sets of small organic and inorganic molecules and compare the results to experimentally obtained values.

A simple choice of frequencies z, for which the dielectric function κ and subsequently W are evaluated, are the true Gauss–Legendre roots used in Eq. (12). For the latter, W is readily available, causing no overhead in the algorithm. In order to stabilize the analytic continuation further, the imaginary frequency grid is extended by a set of reference points along the real axis.¹⁹ Duchemin and Blase employed a set of additional supporting points, linearly spaced along the real axis, and shifted by η into the imaginary plane, where η was chosen between 0.1 and 1.5 eV.¹⁹ In this work, we, instead, choose to use adaptive grids on the real axis as outlined in the following Sec. II C.

We note, in passing, that while for time-reversal symmetry breaking GKS references W^c(z) ≠ W^c(−z), it is still sufficient to construct Π_PQ(ω) only for those points of Z_s with $R(z)>0$ and then exploit the known symmetry relations of Π_PQ(ω).¹⁵ By constructing adaptive grids for individual systems, a balanced description is obtained for each system using a minimal number of actual frequencies, for which Eq. (16) needs to be evaluated.

Core excitations or regions with a high spectral density are usually inaccessible using standard linear response solvers that extract each excitation separately, that is, in a pole-by-pole manner. For molecules with more than a few atoms, core excitations are buried deep in the optical response spectrum, with hundreds of thousands of excitations being energetically lower and, therefore, preferably extracted by standard solvers based on the Davidson algorithm or similar eigenvalue solvers. Accordingly, special strategies need to be employed to tackle core excitations. Similar things hold true for spectral dense regions, as, for example, occurring in metal clusters or plasmonic systems, where even the visible range is already occupied by tens of thousands of single excitations that compose the optical response. In these cases, pole-by-pole solvers become inefficient, and one has to resort to different techniques.

The fsCD-G₀W₀ and fsCD-evGW methods outlined in Sec. II B and CVS-BSE method outlined in Sec. III A were implemented in a locally modified version of Turbomole 7.7.^60,61 The implementation makes use of OpenMP parallelism as well as being capable of running on modern graphical processing units (GPUs).^62–64 

To validate the implementation and accuracy of the fsCD-GW approach, the 1s ionization energy as well as the HOMO ionization energy of the trivalent lanthanoids from Ce³⁺ to Yb³⁺ were investigated. The PBE0 functional was chosen, in combination with the all-electron relativistic x2c-QZVPPall-2c basis set. Suitable auxiliary basis functions have been used for the RI approximation.^65,66 The GW@PBE0/x2c-QZVPPall-2c method was shown to perform well for these ions,¹⁵ and PBE0 is a broadly available density functional. Convergence thresholds of 10⁻⁹ hartree for the energy and 10⁻⁸ for changes in the atomic orbital density matrix were used. A grid of size 4 was used.⁶⁷ Calculations were done within a 2c formalism, including spin–orbit effects. Furthermore, a finite nucleus model was used.^68,69 For the evGW calculations, the spinors 55–80 were optimized in each iteration using (fs)CD-evGW, while the remaining spinors are being shifted accordingly (“scissoring”). Reference analytical results for the 1s spinors are obtained from the spectral $O(N6)$ representation using η = 10⁻³.⁷⁰ 

To test the scaling behavior of the fsCD-GW algorithm for core quasiparticle energies, the CD-G₀W₀@PBE0 and fsCD-G₀W₀@PBE0 quasiparticle energies of the 1s states are calculated for small cadmium clusters ranging from 4 to 18 atoms. The all-electron x2c-TZVPPall-2c basis set was used, in conjunction with a suitable auxiliary basis set.⁶⁵ Geometries of the cadmium clusters were taken from Ref. 71. All timings were measured on an Intel E5-2687W CPU, using eight threads. Kramers symmetry was exploited in 2c calculations.

To validate the different fsCD-GW approaches themselves, the core ionization energies of the CORE65 test set were evaluated.⁷² This set is composed of core ionization energies associated with 30 C, 21 O, 11 N, and 3 F 1 s states for a set of 32 small organic and inorganic molecules. For all molecules, core ionization energies were calculated from the fsCD-G₀W₀@DFT and the fsCD-evGW@DFT methods. All occupied states and the 20 lowest lying virtual states were optimized in the evGW treatment. The quasiparticle energies were evaluated for various density functional approximations. Specifically, they were evaluated for the pure DFAs BP86,^73,74 PBE,⁷⁵ TPSS,⁷⁶ and r²SCAN,^77,78 the hybrid functionals B3LYP, BH&HLYP,^74,79,80 PBE0,^75,81 and TPSSh,^76,82 the range-separated hybrid functionals CAM-B3LYP,^83, ωB97X-D,⁸⁴ and LC-ωPBE,⁸⁵ and the local hybrid functionals Lh12ct-SsirPW92,⁸⁶ Lh14t-calPBE,⁸⁷ Lh20 t,⁸⁸ and TMHF.⁸⁹ The triple-ζ x2c-TZVPPall-2c and quadruple-ζ x2c-QZVPPall-2c basis sets were used, respectively, in combination with a suitable auxiliary basis set.^65,66 Convergence thresholds of 10⁻⁹ hartree for the energy and 10⁻⁸ for changes in the density matrix were used. A grid of size 4 was used for pure, hybrid, and range-separated hybrid functionals,⁶⁷ while an ultrafine grid was used for local hybrid functionals.⁶² As only 1s states are investigated, spin–orbit effects are estimated to not be too important, and therefore, the spin-separated scalar relativistic exact two-component approach (X2C) is used. The finite nucleus model was used in all calculations. Note that unlike in Ref. 72, we do not add further corrections to the quasiparticle energies, as our model accounts for the most significant relativistic effects already. To assess the statistical error of each DFA, we use the same geometries and reference ionization energies as in Ref. 72.

To assess the performance of the CVS-GW-BSE approach, 40 core excited states of eight molecules associated with edges of the central metal (or metalloid) atoms, namely the L-edges of CrO₂Cl₂,⁹⁰ the L-edges of VOCl₃,⁹⁰ the K and L-edges of PdCl₂,⁹¹ the K and L-edges of SiCl₄,⁹² the K⁹³ and L-edges⁹⁴ of FeCp₂, and the K⁹⁵ and L-edges⁹⁶ of TiCl₄, TiCl₃Cp, and TiCl₂Cp₂, were investigated and compared to the respective experimental values. Convergence thresholds of 10⁻⁹ hartree for the energy and 10⁻⁸ for changes in the density matrix were used. A grid of size 5 was used for pure, hybrid, and range-separated hybrid DFAs,⁶⁷ while an ultrafine grid was used for local hybrid functionals.⁶² For the CVS approximation, all non-core spinors with an energy higher than the spinor of interest were frozen and projected out. A 2c formalism was used, including spin–orbit effects through the full X2C Hamiltonian. A finite nucleus model was used, and SOC corrections to two-electron integrals were estimated using the modified screened nuclear spin–orbit (mSNSO) approach.^97–99 As including spin–orbit effects via the 2c formalism, using spinors instead of orbitals, adds significant computational strain, we limit the investigation of this set of molecules to the G₀W₀-CVS-BSE@DFT and evGW-CVS-BSE@DFT levels of theory. The triple-ζ x2c-TZVPPall-2c basis set was used, in combination with a suitable auxiliary basis set.^65,66 The PBE0 and TMHF functionals were chosen as a suitable basis for testing, as PBE0 is again broadly available and performs reasonably well, while TMHF was found to be an excellent performer for the CORE65 test set, yet with currently limited availability. Finally, damped response BSE is used to get a grasp of the accuracy of the CVS approximation. For this purpose, we calculate the core excitation spectra in the energy range of the well separated L-edges of CrO₂Cl₂, VOCl₃, PdCl₂, and TiCl₃Cp using damped response BSE as outlined in Sec. III B and compare the obtained results with those from the CVS-BSE method.

To validate the frequency-sampled contour deformation based GW algorithm outlined in Sec. II B, the self-energies of the 1s orbital of open-shell lanthanoids are investigated. The corresponding correlation parts of the self-energies are listed in Table I.

As outlined in Table I, fsCD-G₀W₀ nearly perfectly mimics the results obtained from canonical CD-G₀W₀. Both also perfectly align with the analytical $O(N6)$ results obtained from the spectral representation.⁷⁰ Notable deviations are only observed for the 1s spinor of Nd³⁺, with (fs)CD-GW deviating by 0.02–0.04 eV from the analytical $O(N6)$⁠. This amounts to a deviation of 0.005%–0.009% of the total correlation part of the self-energy for the 1s spinors. For Nd³⁺, furthermore, the correlation part of the self-energy is rather large, especially when compared to the values found for the other lanthanoid ions. This is caused by a singularity being located in the vicinity of the quasiparticle peak. The latter leads to a rather steep slope of the correlation part in the region of interest and also to a stronger than usual dependence of the reference results on the chosen parameter η. A corresponding figure outlining this issue and further discussions of this case can be found in the supplementary material. Overall, fsCD-GW is in near perfect agreement with canonical CD-GW,^14,15 and both are in excellent agreement with spectral $O(N6)GW$⁠.⁷⁰ As time-reversal symmetry breaking GKS determinants correspond to the most challenging cases, it can safely be assumed that this excellent agreement will also be achieved for simpler non- and scalar-relativistic determinants. Furthermore, 2c and 4c GKS determinants that obey time-reversal symmetry are simply special cases of non-time reversal symmetric determinants in the framework of the corresponding fsCD-GW equations, and therefore, performance will be similar.

For valence ionization potentials of trivalent lanthanoid ions, which have been shown to be rather accurate when evGW@PBE0 is used,¹⁵ our findings are similar. Deviations of fsCD-evGW from CD-evGW are insignificant even when iteratively updating the quasiparticle energies as outlined in Table II. The only deviation being larger than 1 meV is found for Ce³⁺, with a deviation of 5 meV when compared to the CD-evGW reference. This is remarkable, as 16 spinors are updated in each iteration, outlining that the deviation of the correlation part of the self-energy from fsCD-GW is vanishing when compared to the reference canonical CD-GW implementation. The frequency sampling algorithm does an excellent job, providing basically the same accuracy as the canonical algorithm at a significantly reduced cost. This is further outlined for the 1s orbital of the water molecule in the supplementary material, for which fsCD-G₀W₀ recovers the behavior of the correlation part of the self-energy for a broad frequency range around the point of interest.

After having shown the validity of the fsCD-GW approach, we turn to investigate the improved scaling behavior of this algorithm. Figure 2 outlines the computational time needed to perform scalar relativistic (1c) and relativistic (2c) (fs)CD-G₀W₀ calculations on a set of cadmium clusters with 4–18 atoms. Only the 1s spinors were optimized, and the corresponding timings are outlined in Fig. 2. Two important conclusions can be drawn from Fig. 2. First, the overall scaling of fsCD-GW is vastly improved, with the algorithm being nearly an order of magnitude faster than the canonical algorithm. This behavior is observed for both references, 1c and 2c, with the latter being approximately eight times more expensive to calculate initially but then exhibiting basically the same scaling behavior. Second, there is no crossover between frequency sampled CD-GW and canonical CD-GW. The fsCD-GW algorithm is always at least as fast or faster than the canonical CD-GW algorithm, making it generally applicable without any drawbacks. Starting from ten cadmium atoms, equaling 480 active electrons, 2c fsCD-G₀W₀ even becomes faster than 1c CD-G₀W₀, outlining the large gain in performance from reducing the scaling by an order of magnitude. Interestingly, for both the standard CD-G₀W₀ and the fsCD-G₀W₀ algorithm, the asymptotic limits of $O(N5)$ and $O(N4)$ are not yet reached, respectively. This can likely be attributed to the time limiting step, i.e., evaluating Eq. (10), not being fully dominating in the investigated systems. For the investigated cadmium clusters, we further note that the fsCD-G₀W₀ results are remarkably accurate. Deviations to the canonical algorithm are found to be 10⁻⁵ eV or lower for all cluster sizes. This outlines that the gain in efficiency does not affect the accuracy of the method, rendering fsCD-GW as an extremely favorable algorithm for non-valence states in general.

The CORE65 set of molecules is composed of mainly organic molecules, for which the 1s core-electron ionization energies have been assessed by Golze and co-workers in Ref. 72. In the latter work, only PBE and a hybrid PBEh variant with 45% HF exchange were investigated, with the latter being found to work reasonably well together with the G₀W₀ method. An ad hoc atom-dependent correction was added to the 1s orbital energies in Ref. 72 to account for relativistic effects. In this work, the X2C Hamiltonian in combination with the mSNSO approach is used to account for relativistic effects, therefore also entering the GKS density, and no further corrections are applied. The resulting statistical deviations of the obtained core ionization energies with respect to experimentally obtained values for the fsCD-G₀W₀ method are shown in Fig. 3. Due to the different treatment of relativistic effects, it is noted that the obtained results may differ from those previously obtained.⁷² 

For the all-electron x2c-TZVPPall-2c basis set, surprisingly, the pure DFAs and DFAs with a small amount of Hartree–Fock (HF) exchange yield lower deviations from the reference results when used together with the G₀W₀ when compared to DFAs with high amounts of HF exchange. This is in contrast to valence shell quasiparticle energies, where DFAs with a rather high amount of HF exchange and especially range-separated hybrid DFAs perform significantly better than pure DFAs or DFAs with a small amount of HF exchange.¹⁰⁰ BH&HLYP yields error comparable to those of range-separated hybrids, with PBE0 being in between those two classes. Range-separated hybrid DFAs with pre-determined parameters fail to provide good starting points for core ionization energies, which can simply be rationalized by the fact that core orbitals are well localized, with no general need for range separation arising for them. Furthermore, it has been shown that optimally tuned range-separated hybrid functionals can far exceed the performance of their pre-fitted counterparts, and we expect this to be the case also for core excitations.^101,102 It is, furthermore, notable that replacing the x2c-TZVPPall-2c basis set by the quadruple-ζ x2c-QZVPPall-2c basis set yields improved core ionization energies for generally all hybrid functionals but worsens the agreement of the pure density functionals. Especially, BH&HLYP exhibits a large reduction in its overall deviation from experimental reference results as shown in Fig. 3(b), outlining that large basis sets are important when core ionization energies are tackled with DFAs that feature a high amount of HF exchange. The findings of Golze and co-workers that PBEh with 45% HF exchange also works well after a basis set extrapolation is, therefore, in good agreement with our data.⁷² In combination with the x2c-QZVPPall-2c basis set, the hybrid DFAs B3LYP and TPSSh and the metaGGA r²SCAN again provide excellent starting points.

For the class of local hybrids, TMHF excels for both basis sets, providing very good core-electron ionization energies. TMHF has high flexibility, allowing for incorporation of various amounts of non-local exchange depending on the density. This flexibility is key to an accurate description, allowing this functional to perform good performance. For fsCD-G₀W₀, it can, therefore, be concluded that r²SCAN, B3LYP, TPSSh, and TMHF provide suitable starting points when combined with triple-ζ and quadruple-ζ size basis sets, hinting that only a moderate amount of HF exchange is important. Quadruple-ζ basis sets generally improve accuracy, especially for DFAs that incorporate larger amounts of HF exchange. As TMHF also provides accurate valence ionization energies when combined with the GW method,^28,89 it provides an excellent starting point for GW calculations in general.

As shown in Fig. 4, partial self-consistency in the quasiparticle energies achieved through fsCD-evGW evens out most differences between the different GKS starting points. The only major outlier is the local hybrid TMHF, which remains significantly more precise compared to all other functionals when combined with fsCD-evGW. Using a larger quadruple-ζ basis set again significantly reduces all errors, with TMHF again providing the best performance by a good margin. Compared to fsCD-G₀W₀ results with the x2c-QZVPPall-2c basis set, the error of pure DFAs is significantly improved, but they still exhibit larger standard deviations than DFAs that include a low amount of HF exchange. Nevertheless, self-consistency in the quasiparticle energies fails to generally improve the accuracy of the GW method for core ionization energies. This is in line with earlier observations, finding errors of up to 1 eV for lower lying states.¹⁰³ Overall, the disappointing performance of evGW can be related to a worsened description of the density from the GKS reference in the vicinity of the nucleus. As evGW does not update the eigenvectors but is otherwise self-consistent, it is prone to amplify density errors. While this can be partly related to the density functional approximation itself, underscreening of the GW in the RPA approximation is another problem that becomes severe for core orbitals and could be corrected by the inclusion of vertex corrections.^103–105 

Compared to valence ionization energies, where hybrid and range-separated DFAs with high amounts of HF exchange perform best,¹⁰⁰  Figs. 3 and 4 give a qualitatively different picture. Pure DFAs and especially hybrid DFAs with moderate amounts of HF exchange perform rather well, while DFAs with high amounts of HF exchange perform comparably poor. Ab initio local hybrid DFAs as, for example, TMHF can bridge the gap and perform outstandingly well for both valence and core ionization energies.^89,106

While the 1s states of the CORE65 test set are not affected by spin–orbit coupling, SOC effects become rather important for certain interesting K- and L-edges in molecular systems containing metal centers. Therefore, we chose to investigate 40 core excited states of eight transition metal molecules as shown in Fig. 5, where the inclusion of SOC is mandatory due to state splitting.

The corresponding core excited states were calculated using the 2c CVS-BSE method outlined in Sec. III A, using PBE0 and a TMHF as GKS reference and quasiparticle energies from fsCD-G₀W₀ and fsCD-evGW. PBE0 is again chosen for its general availability and good performance, while TMHF is also compared due to being the single best density functional as outlined in Sec. V C. The resulting core excited states are compared to experimental values for which the statistical errors are provided in Table III. A detailed overview of all core excited states with the corresponding references can be found in the supplementary material. Overall, the results from the statistical evaluation presented in Table III are comparable to those obtained in Sec. V C, with fsCD-G₀W₀@TMHF being the best performing combination by a good margin. fsCD-G₀W₀@PBE0 performs significantly worse, especially for the high-energy K + L₁ edges. This is especially interesting given the fact that the initial starting GKS spinor energy of TMHF is assumed to be a worse approximation than the PBE0 one, as outlined in Ref. 15. As reference GKS determinant for a GW calculation, however, TMHF and other possible first principle local hybrid functionals are well suited, hinting at providing a reasonable description of the behavior of the electron density near the nucleus. The partially self-consistent fsCD-evGW variant again evens out the differences, with TMHF still performing superior to PBE0, but an overall reduced difference between both GKS references. The large maximum deviations found for all methods, but especially for PBE0 and fsCD-evGW@TMHF, can be traced back to the TiCpCl₃ complex, for which only fsCD-G₀W₀@TMHF yields reasonable values. The latter has a single outlier found for the L_2,3-edges of PdCl₂, which are, however, also present in all other combinations. Palladium L_2,3-edges are, therefore, by far the most troublesome in the investigated test set with deviations of up to 6 eV, which is an order of magnitude larger than the deviations found for other species. Subsequent testing hints that this error is a general failure of the GW-BSE method, as it is neither fixed by using the larger x2c-QZVPPall-2c basis set nor by dropping the CVS approximation as shown in Sec. V E. For palladium compounds, therefore, special care needs to be taken when those are tackled by the GW-BSE method, and an ad hoc shift parameter might be necessary. The remaining molecules, especially fsCD-G₀W₀, do a very good job predicting K- and L-edges, with mean absolute deviations of roughly 1 eV. The percentwise deviations outline that this is rather acceptable, given that core excitations are located in energy regimes of several 100 or even thousand eV. This is also confirmed when comparing our 2c fsCD-GW-BSE results to those of four-component (4c) damped response TD-DFT.³ In Ref. 3, Konecny and co-workers have tested various DFAs for their reliability in obtaining core excitations. Table IV compares the 4c damped response TD-DFT results with the fsCD-GW-BSE results obtained in this work.

When inspecting the values obtained for VOCl₃ and CrO₂Cl₂, it is notable that fsCD-GW-BSE has a significantly reduced sensitivity to the chosen starting point when compared to TD-DFT. The latter can vary by as much as 20 eV for the same excitation, depending on the chosen DFA. This is also in agreement with the previous findings of the authors.² Contrary, fsCD-GW-BSE varies by a maximum of 2 eV for the core excitations listed in Table IV, reducing this variance by an order of magnitude. The overall agreement of fsCD-GW-BSE with experimental reference values is excellent, especially for the fsCD-G₀W₀@TMHF method, which was already found to provide the best accuracy overall. GW-BSE should, therefore, be preferred over TD-DFT for core excitations whenever possible.

In this work, a relativistic frequency-sampled contour deformation based GW approach, labeled fsCD-GW, was outlined and implemented. fsCD-GW is an accurate yet efficient method to assess the quasiparticle energies of orbitals and spinors far away from the Fermi level, while fully taking spin–orbit coupling into account if needed. Computationally, fsCD-GW reduces the scaling by an order of magnitude when compared to canonical CD-GW and by two orders of magnitude if compared with spectral algorithms.^70,110 This computational advantage is accompanied by an intrinsically high accuracy, with results practically identical to those of the canonical algorithms. This high accuracy can even be maintained for the core and valence states of open-shell relativistic lanthanoid ions, which are among the electronically most complicated systems tackled by the GW approximation to date. Using fsCD-GW, the performance of various density functional approximations as reference GKS determinants for the GW approximation has been tested. Contrary to valence excitations, DFAs with a low amount of HF exchange, as well as the local hybrid functional TMHF, were found to be best suited to perform GW calculations for core ionized states. Furthermore, the fsCD-GW-BSE method was used to assess the core excited states of small molecules with metal centers. The core–valence separation approximation for the BSE was shown to yield reliable excitation energies and absorption cross sections for these core excited states. The combined fsCD-GW plus damped response BSE approach will also be a highly valuable tool for valence spectra with a very high density of states, as, for example, in nanoparticles.¹¹¹ Furthermore, recent developments in multiscale modeling of optical materials also rely on a similar approach,^112,113 and thus, the efficient fsCD-GW-BSE approach will make ab initio many-body multiscale simulations feasible with tools other than time-dependent DFT.

See the supplementary material for a discussion and graphical representation of the correlation part of the self-energy of the Nd3+ ion, the correlation part of the self-energy of the water molecule, individual core ionization energies of the CORE65 test set for all 15 tested density functional approximation, the individual values of the core excited states calculated for the eight small molecules in Sec. V D, and the corresponding optimized geometries of these eight molecules as ASCII coordinate files.

M.K. and W.K. gratefully acknowledge financial support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through the Transregional Collaborative Research Center “Cooperative Effects in Homo- and Heterometallic Complexes (3MET)” (Grant Nos. CRC/TRR 88 and Project C1). C.H. gratefully acknowledges the Volkswagen Stiftung for financial support.

The authors have no conflicts to disclose.

Max Kehry: Conceptualization (lead); Data curation (lead); Investigation (lead); Methodology (lead); Software (lead); Validation (lead); Writing – original draft (supporting); Writing – review & editing (supporting). Wim Klopper: Conceptualization (supporting); Supervision (lead); Writing – review & editing (supporting). Christof Holzer: Conceptualization (supporting); Data curation (supporting); Investigation (supporting); Methodology (supporting); Software (supporting); Validation (supporting); Writing – original draft (lead); Writing – review & editing (lead).

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