Assessing MP2 frozen natural orbitals in relativistic correlated electronic structure calculations

The high computational scaling with the number of correlated electrons and the size of the basis set is a bottleneck which limits applications of coupled cluster (CC) algorithms. This is particularly so for calculations based on 2- or 4-component relativistic Hamiltonians, which often employ uncontracted basis sets and lead to large virtual molecular orbital (VMO) spaces. This problem may be alleviated by employing a more compact set of virtual spinors than those provided by the canonical Hartree-Fock (HF) set, such as those based on natural orbitals (NOs). In this paper we describe the implementation of a module for generating NOs for correlated wavefunctions, and in particular MP2 frozen natural orbitals (MP2FNOs), as a component of our novel implementation of relativistic coupled cluster theory for massively parallel architectures [J. Pototschnig et. al., J. Chem. Theory Comput. 17, 5509, 2021]. Our implementation is capable of manipulating both complex and quaternion density matrices, thus allowing for the generation of both Kramers-restricted and Kramers-unrestricted MP2FNOs. Furthermore, NOs are re-expressed in the parent atomic orbital (AO) basis, so that the code also makes it possible to generate CCSD natural orbitals in AO basis for further analysis. By investigating the truncation errors of MP2FNOs for both the correlation energy and molecular properties at CCSD level such as the electric field gradients at the nuclei (EFGs), electric dipole and quadrupole moments for hydrogen halides HX (X=F-Ts), and parity-violating energy differences (PV) for H$_2$Y$_2$ (Y=O-Se), we find that MP2FNOs accelerate the convergence of the correlation energy in a roughly uniform manner across the periodic table and that, with VMO spaces truncated to around half the size of the full spaces ones, it is possible to obtain reliable estimates for both energies and all molecular properties considered.


I. INTRODUCTION
Understanding the electronic structure of large molecules or complexes containing heavy elements such as lanthanides or actinides is a problem of relevance for many technological applications. Examples are the nuclear fuel cycle, such as in the development of new extractants for separation processes [1][2][3][4][5][6] , and the use of lanthanides as exceptionally strong single molecule magnets 7 . To model such materials, Density Functional Theory (DFT) has become the most widely used correlated electronic structure theory approach 8 , even though it is difficult to systematically approach exact results with the currently available density functional approximations 8,9 . In the particular case of relativistic electronic structure calculations, DFT energies may even for closed-shell species strongly deviate from experimental or accurate theoretical results 10-12 . This also holds for molecular properties, recently Sunaga and Saue 13 reported that the performance of DFT for parity violation energy shift (PV) calculations -a property requiring a very accurate description of the electronic wave function near the nuclei -is somewhat disappointing, with deviations to CCSD being as large as 10%. These uncertainties in the performance of of DFT for heavy elements, especially for cases in which experimental values are absent or difficult to generate, calls for the use of state-of-the-art wavefunction methods to either provide accurate reference data, or be applied directly if the model sizes are small enough.
Among currently available approaches, coupled cluster a) Electronic mail: xiang.yuan@univ-lille.fr b) Electronic mail: l.visscher@vu.nl c) Electronic mail: andre.gomes@univ-lille.fr (CC) theory serves as a "gold standard" 14,15 for its ability to give results approaching chemical accuracy for both correlation energies and properties. However, the main difficulty in employing CC theory in large-scale applications (such as molecules containing several hundred electrons) is its high computational scaling with respect to the size of the system (N). For example, without approximations CCSD and CCSD(T) approaches scale with O(N 6 ) and O(N 7 ), respectively, and also with approximations such as Laplace transforms 16 , distance screening 17 and density fitting 1819 the scaling and prefactors are still significantly higher than in a mean-field approach like DFT.
Another key ingredient in modeling heavy element species is the treatment of relativistic effects [20][21][22] . For cases in which a given molecular property is not particularly sensitive to effects such as spin-orbit coupling (SOC) (which is often so for molecular structural parameters such bond lengths and angles), or to contributions from electrons other than those in the valence (e.g. dipole moments), approximate treatments of relativity can be employed, through pseudopotential approaches 23,24 , or by including only scalar relativistic effects. For more challenging applications, or higher precision, a more general framework can be based on the solution of the 4component (4C) molecular Dirac equation [25][26][27][28] . Nowadays this approach is often made more tractable by solving the exact 2-component (X2C) equation that can be derived from the Dirac equation after a basis set discretization [29][30][31][32][33][34][35][36][37][38][39] . Both the original 4C approach as well as its X2C approximation can deliver accurate molecular properties across the periodic table, also for properties involving core electrons 10,13,[40][41][42][43][44][45] .
The one-electron functions (molecular spinors, for simplicity also referred to in the following as molecular orbitals) obtained from solving the 4C or X2C matrix equations, serve then as a basis for a correlation treatment in the so-called no-pair approximation in which contributions to the correlation energy due to admixture of states with explicit electronpositron pairs are neglected 46 . As the computational cost of 2C approaches is lower than 4C ones, the choice between which treatment to use prior to the correlated treatment will depend on a case-by-case analysis of whether the additional cost of the latter will be offset by improvements in accuracy over the former.
Recently, Pototschnig et al. 47 described a new, efficient relativistic coupled cluster implementation based on ExaTEN-SOR 48 , a distributed numerical tensor algebra library for GPU-accelerated HPC platforms. This code enabled the calculation of molecular properties with CCSD wave functions for systems such as [(UO 2 )(NO 3 ) 3 ] − , for which 200 electrons and around 1000 virtual molecular orbitals were included in the correlated treatment. Compared to nonrelativistic implementations, it is in particular the introduction of SOC which increases the computational cost (though not altering the overall scaling) 49 . Not only does this necessitate the use of complex algebra, it also makes use of contracted basis sets more difficult as one needs to be able to describe the differences in radial extent of spin-orbit split orbitals. This effectively doubles the number of functions that is needed to describe the p-, d-and f -type core orbitals 50 in the case of implementations based on expressing spinors in scalar basis sets 26,28,50 . In this framework another complication arises in how to define contractions for the small component part of the spinors that would respect the kinetic balance condition. We note that these problems do not arise in implementations capable of directly evaluating 2-electron integrals over 2C spinor basis functions 27,51 .
The lack of efficient contraction schemes for scalar basisderived implementations, combined with the fact that in many potential applications it may be required that all, or a significant part of core electrons (e.g. in core spectroscopies 45 ) are treated explicitly, including the deep core orbitals that require many Gaussian type functions to be described correctly, leads to calculations having large virtual spaces. This directly affects the performance of CC-based approaches as these scale with the fourth power of the number of virtuals (O(V 4 )). In practice one therefore often alleviates the computational effort to some extent by leaving out high energy virtuals that are mostly localized in the core, but convergence with respect of such an energy cut-off threshold is slow, especially when semi-core or core correlation need to be described as well.
In non-relativistic electronic structure theory it has long been recognized that, in contrast to canonical Hartree-Fock molecular orbitals (CMO), natural orbitals (NOs)-the eigenvectors of the one-body reduced density matrix (1RDM) 52,53provide a more compact and quickly converging orbital representation for describing post Hartree-Fock wavefunctions. Based on this observation, the idea of replacing the CMOs by NOs to reduce the size of the virtual space and thereby computational cost was introduced 54 .
Rather than using an energy threshold as is done by CMOs, one may instead omit NOs which are likely to not strongly contribute to the total correlation energy from the virtual space. This is done on by considering the magnitude of natural occupation numbers of approximate NOs, obtained from a method that quickly gives access to a reasonable approximation of the 1RDM of the correlated wave function. For this purpose, second order Möller-Plesset perturbation theory (MP2) is a particularly appealing approach to obtain the 1RDM and the approximate NOs because of its low scaling, non-iterative O(N 5 ), and the ability to recover most of correlation effects. Within the natural orbital family of methods, the virtual frozen natural orbitals (FNOs) approach 54-56 has therefore become popular because of its clear concept and simple implementation. FNO theory has been thoroughly discussed for non-relativistic models such as configuration interaction (CI) 54 , multi-configuration self-consistent field (MCSCF) 57,58 and coupled cluster 56,59 . Recently, Verma et al. 60 furthermore extended the FNOs algorithm to quantum computers.
The ability of virtual FNOs to reduce computational cost compared to CMOs is even more appealing in relativistic electronic structure calculations, in particular in the case of contracted basis sets and/or when correlating sub-valence electrons. The main goal of this paper is therefore to describe and showcase the implementation of MP2-based FNOs (MP2FNOs) within the framework of the new relativistic coupled cluster implementation for massively parallel, GPUaccelerated platforms 47 . While our code primarily aims at describing cases in which spin-orbit coupling is taken into account, we also demonstrate its use with a non-relativistic Hamiltonian.
Our second aim is to discuss the performance of MP2FNOs across the periodic table, by treating model systems containing elements ranging from first-row such as fluorine, to superheavy elements such as tennessine. Here, we shall focus in how well truncated FNO spaces describe both correlation energies as well as first-order properties such as the electric field gradient at the nuclei (EFGs), parity-violation energy shifts, and electric multipoles (dipole and quadrupole). The first two properties are chosen as representatives of properties for which a good description of wavefunctions in the region close the nuclei is important (even if only valence electrons are correlated), while the electric dipole and quadrupole moments are an example of properties for which the major contributions arise from valence electrons.
Besides being able to improve the efficiency of a calculation, natural orbitals are also interesting as tools for analysis, such as in estimating the multi-reference character of a system 61 . Their visualization in real space can furthermore provide insight into correlation effects, since by their one-particle nature they are easier to interpret than the wavefunction itself. Thus, another aim of our implementation was to provide a tool for obtaining natural orbitals for any correlated wavefunction, from a 1RDM. In this paper we shall make use of the analysis of the CCSD natural orbitals, and in subsequent work we plan to further explore the use of natural orbitals in visualization.
This article is organized as follows. In Sec. II the background of MP2FNOs theory is summarized. In Sec. III we described the details of implementation. All the sample calculations are presented and discussed in Sec.IV-V. And finally a brief summary is given in Sec. VI.
We note that upon completing this manuscript, we have become aware of another implementation of the MP2FNOs approach for relativistic correlated methods, in the BAGH code 62 . While the main features of the MP2FNOs method are the same in both implementations we first note that our implementation fully exploits ExaTENSOR's single-node or distributed memory (multi-node) and GPU acceleration capabilities, and as such can be efficiently employed in systems ranging from local clusters to latest-generation supercomputers. Second, as it will be outlined below our implementation is capable of manipulating both complex and quaternion density matrices, thus allowing for the generation of both Kramersrestricted and Kramers-unrestricted MP2FNOs. Finally, it allows for re-expressing NOs in atomic orbital (AO) basis for further analysis.

II. THEORY
As the MP2FNO approach is well-known in a non-relatistic context and requires essentially no modification for application in a no-pair relativistic context, we will only provide a brief description. We apply the orbital-unrelaxed MP2 approximation (for working equations for the orbital-relaxed formalism in a relativistic context see reference 63). The second-order contribution to the occupied-virtual block of the density matrix from single excitations is zero for canonical orbitals 53 while the relaxation contributions to this block are ignored in orbital-unrelaxed MP2. This approximation thus decouples the occupied and virtual spaces and allows us to focus only on obtaining the virtual-virtual block of the density matrix. We want to keep occupied orbitals in their canonical form and therefore will not determine the occupied-occupied part of the density matrix, for which we simply retain the diagonal Hartree-Fock form D (0) oo = 1 oo . With these approximations the second-order FNOs density matrix is given as with the relevant matrix elements given as Here, following usual conventions I, J denote occupied spinors and A, B virtual spinors, and ε AB IJ = ε I + ε J − ε A − ε B , is the energy difference between canonical Hartree-Fock spinors. In our implementation, these matrices are obtained without imposing any time-reversal symmetry on the spinor set and are thus valid for both Kramers-restricted and Kramers-unrestricted Hartree-Fock procedures. If timereversal 64 and spatial symmetry 65 are not present, these matrices are in general to be represented in complex algebra.
The simplest procedure to obtain natural spinors is to apply a diagonalization in complex algebra after D (2) vv is formed. However, for closed shell systems in which time-reversal symmetry is enforced for the orbitals at the mean-field level, as is the case in the DIRAC code, one would like to retain such Kramers pairing also for the natural spinors. This requires careful attention when degeneracies beyond the twofold Kramers degeneracy are present, as complex diagonalization will arbitrarily mix the degenerate solutions such that Kramers pairing is not guaranteed.
In such cases, Kramers pairing of the natural spinors can be enforced by first transforming D (2) vv from complex to quaternion algebra D where the tranformation matrix U is given by 66 and I is a unit matrix of dimension n × n, n being the number of Kramers pairs. This transformation will block diagonalize D (2) vv , leading to two decoupled matrix problems of half the original dimension in quaternion algebra.
Diagonalization of one of the two blocks provides a unique set (V) of quaternion FNOs in MO basis and their respective occupation numbers (O) from which the solutions for the other block can be generated via Kramers symmetry.
At this stage, we can employ the information from the occupation numbers to reduce the original virtual space V toV by removing from V orbitals with occupation numbers lower than a user-defined threshold. The larger such threshold is, the smaller the remaining set will be. In the numerical examples presented in the following we will try to find the optimal balance between efficiency (smaller sets) and accuracy (recovering more of the correlation energy and better reproducing first order properties).
After the basis is truncated, it is convenient to recanonize the remaining orbitals. To this end, we perform two operations: first we transform the virtual-virtual block of the Fock matrix (F vv ) into the truncated FNOs basis (F vv ) and then do a (quaternion re)diagonalisation to obtain a new set of canonical orbitalsW and orbital energies¯ .
The productVW of these two transformation matrices gives the transformation that expresses the recanonized truncated natural orbital set in the atomic orbital(AO) basis. These orbitals can then be used in any subsequent correlated calculation as replacement to the original Hartree-Fock orbitals. In summary the transformation of the original Hartree-Fock orbitals U to the new setŪ is thus given by whereŪ U occ = U occ (10) and the dimension of the rectangular matrixV and that of the square matrixŪ depends on the truncation threshold.

III. IMPLEMENTATION
The aforementioned algorithm has been implemented in the relativistic quantum chemistry package DIRAC 26 , as part of the ExaCorr code 47 . Our implementation allows for calculations to be carried out either using a single-node configuration (for which ExaTENSOR provides OpenMP parallelization, on top of GPU offloading) or employing distributed memory in the case of multi-node runs.
Our implementation is structured following its three main tasks: one module deals with the construction of density matrices in the MO basis, another carries out the complex or quaternion diagonalization of density matrices, and the final module takes care of the construction of the recanonized MP2FNOs. The final step should be repeated if the truncation treshold is changed, the first two are independent of this threshold.

A. Density matrix construction
In Exacorr all computationally expensive operations such as the tensor contractions used in the determination of amplitudes or the construction of density matrices are offloaded to ExaTENSOR library.
In the case of MP2FNOs, we have created a dedicated driver for MP2 calculations, in which we (i) calculate and store in memory two electron repulsion integrals (ERIs) in AO basis with the aid of the InteRest library 67 ; (ii) employ the standard Yoshimine scheme 68 (which scales as (O(N 5 )) to transform the AO integrals to the MO basis. Only the direct and exchange contributions to the i j| |ab (OOVV-type) integrals are calculated, which makes the integral transformation step much faster than the generation of a complete set of MO integrals as is normally done in ExaCorr. After antisymmetrization the i j| |ab integrals are stored (in memory) as ExaTENSOR tensors, which may reside on a single node or be distributed over several nodes. Finally, (iii) the MP2 amplitudes tensor and the MP2 energy are determined and the density matrix D (2) vv is constructed according to equation 2.

B. Complex-quaternion transformation and diagonalization
The density matrix computed in Exacorr is generated in complex algebra. The complex-quaternion transformation is carried out following Eqn. 27 in ref. 69 Q γ pq = Re(γ pq ) +ǐIm(γ pq ) +ǰRe(γ pq) +ǩIm(γ pq ), (11) in which lowercase symbols with (without) bars indicate the Kamers pairing of the original MO basis and the one-particle reduced density matrix (1RDM) is now indicated by γ for consistency of notation with that of ref. 69. This quaternion form can be diagonalized using the quaternion diagonalization routine provided by DIRAC and be back-transformed to complex representation by the routines provided in this module. If the original basis did not possess Kramers symmetry, the diagonalization is directly carried out in complex algebra. The resulting full set of FNOs is then stored on file for analysis and processing by module (C) described below. The module responsible for the quaternion transformation and diagonalization is constructed in such a way that it can also be used with wave function models other than MP2.

C. Selection and recanonization
The third module retrieves FNOs that have occupation numbers above treshold and transforms the original Fock matrix to the truncated virtual space. After recanonization (Eqn. 7) and transformation to AO-basis (Eqn. 10) the final reduced set of FNOs is stored on file. As this step takes virtually no time compared to the other steps in the procedure, this may be easily repeated to test out the effect of varying the treshold value on the generated recanonized orbitals and their energies.  The MP2FNO workflow is schematically represented in Figure 1, and consists of the following steps: 1. A calculation in the full basis set is performed to get a set of occupied and virtual molecular orbitals.
2. The virtual-virtual block of MP2 density matrix is generated. We note that in this step is possible to either consider the full virtual set, or already employ a preliminary screening in which for instance very high-energy virtuals are neglected.
3. Optional, only for restricted starting orbitals: Transform the complex D (2) vv into the quaternion representation, D vv to determine FNOs and occupation numbers.
5. For a given threshold: truncate the FNO space and transform the Fock matrix to recanonize and write the orbitals out in the AO basis.
6. Use the recanonized FNOs in higher level wavefunction models such as coupled cluster.
Here we make use in all calculations of the exact two component Hamiltonian (X2C) 32 , in which we include twoelectron spin-orbit contributions via to the untransformed twoelectron potential via atomic mean-field contributions calculated with the AMFI code [81][82][83] . As part of the supplementary information, we provide in table S2 results with different Hamiltonians (Dirac-Coulomb, X2C, spinfree X2C and Non-relativistic) for expectation values for the HTs molecule. We note the code is equally capable of handling orbitals obtained with the so-called molecular-mean-field approach, in which the transformation to 2C is carried out after a 4C meanfield calculation 39 . In the case of the HCl molecule, we have also employed the non-relativistic Hamiltonian (as specified by the .NONREL keyword), both with contracted and uncontracted bases sets.
The molecular structures employed in all calculations have been taken from the literature: In the case of hydrogen halides, from Huber 84 for HF to HI, from Gomes and Visscher 41 for HAt and from Thierfelder et al. 42   Besides calculating MP2, CCSD and CCSD(T) energies, we have obtained electric dipole moment (EDM), electric quadrupole moment (EQM) and electric field gradient (EFG) for the HX systems, and parity violation energy differences (PV) for H 2 Y 2 systems. These properties have been obtained analytically for CCSD wavefunctions using the implementation described in ref. 69. In the calculations, the size of the complete virtual spinor spaces is 156 (HF), 182 (HCl), 310 (HBr), 382 (HI), 544 (HAt) and 588 (HTs), 314 (H 2 O 2 ), 366 (H 2 S 2 ) and 622 (H 2 Se 2 ). Unless otherwise noted, all electrons were correlated in the calculations. For the smaller systems HF, HCl, HBr, H 2 O 2 and H 2 S 2 we were able to perform calculations on a single node of the laboratory compute cluster in Lille. For HI, HAt, H 2 Se 2 and HTs we employed, respectively, 32, 32, 49 and 64 nodes of the Summit supercomputer. For HI calculations with double-zeta quality basis sets, we were able to perform single-node calculations on the Jean Zay supercomputer.
The data presented in this manuscript is available at the Zenodo repository of reference 86.

V. RESULTS AND DISCUSSION
In this section we discuss the performance of the MP2FNOs approach for the CCSD correlation energy and molecular properties, and the impact of employing different occupation number thresholds on the results. In order to minimize computational cost on our systematic studies, we carried out CCSD(T) calculations only when discussing the effect of FNO and CMO truncation on HCl bond lengths.
For minimizing the bias in the CMO and FNO comparison, we report results in which we employ the closest possible number of CMO to that of FNO, such that for CMO we avoid truncations that would remove close-lying orbitals such as those belonging to a same atomic shell. In practice, this means CMO calculations generally contain a few more orbitals than FNO ones.
In figure S1 in the supplemental information, we present for the EFG of HI the difference between the approach taken here for CMO, and that of considering striclty the same number of CMO and FNO. We can see that arbitrarily truncating the virtual space for CMO can lead to large oscillations in the expectation values, a point that will be addressed in detail in the following. Additionally, from the comparison of Hamiltonians for HTs in table S2 in the supplemental information, we see that our choice of the X2C Hamiltonian is a suitable one for our purposes; while total energies obviously are very different between Hamiltonians, we see that for those taking into account SOC, relative energies (HOMO-LUMO gap) and expectation values are very close to each other. This is in stark contrast with scalar relativistic and non-relativistic calculations, which still show fairly good agreement with 2C/4C approaches for the quadrupole moment and even the HOMO-LUMO gap, but are completely off the mark for the EFG and dipole moment. This is due to the very strong spin-orbit splitting in the valence p-shell of the superheavy tenessine.

A. Correlation Energy
CCSD correlation energies obtained with CMO and FNOs are displayed in Figure 2 for the hydrogen halides. In these, we show values computed for four values of natural occupation number threshold (1.0 × 10 −3 , 1.0 × 10 −4 , 1.0 × 10 −5 , 1.0 × 10 −6 ), as well as the two extrema: a point without any virtual orbitals (corresponding to the Hartree-Fock solution), and another without any truncation of the virtual space. Since the number of virtuals varies for each system, we provide relative measures in terms of the percentage of the virtual space included in the calculations for each point, and likewise for the amount of correlation energy recovered at each point.
From the plots, it is clearly the case that FNOs show a more rapid convergence than CMO across the periodic table, as could be expected from the non-relativistic literature 56,59,87 . For example, to recover only 50% correlation energy already 40% to 50% of the virtual CMO space is required, while only 20% of the FNO space suffices. Also for the more realistic goal of attaining at least 90% of the correlation energy use of FNOs can reduce the required space by 10s of % thereby introducing significant savings in (memory) storage of tensors with virtual indexes, as well as in operation counts for contractions involving virtual indexes in the CCSD equations.
An interesting finding is that the area enclosed by two lines slightly decreases as we go down the periodic table from F to Ts. This means that the FNOs approach recovers less correlation, for a given truncation, as the number of electrons increases. For example, a threshold of 1.0 × 10 −5 (the 4th point in Figure 2) means, for HF and HCl, that more than 60% of the virtual space is included in the calculation, whereas for HTs the same threshold only includes 43%; this means that for HTs virtual orbitals with low occupation numbers are still important for the correlation treatment.
Apart from testing the performance for single point at the potential energy surface, it is also important to verify the performance of FNOs at different geometric structures. To this end, in Figure 3 we present the potential energy curves around equilibrium for HCl at CCSD(T) level, comparing sets of FNOs and CMO virtuals that are truncated to the same number of virtuals (corresponding to about 50% of the complete virtual space). We see that, except for a global energy shift, the FNOs follow the curvature of full virtual space potential energy curve more closely than the CMO curve does.
The better agreement of FNOs relative to CMOs can be quantified by a comparison of spectroscopic constants for the three curves, shown in Table I. There, we report the equilibrium distance (R e ) and vibrational constant (ω e ), calculated with the LEVEL program 88 . Taking the full valence space result as a reference, the truncated orbitals overestimate the R e but underestimate the ω e . However, the error of truncated FNOs is 0.0037Å(R e ) and 17 cm −1 (ω e ), which is only half that of truncated CMOs. As we correlated all electrons, but employed basis sets without specific core correlating functions, one may ask what would be the effect of including such functions. In table S1 of the supplemental information, we provide a comparison of correlation energies and expectation vales between the uncontracted aug-pVTZ and aug-pCVTZ basis sets for the HCl molecule without using virtual space truncation. These results clearly indicate that, although the additional high angular momentum core correlating functions are quite important for increasing the amount of correlation energy for the core and core-valence interactions, the two types of uncontracted basis sets agree well for the computed property values. Another interesting point to note is that, as can be seen from table S3 of the supplemental information, the percentage of correlation energy recovered at each virtual space truncation point is nearly the same for MP2 and CCSD. This indicates that the information obtained from the full virtual space MP2 calculation preceding the truncation can potentially be used to correct the CCSD correlation energy for the effect of truncation. While going beyond the scope of this paper, this point merits further investigation, in which also correction for truncation errors in higher order methods such as CCSD(T) could be investigated.

Electric Dipole and Quadrupole Moments
To asses performance of FNOs for molecular properties we first considered the molecular EDM and EQM because these do well characterize the overall electronic charge distribution within molecules. These properties sample the regions away from the nuclei, as is clear from their operator forms in Eqs. 12 and 13, respectively: Figure 4 shows how the EDM and EQM correlation correction varies with the truncation of the virtual orbital spaces, with dotted and solid lines representing EDM and EQM, respectively. As before FNOs results are plotted with square markers. For these properties, the convergence is nonmonotonic, unlike energy that we considered before. For both CMO and FNO truncation, the use of a high truncation value Figure 2. Convergence of the CCSD correlation energy with respect to the size of the virtual orbital space, for the X2C Hamiltonian. The X axis indicates the fraction of the virtual space retained, while the Y axis gives the fraction of the correlation energy recovered with respect to the value obtained with the untruncated virtual space. and consequently small virtual space (less than 30% virtual orbitals) may lead to a strong overestimation or even a wrong sign of the correlation correction to the molecular property. Especially for CMO truncation, it is almost impossible to estimate the truncation error from a sequence of results for small virtual spaces. For FNO truncation oscillations are much less pronounced than for CMO, which indicates that also for these properties it is advantageous to use FNO truncation rather than CMO truncation. The similar performance for the correlation energy could thereby be used as a guideline. Taking the HF molecule as example, we note that with a threshold of 1.0×10 −4 , the FNOs recover 72% of the correlation contribution to EDM while recovering 85% of the correlation energy. This indicates that the simplification of using a common cutoff for both properties could be a suitable strategy. In addition, we note that the performance of FNOs in EDM also applies to the heavier elements in the studied series. Even for HTs, setting the truncation threshold to 1.0 × 10 −4 (using only 23% of the orbital space) already recovers 91% of the correlation contribution to the EDM.
The EQM shows a similar behaviour as the EDM, albeit with stronger oscillations. The largest oscillations for the EQM occur at the same position as for the EDM, but the amplitude thereof is much larger. This is again most pronounced in the CMO case. Taking HI as an example, retaining 12% of the CMO virtual space overestimates the correlation contribution to EDM by 150%, while the EQM contribution is even seven times too large. As for the EDMs, the oscillations resulting from FNO truncation are smaller and convergence is more smooth, which should make it possible to use a more aggressive cut-off strategy than is possible with CMOs.

Electric field gradient
To also consider the effect of truncation on properties that probe the regions close to atomic nuclei we now turn to the electric field gradients (EFGs) at the halogen nuclei. EFGs couple with nuclear quadrupole moments 89 for nuclei with spin greater or equal to one and are important in the analysis of nuclear magnetic resonance (NMR) experiments 90 . For EFGs both (semicore) correlation and SOC effects are known to be of importance 91 and they form therefore a good second test for the applicability of FNO truncation in relativistic calculations. The EFG is defined as the second derivative of the electric potential V(R) with respect to the nuclear position vector, taken at the nuclear position R N By introducing the EFG tensor operator 92 the EFG can be expressed as the expectation value of a onebody operator which makes evaluating its value similar to computing the EDM and EQM. For linear molecules it suffices to compute only the zzcomponent (with the z-axis chosen along the molecular bond) as the other non-zero parts of the tensor can then be determined by symmetry. Figure 5 plots this zz component of the EFG and shows for both FNO and CMO oscillations upon truncating the orbital spaces. In this case FNO truncation can lead to similar and even more rapid oscillations than with CMO truncation. Recalling that the only difference between computing the EFG and EDM properties is in the operators used, it is of interest to consider their differences. The EFG tensor operator scales as r −3 in contrast to the r 1 and r 2 scaling of the EDM and EQM operators. Density changes in the core region due to correlation are thus magnified by the EFG operator while they are hardly of influence for the EDM and EQM. Such changes may both come from core correlation as well as from correlating the valence electrons (through the tails of the valence orbitals in the core region). As core correlation may be harder to converge than valence correlation, we will separate the two effects by comparing all-electron and valenceonly electrons calculations in the next subsection.

Convergence analysis for HCl and HI
The issue of oscillating convergence can be conveniently analyzed for the HCl molecule as chlorine is large enough to investigate the effect of core correlation, yet small enough to allow for quick calculations. Rather than taking only a few truncation values, we now systematically extend the size of the virtual space by adding individual FNOs and show the effect thereof on the EDM, EQM and EFG of HCl in Figure 6 for both all-electron and valence-only correlation calculations.
For the all-electron case, we see in the EDM and EQM plots for CMO truncation peaks when 13 and 23 orbitals are used. This is probably due to quasi-degeneracies, as the orbital energies of the 12th, 13th, 14th orbital are 0.42418, 0.43657 and 0.43673 Hartree, respectively. This set of virtual orbitals is primarily a diffuse chlorine p-type orbital shell that is split due to the formation of the hydrogen-chlorine bond as well as due to spin-orbit coupling. Taking only one of the two almost degenerate π-type orbitals in the correlation space, gives an unbalanced description and impacts these valence properties significantly. For the EFG this particular set of three CMOs is less important because diffuse orbitals only contribute indirectly to this property. For the EFG one may notice ( Figure 6, lower left panel, blue line) a large oscillation around virtual orbital number 70 in the all-electron calculation. This is due to the more core-like p-orbitals 70, 71, and 72 with energies of respectively 35.63, 35.73 and 36.10 Hartree, that are impacting the correlation of core electrons. Restricting the correlation treatment to valence electrons only strongly reduces the effect of these virtuals ( Figure 6, lower right panel).
The FNO curves for the EDM and EQM show only an initial oscillation that is followed by rather smooth convergence. This initial oscillation can also be traced to quasi-degenerate orbitals (but now in terms of occupation numbers). For the FNOs, the occupation number of the 2nd, 3rd, 4th, 5th orbital are 0.008692, 0.008676, 0.008508, 0.008491, respectively, and like with near-degenerate CMOs it appears recommendable to either include all or none of a degenerate set. For the EFG, the situtation is unfortunately more complicated, also with FNOs.
In the all-electron calculation values only stabilize after inclusion of about 50 orbitals, while in the valence-only calculation stable convergence is reached after addition of about 20 FNOs. Compared to the CMO truncation scheme, the advantage of FNO truncation appears to be absence of "late" oscillations (at low threshold values) that could cause artefacts in the CMO truncation schemes. Such oscillations are typical for EFGs, for which an indivdual orbital may provide a significant contribution but where the contribution of full, spherically symmetric shells of orbitals add up to zero due to symmetry. This is well-known for CMOs, but also holds for FNOs. An example are the EFG integrals of the 13th, 14th, 15th FNOs which evaluate to -78, 44 and 26 au. These large values get multiplied by very similar occupation numbers (resp. 7.9 10 −4 , 7.0 10 −4 , 7.0 10 −4 ) so that their contributions to the total EFG nearly cancel.
A similar argument can be used to rationalize the difference between the valence and all-electron calculations. In the correlation of core orbitals one effectively changes the relative occupation of the three 2p orbitals in the chlorine core by making their occupations slightly smaller than 2. Imbalances in the correlation get thereby for the EFG multiplied by large integral values (about 80 au for the 2p 3/2 orbitals) which creates the large oscillation seen in these calculations. For EDMs and EQMs the integral values are much more alike and smaller in magnitude, leading to the observed smoother convergence.
To check if such oscillation problems are particular for the relativistic domain, or can also be found in the non-relativistic case, we have carried out the same analysis as above by employing a non-relativistic Hamiltonian, employing both contracted and uncontracted (valence) basis sets and correlating both valence and all electrons. The results are shown in Fig. 7.
From these, it can be seen that for the electric dipole and quadrupole moments, the non-relativistic calculations basically show the same oscillations as for X2C, and results for contracted and uncontracted basis sets are nearly indistinguishable from each other. For the EFG there is also no noticeable difference for the convergence patterns between X2C and non-relativistic results for uncontracted basis sets.
On the other hand, we observe significant differences between EFG calculations employing contracted and uncontracted basis sets. The first one is that, already at the Hartree-Fock level, there is a significant difference (nearly 10%) on the absolute value of the EFG, with contracted values underestimating uncontracted ones. Second, we note that with con- This can be understood as a manifestation of the degree of atomic symmetry that is imposed by the contraction for the different orbital shells, particularly for the occupied core orbitals. In the contracted sets these are forced to maintain their atomic-like nature during the molecular calculation because electron correlation cannot as easily deform the orbitals as is possible when the basis set is uncontracted. This is underscored by the nearly identical behavior of the contracted all electron and valence calculations, as in the latter core orbitals are kept frozen.
To check if the convergence behavior mentioned above is also found for a heavier system, we carried out similar calculations for HI, but using an uncontracted double zeta basis set. It can be seen from Fig 8 that although HI has many more electrons, it still shows a very similar convergence as observed for HCl. For instance, to get a converged EFG value of HI, in the valence only computation, one just needs to correlate 20 FNOs.
Looking at these tests, we may conclude that for properties like EDM and EQM, which are dominated by valence electron contributions, one can safely use FNOs and their occupation numbers to cut off the orbital space, both in valence-only and all-electron correlation calculations. For sensitive properties for which core electrons may give large and nearly cancelling contributions, it is also with FNOs numerically more stable to correlate only the valence electrons.

C. Parity violation
Detection of parity violation (PV) effects, associated with weak force in atoms and molecules is an active field of research 13,93,94 . While this property can also be computed with perturbation theory starting from nonrelativistic theory 95 , it is advantageous to use a relativistic quantum chemistry framework because the PV energy can then be formulated as an expectation value of an effective one-body operator 85 : and  θ ω is Weinberg mixing angle, that is set to 0.2319 for sin 2 θ ω . ρ A and γ 5 i are normalized nucleon density and 4-dimensional chirality operator, respectively.
Sunaga et al. 13 already investigated use of CMO trunca-tion at CCSD level with two different threshold values and showed that truncation is well possible for this property. In current work, we also employ FNO truncation and test more thresholds for CMO truncation and its effect on the PV value.
The results of the two truncation schemes for the H 2 Y 2 (Y=O, S, Se) molecules are displayed in Figure 9. We find that both FNO and CMO truncation leads to strong oscilla- tions like observed in the EFG case. Again, these are most pronounced for aggressive truncation in which more than 50% of the virtual orbital space is removed. Looking more closely at these systems we find that orbital energies (for CMO) or occupation numbers (for FNO) show a number of near-degeneracies. These orbitals give individually large but partially cancelling contributions. We have not investigated valence-only correlation explicitly for this property as core contributions are probably more important than for the EFG and it is more difficult to define a representative small test case. For our current tests we note that the convergence behaviour is consistent with results of reference 13, in which rather conservative CMO truncation tresholds of 100 and 500 Hartree were used for H 2 Se 2 . With both CMO and FNO truncation a 50 % reduction appears realistic, leading to a significant reduction of the computational effort. For instance, in Figure 9. Convergence of the CCSD PV energy with respect to the size of the virtual orbital space, for the X2C Hamiltonian. The X axis indicates the fraction of the virtual space retained, while the Y axis gives the fraction of the expectation value recovered with respect to the value obtained with the untruncated virtual space. the current test the H 2 Se 2 calculation on 49 nodes of Summit with the full orbital space took 100 minutes while the 50 % reduced calculation too only 20 minutes.
Nevertheless, for the particular case of PV energy evaluation, FNO truncation does not appear to give a significant advantage over the simpler CMO truncation scheme. Apart from the sensitivity of PV operator on core correlation already observed in the EFG evaluation, another possible reason is our use of the unrelaxed MP2 density matrix. Shee, Visscher, and Saue 69 found for H 2 S 2 that the contribution of orbital relaxation to the density matrix is significant when used to evaluate this property. It would be interesting to investigate how these contributions impact the generated FNOs are thereby the convergence of the FNO truncation scheme.

VI. CONCLUSIONS
In this work we describe the formulation and implementation of the MP2 frozen natural orbitals (MP2FNOs) method for relativistic electronic structure calculations, which is an appealing approach for truncating the large virtual orbital spaces typically associated with relativistic calculations, while retaining the high accuracy expected of wavefunction-based approaches such as coupled cluster.
This implementation was carried out in the massively parallel coupled cluster module of the DIRAC program, with the help of a framework to manipulate 1RDMs obtained from correlated wavefunction calculations. A particularity of our code is its ability to generate a set of canonical occupied and truncated virtual orbitals as well as natural orbitals in the atomic basis, so that these can be conveniently employed in post-Hartree-Fock calculations, and also used for analysis.
We employed our code to investigate the performance of MP2FNOs for the calculation of correlation energies and ground-state first-order properties such as dipole and quadrupole moments, electric field gradient at the nuclei, and parity violation energy shifts. As model systems, we considered species containing elements from the first row up to and including the superheavy element tennessine.
We have found that though MP2FNOs are always capable of recovering more correlation energy than their canonical counterparts for a given truncation, this advantage is slightly diminished for the heaviest systems considered, containing Astatine and Tennessine. In spite of that, it is generally found that MP2FNOs can reduce virtual spaces to about 50% of their original size without significant errors in the energy.
Among the properties considered, a truncation of about 50% of the original virtual space has also been shown to provide values that are nearly converged to the results without truncation. That said, different properties exhibit very different convergence behavior; for valence properties such as the electric dipole moment MP2FNOs show a fairly smooth convergence to the reference values, whereas for the electric field gradient and in particular the parity violation energy shifts, for which regions close to the nuclei are important (and thus higher-lying virtuals are more important in the correlation treatment), significant variations on the calculated values are found for virtual spaces smaller than 50% of the original virtual space.
We have found that for EFGs, correlating only valence electrons provides a pragmatic solution that recovers the smooth convergence seen for valence properties by removing difficult to converge, individually large but in total nearly cancelling, contributions from the core electrons. Whether such a strategy can also work for the PV energy shifts remains to be investigated.
As a final point and perspective we note that improvement of the efficiency of the current scheme is well possible by implementing efficient approximate schemes like Cholesky decomposition and Laplace transforms to generate the MP2 1RDM. Work along these lines is in progress and should enable treatment of larger systems in the future. Figure S1: Convergence of the CCSD EFG at the I atom with respect to the size of the virtual orbital space, for the X2C Hamiltonian, comparing FNO and CMO, the latter being truncated at exactly the same number of FNO and approximately respecting atomic shell boundaries as shown in the article body. The X axis indicates the fraction of the virtual space retained, while the Y axis gives the fraction of the correlation energy recovered with respect to the value obtained with the untruncated virtual space.