Challenges with relativistic GW calculations in solids and molecules

For molecules and solids containing heavy elements, accurate electronic structure calculations require accounting not only for electronic correlations but also for relativistic effects. In molecules, relativity can lead to severe changes in the ground-state description. In solids, the interplay between both correlation and relativity can change the stability of phases or it can lead to an emergence of completely new phases. Traditionally, the simplest illustration of relativistic effects can be done either by including pseudopotentials in non-relativistic calculations or alternatively by employing large all electron basis sets in relativistic methods. By analyzing different electronic properties (band structure, equilibrium lattice constant and bulk modulus) in semiconductors and insulators, we show that capturing the interplay of relativity and electron correlation can be rather challenging in Green's function methods. For molecular problems with heavy elements, we also observe that similar problems persist. We trace these challenges to three major problems: deficiencies in pseudopotential treatment as applied to Green's function methods, the scarcity of accurate and compact all-electron basis-sets that can be converged with respect to the basis-set size, and linear dependencies arising in all-electron basis-sets particularly when employing Gaussian orbitals. Our analysis provides detailed insight into these problems and opens a discussion about potential approaches to mitigate them.


I. INTRODUCTION
Computational methods in ab initio electronic structure theory have become indispensable tools in the design and study of new functional materials and molecules.For many years, mean-field methods such as density functional theory [1,2] (DFT) and Hartree-Fock (HF) were the commonly employed simulation tools mostly due to their low computational cost and reasonable accuracy.
However, recently, due to advances in computational resources as well as algorithms, sophisticated methods capturing electronic correlation beyond the mean-field approximation have gathered much more attention.
The development of new electronic structure methods is primarily driven by one or both of the following objectives: (i) accurate description of electron correlation, (ii) applicability to realistic systems with large number of electrons.Addressing both of these challenges satisfactorily is usually very difficult.For instance, low-scaling mean-field methods such as DFT and HF (with O(N 3 ) and O(N 4 ) computational costs, respectively, N being a measure of system size) can be applied to very large systems.However, they often lack quantitative accuracy.On the other hand, coupled cluster with single, double and perturbative triple excitation, i.e., CCSD(T), is considered as the gold standard for description of dynamical correlation, but has limited applicability due to its high O(N 7 ) scaling.Consequently, greater recovery of correlation effects usually comes hand-in-hand with an increased computational cost.
Even for elements in third and fourth row of periodic table, relativistic effects are important.[46] While for a full relativistic description, the 4-component Dirac equation, instead of the regular non-relativistic Schrödinger equation, should be solved, solving the 4-component equation is a formidable task and not strictly necessary for most systems if only the electronic structure properties are of interest.Several approximations have been proposed to introduce the full relativistic physics into an effective Hamiltonian that describes only the electronic degrees of freedom.The most widely used approximation follows the work of Breit, [47] where the Dirac equation is reduced to Schrödinger-like equation with an effective Dirac-Coulomb-Breit Hamiltonian that contains the relativistic terms.Theories such as the Douglas-Kroll-Hess [48][49][50][51][52] and the exact two-component [53][54][55][56][57] (X2C) introduce further approximations to make relativistic calculations achievable.
Several ab-initio implementations for realistic systems of the relativistic GW theory have emerged in the condensed matter community, [58][59][60][61] and in the quantum chemistry community, [37,62,63] where the latter are based on the X2C Hamiltonian.
Calculations in the X2C approach are generally performed using all-electron basis sets since current pseudopotentials (PPs) or effective core potentials (ECPs) do not recover relativistic effects such as spinorbit coupling (SOC).Most modern implementations or ab-initio relativistic codes are capable of handling relativistic elements, at least up to 5th or 6th row, including the Lanthanides, in the periodic table.Yet, performing reliable GW calculations for accurate electronic structure properties still remains an arduous task.Considering that the GW approximation is the simplest in the hierarchy of Hedin's perturbation theory, it offers an understandably limited accuracy in describing electron correlation.However, factors other than the accuracy of the GW method, e.g., poor quality and availability of basis sets, lack of effective PPs, etc., can have a larger impact on the quality of results, particularly in heavy elements.
In this discussion article, we analyze the factors that hinder the applicability of the sophisticated machinery of GW to both solids and molecules containing heavy elements.In particular, we first highlight the need for all-electron relativistic calculations due to the inability of PP to account for many important relativistic effects.We then investigate problems with achieving the basis set convergence in GW calculations.This is followed by results that highlight arising linear dependencies for larger basis sets.Such problems often aggravate the problem of slow convergence with respect to basis set size.Finally, we also comment on possible improvements of experimental results that are ultimately necessary for benchmarking and validating new theoretical methods.By investigating multiple properties of solids and molecules, we perform a comprehensive analysis of how several external factors, other than the accuracy of the method itself, influence the results in self-consistent GW (scGW ).
This article is organized as follows: in Section II, we recap some of the theoretical concepts on which majority of this paper is based.This is followed by computational and implementation details in Section III.
In Section IV, we analyze various hurdles, one by one, for the application of scGW to relativistic systems.A general discussion and conclusions on possible factors behind these challenges, along with the current and future developments that can address these problems, is presented in Section V.

A. Relativistic Hamiltonian: exact two-component formalism
The X2C formalism is widely used among electronic structure relativistic methods.[55][56][57]81] In this formalism, the four-component one-body Dirac-Coulomb Hamiltonian H 4C1e is used as a starting point on which a unitary transformation is performed that decouples the eigenvectors with large and small eigenvalues, generally labeled as large and small components, Here, V is the local Coulomb potential, T is the kinetic energy matrix, and the L operator is defined as The assumption invoked in the X2C formalism is that after the unitary transformation, the largecomponent can be effectively considered as electron, even though in principle, it mixes the electronic and positronic components.Therefore, the large component Hamiltonian H X2C1e + is used as the one-body Hamiltonian in electronic structure calculations, along with the usual two-electron integrals.We should note that no relativistic correction is used for the latter.
The formalism described so far is called the X2C1e approximation.
By definition, this approximation breaks the S z spin-symmetry and necessitates the use of generalized spin-orbitals in subsequent calculations.Another approximation, called as the spin-free X2C1e (sfX2C1e) is also popular.In sfX2C1e, only the spinindependent contribution arising form the L-operator in Eq. ( 2) is considered.Given that sfX2C1e leads to a symmetry-adapted Hamiltonian, for several systems, where negligible or no SOC is present, using sfX2C1e greatly reduces the computational cost and memory requirements in a calculation, without any accuracy loss.Even when SOC cannot be neglected, total energy trends in sfX2C1e do not differ significantly from X2C1e.[39] B. Self-consistent GW Hedin's equations [5] define a perturbation theory where the one-particle Green's function G, the vertex function Γ, the irreducible polarizability Π 0 , the screened Coulomb interaction W , and the self-energy Σ are related through a set of integro-differential equations.This is traditionally represented as Hedin's pentagram diagram.In the GW approximation, higher-order corrections to the vertex arising from G are ignored and the vertex Γ is approximated as a Dirac delta function in space-time.As a result, the self-energy is expressed as a product of the interacting Green's function G and the screened Coulomb potential W .
Hedin's equations or the GW approximation can be evaluated either on the real or imaginary axis.Here, we work with the finite-temperature (or Matsubara imaginary axis) formalism described in Refs.63 and 64.The one-particle imaginary-time Green's function is defined as where β = 1/k B T is the inverse temperature (k B being the Boltzmann constant), µ is the chemical potential, H and N are the Hamiltonian and particle number operators, c k pσ and c k, † pσ are the second-quantization annihilation and creation operators for the pth Bloch orbital with the spin σ and momentum label k, and τ ∈ [0, β] is the imaginary time.Finally, Z is the partition function, defined as In the GW approximation, the self-energy on the imaginary time/frequency axis is approximated as where Σ k ∞ is the static (frequency independent) Hartree-Fock self-energy, and Σ k c (iω n ) is the dynamic contribution, defined in the imaginary time (τ ) axis as where W is the effective screened Coulomb interaction, N k is the number of k-points sampled in the Brillouin zone (BZ).The summations are performed over momentum vectors q in BZ, as well as the atomic orbitals a and b in the unit cell.The new Green's function is then defined by the Dyson equation, where the chemical potential µ is fixed to ensure a correct particle number, S k is the overlap matrix, and H k 0 is the one-electron Hamiltonian.We refer the readers to Ref. 64 for exact expressions for the screened Coulomb operator W .The total electronic energy is calculated using the Galitskii-Migdal formula, where F k is the Fock matrix, defined as and ) is the one-particle density matrix.In our implementation of scGW , Eqs. ( 5) to (7) are solved iteratively until self-consistency is achieved in both the one-body contribution to the energy E 1b as well as the total energy E. Finite-size corrections for the static and dynamic contributions to the self-energy have been included, as described in Ref. 64.Relativistic effects at the level of both spin-free and full X2C formalism with one-electron approximation (sfX2C1e and X2C1e) have been used.

C. Equation of state for solids
For periodic solids, the relationship between energy and volume in a unit cell is called as the equation of state.Several different parametrizations of this equation are available.Fitting the energy-volume data to these forms allows us to estimate the equilibrium volume (and in turn, the lattice constant), the bulk modulus, and other related quantities.One of the most commonly used parametrizations, the Birch-Murnaghan [66] equation of state is defined as where E 0 , V 0 and B 0 denote the equilibrium energy, volume and bulk modulus, respectively, while B ′ 0 denotes the derivative of the bulk modulus with respect to the cell volume V , taken at V 0 .The Murnaghan [65] and Vinet [82] equations are other notable parametric forms.For the purpose of this paper, we fit the energy-volume data to Eq. ( 10) using the least-squares optimization.

III. COMPUTATIONAL DETAILS
All of our calculations are performed with Gaussian type orbitals (GTO) or their periodic generalizations.[83,84] For initializing all the calculations, we use PySCF [84,85] where we generate the one-and two-electron integrals, and perform mean-field calculations such as Hartree-Fock (HF) or DFT.For the density-fitting in twoelectron integrals, we employ even-tempered Gaussian orbitals [86] with progression factor 2.0 for both molecules and extended systems.For the relativistic calculations, we use x2c-nZVPall (n = S, T, Q) family of basis sets.[87] In contrast to the traditionally used uncontracted bases for the relativistic calculations, the x2c family of basis sets contains contracted bases minimizing the number of basis functions without a significant sacrifice of overall accuracy.To construct the initial Dirac Hamiltonian, the basis set is completely uncontracted (as is also programmed, by default, in PySCF).
All the dynamical quantities, i.e., imaginary-time and frequency dependent objects, in scGW are stored using sparse grids.[88][89][90][91][92][93] Sufficiently large grids are used to avoid loss of information in Fourier transforms performed over successive iterations of the GW cycle.To obtain spectral functions, IPs and band gaps, the Matsubara Green's function in our calculations is analytically continued to the real-frequency axis.
We use the Nevanlinna analytical continuation to accomplish this task.[94] To avoid poles on the real axis, a small broadening, a positive imaginary term, is added to the real frequency grid, i.e., For IP and band gap calculations, we employ η = 0.01 a.u. and η = 0.001 a.u., respectively.For photoelectron spectra of HgCl 2 , we have adjusted the chosen a broadening of η = 0.002 a.u. to match experimental peak widths, while in Fig. 1, η = 0.005 a.u. is used.All our calculations are converged with rigorous thresholds.
In particular, we ensure that the particle number is converged to 10 −8 and the total energy to at least 10 −5 a.u.The finite-temperature scGW code used here is accessible via GitHub.[95] A. Need for relativistic corrections Before we begin discussing the challenges associated with the application of relativistic scGW , it is important to justify the necessity of incorporating relativistic effects into correlated all-electron (AE) calculations.We accomplish this by studying the band structure of germanium.In Fig. 1, we plot three different band structure results evaluated using both PBE and scGW for germanium with diamond lattice, and lattice constant a = 5.657 Å.For all these calculations, we used a 6×6×6 k-mesh sampling in the BZ.
In the top left panel, we employ the GTH basis set and PP, where the relativistic effects are built directly into PP.Here, we find that PBE predicts a semiconducting phase with 0 eV band gap at the Γ point.The top middle and right panels show PBE calculation with and without the inclusion of scalar relativistic effects, respectively.Both were performed using an AE basis, namely x2c-SVPall.Right away, one finds that nonrelativistic AE-PBE yields a qualitatively incorrect result with an indirect and sizeable band gap of 1 eV between the Γ and L high-symmetry points.By including scalar relativistic effects at the sfX2C1e level, the band gap is closed resulting in a correct DFT band structure, similar to the GTH results.For germanium, SOC effects are known to be negligible and therefore, not considered here.
While both top left and top right panels of Fig. 1 display correct band structure of germanium at the PBE level, the experimental results yield band structure with an indirect band gap of 0.7 eV between Γ and L points, which only fortuitously somewhat resembles the band structure from the middle panel evaluated with AE basis and non-relativistic PBE.Consequently, to recover the experimental result for the right reason, both the relativistic treatment and inclusion of correlation beyond the DFT level is necessary.This is confirmed by the scGW results, shown in the lower panels of Fig. 1.In scGW , the band structure for AE basis set without any relativistic effects is further distorted such that the conduction band, specially near the Γ-point, is pushed upwards and an indirect band gap larger than 2 eV is predicted between Γ and L points.With appropriate inclusion of relativity, the AE scGW results provide a better quality band structure, with an indirect band gap of ∼1.3 eV, which agrees far better with the experimental value of 0.7 eV.

IV. RESULTS
We are now well-equipped to explore the challenges with relativistic GW calculations.
In each of the subsections, we describe a specific issue with supporting data.Further analysis and discussion on these results is presented in the Section V

A. Deficiency of pseudopotentials
Modern PPs are designed to account for relativistic effects and, in doing so, eliminate the need for including relativity in subsequent calculations.The GTH PPs and associated basis sets by Goedecker-Teter-Hutter [102,103] constitute one such example.While the use of PPs generally provide a good description of valence shell band structures, as can be observed from Fig. 1, they are far from adequate for other properties, generally derivable from total energy.We investigate the performance of AE and PP-based basis sets for three materials with diamond-lattice: Si, Ge, and αSn.With an increasing atomic number Z and core size, they serve as excellent  Si GTH-DZVP-MOLOPT-SR 5.499 5.456 5.431 [96]   84.8 96.6 97.9 [97]  x2c-TZVPall (sfX2C1e) 5.494 5.413 83.4 102.5 def2-TZVP (non-relativistic) 5.491 5.411 84.0 100.9 Ge GTH-DZVP-MOLOPT-SR 5.812 5.733 5.657 [98]  52.8 66.4 75.0 [99]  x2c-TZVPall (sfX2C1e) 5.769 5.609 56.9 83.9 αSn GTH-DZVP-MOLOPT-SR 6.711 6.649 6.489 [100]  52.8 32.6 53.1 [101]  x2c-TZVPall (sfX2C1e) representative examples to validate the importance of a proper description of the core orbitals.In Table I, by comparing PBE and scGW predictions for lattice constants and bulk moduli against their respective experimental values, we demonstrate that as the atomic number Z of compounds requiring relativistic treatment increases, and consequently the size of their core (approximated by PP) increases, the accuracy of PP-based results decreases, necessitating the use of AE basis sets to achieve experimental agreement.
For Si, which does not require relativistic corrections using PP, the AE basis set has virtually no impact on the quality of results.On the other hand, for αSn, AE results provide a significant improvement, which is even more pronounced in scGW than in PBE.We use GTH-DZVP-MOLOPT-SR basis set with GTH-PBE PP for the PP-based calculations, and x2c-TZVPall basis set, which is the largest basis where we could reliably perform equation-of-state calculations, for AE calculations.
For Si and other elements where the size of atomic core is small and the relativistic effects are negligible PPs are usually capable of yielding very good results.However, as the atomic number Z increases the problems with PPs become more noticeable as illustrated for Germanium and subsequently α-Sn.In Fig. 2, we showcase more examples where scGW evaluated in AE basis sets provides highly accurate results for both equilibrium lattice constants and bulk moduli.In contrast, the scGW results evaluated using GTH PPs are significantly worse than those obtained with PBE based on GTH PP.Particularly, in ZnX (X=S, Se, Te), AE relativistic scGW is highly effective in contrast to the GTH results.This poor performance of GTH PPs with scGW can be attributed to missing correlation contribution arising from the overlap of zinc's inner p and d orbitals with the inner p orbitals of the chalcogen.
In Green's function calculations such as scGW , when AE basis sets are employed, the core orbitals are described by self-energy matrix elements that have both static and dynamic parts.In contrast, when PPs are used, the dynamic part of self-energy is entirely missing from the description of core orbitals.For light elements such as Si, the magnitude of the dynamic self-energy for the core orbitals is small, and in addition, there are also only a few core orbitals.Consequently, the error of neglecting these self-energies is overall insignificant.On the other hand, given the large number of core orbitals in heavy elements, the number of dynamic core-core and core-valence self-energy matrix elements neglected by PPs is significant.Moreover, the mixing between valence and outer-core orbitals becomes important, resulting in sizable values of the outer core and core-valence selfenergy matrix elements.Since the calculations involving PPs, by design, cannot include any dynamic self-energy for the core, it is reasonable to expect that for elements where such self-energy contributions are significant, the use of PPs will result in sizeable errors.These errors can be additionally compounded by a lack of good, well optimized PPs for heavier elements.Therefore, it is worth noting that for Green's function methods, the use of PPs poses a different set of problems than for DFT.The DFT results can only be affected by a lack of well optimized PPs, while Green's function methods would always be penalized by the absence of dynamic self-energy even for the best optimized PPs.For some lighter elements, one can only hope that these dynamic self-energy matrix elements have insignificant values, however, this assumption cannot be generally fulfilled.

B. Convergence issues in basis sets
For heavy elements, the very large size of the one-electron basis sets effectively hinders the ability to investigate convergence with respect to basis set size in AE calculations.While for molecules, this formidable challenge still can be overcome, it becomes insurmountable in calculations of solids.The problem is relatively mild, often almost non-existent, at the level of DFT or HF, but quickly becomes aggravated for correlated calculations such as scGW .This can also be understood by the fact that basis sets are usually optimized at the atomic level using SCF calculations.Such behavior has also been observed by Kutepov in FIG. 2. Errors in Top: lattice constants a, and Bottom: bulk moduli B0 for selected materials.All compounds are calculated with GTH-DZVP-MOLOPT-SR, x2c-SVPall and x2c-TZVPall basis sets.Both AE-PBE and AE-scGW results are reported using the sfX2C1e Hamiltonian, with a 4 × 4 × 4 k-mesh sampling.Experimental values for lattice constants and bulk moduli can be found in Refs.[104,105] and [106][107][108], respectively.TABLE II.PBE0 and scGW results for ionization potentials (in eV) for selected silver halids.The results were evaluated using different basis sets and hten extrapolated to the basis set limit.Notice that the mean-field method appears to be converged already at the triple-ζ level.Ref. 44, where DFT calculations converged far more rapidly than scGW while using linearized augmentedplane-wave (LAPW) basis functions.Here, we examine such convergence issues in detail.First, we look at IPs for AgCl, AgBr and AgI, which are representative examples of closed-shell relativistic molecules with an increasing SOC as we go from Cl to I.These systems have been widely used to benchmark recent implementations of relativistic GW methods.[37,39,60] For these molecules, in Table II, we show IPs evaluated in PBE0 and scGW with x2c-TZVPall and x2c-QZVPall basis sets.
Using these two IP values, an extrapolated IP value (also reported in Table II) is calculated assuming an inverse relationship between IP and the number of AOs in the respective basis set.We use the scalar relativistic Hamiltonian in all these calculations.
Comparing these trends  II.All IPs are reported relative to the respective extrapolated values, which shows that scGW exhibits a much slower convergence than PBE0.
graphically in Fig. 3, we can clearly observe that for PBE0, the extrapolated IPs are within 0.1 eV of QZ, displaying a fast convergence.However, for scGW , the quadruple-ζ results are far from converged, with extrapolation leading to changes of more than 0.2 eV.It is worth mentioning that an extrapolation based on two calculations, triple-and quadruple-ζ, is already far from ideal, specially when the resultant corrections are so large.Having quintuple-ζ results would greatly enhance the quality of results.However, even not considering the computational challenge that such a calculation poses, in many relativistic basis-sets, including x2c, quintupleζ bases are simply not available.
Next, we look at the band structure in relativistic solids, with focus on AgBr and CdSe.For AgBr, rock salt crystal structure with a = 5.774 Å is considered, while the zincblende phase in CdSe with a = 6.05Å is used.All calculations were performed using a 4 × 4 × 4 k-mesh sampling, using x2c-SVPall, x2c-TZVPall and x2c-QZVPall basis sets.An inverse temperature of β = 300 a.u.−1 was used.For the QZVPall basis in AgBr, diffuse s and p orbitals with Gaussian exponents smaller than 0.05 were removed to avoid linear dependencies in the basis set.High angular momentum g-orbitals were removed from QZVPall basis set in both systems.
For in Ref. 63, Yeh et al. noted poor convergence of band gap values along the path containing special symmetry points.Here we present further data contrasting scGW convergence trends against PBE.
In Table III, we present direct and indirect band gaps in AgBr and CdSe, with x2c-SVPall, x2c-TZVPall and x2c-QZVPall basis sets.Already at the triple-ζ level, for both systems, PBE band gaps are well converged to within 0.05 eV when compared with quadruple-ζ basis.However, scGW once again shows much slower convergence trends, and going from triple-to quadrupleζ still introduces corrections larger than 0.1 eV for CdSe and 0.2 eV for AgBr.
Finally, in Table IV we show convergence trends for lattice constants for Si, Ge, GaAs, αSn, and InSb.All these systems have a diamond (zincblende for GaAs and InSb) structure, with an increasing size of atomic core as we go from Si to InSb.Similar to previous observations, the convergence of scGW values is more than 2-3 times slower than for PBE.At the level of DFT, it seems that only for Si the larger basis set marginally improves the agreement with experiment, while for all the other materials, the results either do not change or become worse.Consequently, one can suspect that DFT optimization of higher level bases is not as effective as lower level bases since the convergence with respect to the basis size is achieved relatively early on.This suspicion is somewhat confirmed when examining the scGW results.Almost all SVPall scGW results are surprisingly close to experiment and do not seem to be improved in a larger basis.For periodic systems, an extrapolation procedure, similar to the one used for molecular IPs, cannot be considered.This is because SVP is not a sufficiently large basis and using it together with only the TZVP basis would lead to largely inaccurate extrapolated data.While for molecular systems, QZ-level calculations remain accessible and make extrapolations possible, for the periodic systems, these calculations are extremely expensive and mostly out of reach.

C. Linear dependencies
For correlated calculations, the quest for converging properties to a complete basis set limit is not only affected by slower convergence and lack of larger basis sets but also by linear dependencies arising in the large basis sets.Particularly in solids, due to close packing of atoms, adding more basis functions introduces strong linear dependencies among the atomic orbitals.As a result, the condition number of the overlap matrix becomes too large even for an easy execution of meanfield methods.The problem of linear dependencies is not unique to GTOs considered here, and has been reported in other kinds of orbital representations, namely Slatertype orbitals (STO), and also in LAPW.[36,[110][111][112] Even though linear dependencies are more severe in ionic compounds, they also remain problematic in covalent compounds (such as silicon) and arise at small bond distances in all compounds.One way to remedy this problem in GTO basis sets is to simply remove the diffuse functions.These orbitals have long tails because of their small Gaussian exponents, leading to an increased overlap with other atomic orbitals.However, when studying bulk properties, such a procedure generally does more harm than good.When the bonds are stretched, the lack of diffuse orbitals leads to a deficient description of atomic bonds and substantially alters the equation of state, causing significant errors in both the equilibrium lattice constant and bulk modulus.This phenomenon is investigated for silicon in Fig. 4, where we have considered a rather simple case of non-relativistic calculation involving the widely used def2-TZVP basis Experiment [96, 97]  5.431 97.8   [96]   Ge 5.779 5.773 -0.006 5.622 5.559 -0.063 5.657 [98]   GaAs 5.747 5.748 0.001 5.600 5.559 -0.041 5.653 [104]   α-Sn 6.620 6.663 0.043 6.456 6.391 -0.065 6.489 [100]   InSb 6.597 6.64 0.043 6.463 6.381 -0.082 6.479 [109]   set.Here, all orbitals with Gaussian exponents smaller than 0.1 were dropped, leading to lattice constants differing by more than 0.1 Å and over 50 % error in bulk modulus.Note that an ad-hoc re-introduction of diffuse functions as we stretch the bonds introduces jumps in the energy-volume curve, precluding a reliable fit to the Birch-Murnaghan equation.Even though the computational cost is still affordable, due to arising linear dependencies, converging DFT and GW calculations for Si with a larger basis set proves extremely challenging and is therefore not pursued here.

D. Lack of accurate experimental reference
To benchmark new theoretical methods that are pushing the limits of accuracy and applicability, reliable experimental references are necessary.Here, we consider the HgCl 2 molecule as an example to highlight a potential area were better experimental data can provide feedback to improve theoretical methods.In Fig. 5, we compare scGW results with two different experimental spectra for HgCl 2 .In one of the experiments by Borgges et al. [113] (labeled Exp-1 in the plot), the SOC splitting is unresolved in both first and second IPs.Results by Eland et al. [114] (Exp-2 in Fig. 5) clearly resolve the SOC splitting in the first IP ( 2 Π g state), but reports difficulty with the second peak ( 2 Π u state).Among GW results, scalar relativistic corrections (sfX2C1e) capture the overall peak structure, producing a spectra similar to Exp-1.However, the spin-free theory understandably misses out on the SOC splitting which, in turn, is captured at the X2C1e level.In fact, for scGW with X2C1e, spin-orbit splitting is visible for both 2 Π g and 2 Π u states.For the first ionization peak, scGW predicts a 0.16 eV splitting for 2 Π 3 2 g , 2 Π 1 2 g states which is close to the experimental value of 0.12 eV.Similarly, for the 2 Π 3 2 u , 2 Π 1 2 u state, a splitting of 0.12 eV is predicted, but the corresponding results are not resolved in experiments.

V. DISCUSSION AND CONCLUSIONS
We have presented ample numerical evidence that when employing Green's function methods, particularly scGW , the current frameworks of PPs and orbital representation encounter challenges in fully realizing the potential of precise quantum chemistry calculations for relativistic systems.One of the challenges is that PPs are not adept in accurately describing properties other than valence-shell electronic structure such as band gaps and ionization potentials.When a material is pressurized, e.g., while studying the equation of state in a material, the correlation effects arising from the overlap between atomic cores cannot be captured by PPs, and highly expensive AE calculations become unavoidable.Effective cores and PPs where a smaller number of electrons are modelled as part of the core offer one way to overcome this hurdle, and while some candidates are available, [115][116][117][118][119] more developments and systematic benchmarks of new ECPs and PPs are desirable.
Additionally, it should be noted that when using PPs in Green's function methods, there is an inherent error that comes from the absence of dynamic self-energy matrix elements for core and core-valence sectors even if the best optimized PPs are employed.Especially for heavier elements, one can surmise that such a PP approximation will lead to significant errors since both the values and the number of outer core dynamic selfenergy matrix elements may be large.Consequently, even with the most optimized PPs, as they are formulated now, we can expect significant errors.
For correlated methods such as scGW understanding the basis set convergence is a bit difficult.
Even when well designed and optimized basis sets are available, [120], correlated methods generally exhibit a slower convergence than mean-field ones.Yet, the fact that basis sets are optimized at the level of meanfield cannot be overlooked and, in fact, is potentially one the reasons behind faster convergence in DFT and HF.[121,122] For lighter elements, basis set convergence is easily achieved and such problems do not require any particular attention.However, same is not true for heavy elements where relativity is also important.The issues with basis sets for heavy elements is well known and has been a topic of recent research.[123,124] Furthermore, most quantum chemistry basis sets have been optimized for molecular applications.When these are applied to solids, not only the molecular level optimization may prove insufficient, additional problems related to linear dependency and unfavorable convergence trends are more prone to arise.Recently, Ye et al. [125] also noticed similar problems and proposed correlation-consistent basis sets optimized for solids.Another way to improve the quality of basis sets is to perform a material-specific optimization.[126] These developments, including other means to eliminate linear dependencies, [127,128] are yet to be applied to relativistic materials and molecules.
It is also likely that a better one-particle orbital representation, with strictly orthogonal atomic orbitals, may be required to simultaneously overcome the issues related to convergence and linear-dependencies. [129][130][131] Consequently, one can expect that in the future explicitly orthogonal single-particle orbitals, such as tensor-train numerical STOs, [132] or orthogonal Gausslets, [129][130][131] will show promise to solve the problem of linear dependency.Constructing PPs and basis functions from correlated calculations also may offer a viable route to improve the reliability of relativistic calculations.
Lastly, we also show that as theoretical methods become more accurate, better reliable experimental benchmarks data with well resolved photo-electron spectra are desirable to both validate theoretical results and push for new theoretical developments.

FIG. 1 .
FIG. 1. Band structure of germanium calculated with Top panels: PBE and Lower panels: scGW theories, using Left: GTH-PBE pseudopotential with GTH-DZVP-MOLOPT-SR basis set, Middle: x2c-SVPall all-electron basis set with no relativistic corrections, and Right: x2c-SVPall basis set with sfX2C1e Hamiltonian.For the GTH result, relativity is intrinsically accounted for in the pseudopotential.Diamond lattice with lattice constant a = 5.657 Å was used, and all calculations were performed with a 6 × 6 × 6 k-mesh sampling in the BZ.

FIG. 3 .
FIG. 3. Graphical representation of convergence trends inTableII.All IPs are reported relative to the respective extrapolated values, which shows that scGW exhibits a much slower convergence than PBE0.

FIG. 4 .
FIG. 4. Left panel: Lattice constants a and bulk modulus B0 for silicon, calculated using def2-TZVP basis with and without diffuse orbitals (Gaussian exponent < 0.1).Both PBE and scGW results, along with experimental values, are shown.Nonrelativistic Hamiltonian was employed for these calculations.Right panel: The corresponding energy-volume curves.

FIG. 5 .
FIG.5.Comparison of experimental spectra with scGW results for HgCl2, obtained with both scalar and twocomponent relativistic effects.For these calculations, we have used the x2c-QZVPall basis.Both experimental spectra, exp-1 by Bogges et al.[113] and exp-2 by Eland et al.[114], are shifted by 0.35 eV to match the associated IP peaks.

TABLE I .
Comparison of the predicted lattice constants and bulk moduli from calculations with AE basis sets against those with PP.The data highlights the deficiency of PP in describing energy-derivable bulk properties as atomic number Z increases resulting in an increase of the PP approximated core.

TABLE III .
PBE scGW PBE scGW PBE scGW Convergence of PBE and scGW band gaps for AgBr and CdSe, with respect to basis sets from the x2c family.The results were calculated with the sfX2C1e-Coulomb Hamiltonian and a 4 × 4 × 4 k-mesh and inverse temperature β = 300 a.u.−1 (for scGW ).The total numbers of GTOs per cell for each basis set are listed in the second column.

TABLE IV .
Basis set convergence trends in lattice constants [ Å] for selected compounds.Results for x2c-SVPall and x2c-TZVPall basis sets are shown, denoted by SVP and TZVP, respectively.For readers' convenience, we also list the difference between them.